Spot and forward curve from interest rate swaps¶

clear all, clc

[ndata, text, alldata] = xlsread('IRS_plain_vanilla.xlsx');

time = ndata(:, 1);
z = ndata(:, 2);   % z = swap rates
z = z / 100;      % note that it is necessary to divide z by 100, because the rates in IRS_plain_vanilla.xls are already in percentage.

% Point 1
    [v_spot] = F_es59(z);

% Point 2
    % see function F_es60.m

% Point 3
    % see function F_forward

    step = 0.5;                 % time in years (one point every six months)
    time_int = step : step : time(end);

    [v_spot_int] = F_es60(v_spot, time, time_int);

    spot_r_int = (1 ./ v_spot_int).^(1 ./ time_int) - 1;    % spot rates (converts the interpolated prices into annual rates for the plot)
    [fwd_p, fwd_r] = F_es55(v_spot_int, time_int);       % forward prices and rates from the semi-annual spot price curve

% Point 4
   % Graph
    figure(1)
    plot(time_int, spot_r_int, 'b', 'LineWidth', 2)
    hold on
    plot(time_int, fwd_r, 'r', 'LineWidth', 2)
    hold off
    xlabel('Time')
    ylabel('Rate')
    title('Spot rate and forward rate')
    legend('Spot Rates', 'Forward Rates', 'Location', 'best')

% Point 5
    % the wrong procedure on purpose: linearly interpolate the market swap
    % rates to get the semi-annual rates (instead of extracting pure rates first)
    [z_int] = interp1(time, z, time_int, 'linear', 'extrap');

    [v_spot_int] = F_es59(z_int);

    spot_r_int2 = (1 ./ v_spot_int).^(1 ./ time_int') - 1;  % spot rates

    [fwd_p, fwd_r2] = F_es55(v_spot_int', time_int);

    % Graph
    figure(2)
    plot(time_int, spot_r_int2, 'b', 'LineWidth', 2)
    hold on
    plot(time_int, fwd_r2, 'r', 'LineWidth', 2)
    hold off
    xlabel('Time')
    ylabel('Rate')
    title('Spot rate and forward rate')
    legend('Spot Rates', 'Forward Rates', 'Location', 'best')

% CONCLUSION: COMPARING THE TWO CHARTS
% The comparison shows why raw swap rates are never interpolated. The first
% chart shows smooth, realistic curves because the interpolation is applied
% to pure zero-coupon yields, extracted beforehand by bootstrapping. The
% second shows strongly unstable (zig-zag) forward rates with out-of-scale
% values: linearly interpolating swap rates introduces small coupon-related
% distortions, which forward rates (behaving like derivatives) amplify. The
% rule is to always extract pure rates before interpolating.

function [v_spot] = F_es59(z)

% output : v_spot := Zero Coupon Bond prices (discount factors)
% input  : z      := Swap Rate vector

z = reshape(z, length(z), 1);                        % force column vector

M = repmat(z, 1, length(z));
A = tril(M);
A = A + eye(length(z));

v_spot = A \ ones(length(z), 1);                     % solve for discount factors

end

function [ v_spot_int ] = F_es60( v_spot, time, time_int)

% inputs:    v_spot   : price of a ZCB
%            time     : the original maturities
%            time_int : the new maturities at which the rates are wanted
% output:    v_spot_int : vector of ZCB prices on the new time grid time_int

yield = -log(v_spot) ./ time;                                     % yield to maturity: converts each ZCB price into its continuous spot rate
                                                                  % rates interpolate well, prices do not

yield_int = interp1(time, yield, time_int, 'linear');             % linear interpolation of the known points (time, yield) at the intermediate maturities

% Now there are spot rates for every maturity (including the semi-annual ones),
% but pricing needs discount factors again. The inverse formula:

v_spot_int = exp(-yield_int .* time_int);                         % discount factors at the intermediate points

end

function [fwd_p, fwd_r] = F_es55(v, t)

% inputs:   v := ZCB prices
%           t := the schedule

% outputs:  fwd_p := forward prices
%           fwd_r := forward interest rates

% v and t must have the same length!

n = length(v);

% Initialization
fwd_p = zeros(1, n); % preallocation
fwd_r = zeros(1, n); % preallocation

% Calculate the first element of fwd price and rate vectors
fwd_p(1) = v(1);
fwd_r(1) = (1 / fwd_p(1))^(1 / t(1)) - 1;

for k = 2:n
    fwd_p(k) = v(k) / v(k-1);
    fwd_r(k) = (1 / fwd_p(k))^(1 / (t(k) - t(k-1))) - 1;
end

end