Fitting a yield curve to bond data using Matlab's Financial Instruments Toolbox

Ken Nyholm, August 2014

Financial market analysis of fixed income markets typically rely on the availability of yield curve data. That is, at the outset, yield observations for the relevant market segments are directly observable at the desired maturities. For example, it is possible to download yield curve data from Bloomberg, the FRED data base, and ECB's statistical data warehouse. Such data comprise the date of observation, the maturities at which yields are observed, and then the actual yield curve data points. To illustrate, Figure 1 shows the Spanish Government zero coupon yield curve for Monday 11 August 2014 (close of business), as displayed by Bloomberg, for maturities ranging from 3 months to 30 years.

y_ES = [0.075; 0.133; 0.334; 0.628; 0.884; 1.216; 1.532; ...
        1.841; 2.178; 2.474; 2.711; 3.499; 3.716; 4.341];

tau_ = [ 3; 12; 24; 36; 48; 60; 72; 84; 96; 108; 120; 180; 240; 360]./12;

plot(tau_,y_ES,'-ob'), ... %title('Figure 1: Spain - Gov Curve on 11 AUG 2014'), ...
                     xlabel('Maturity in years'), ylabel('Yield in pct')

%saveas(gcf,'Spain_Gov.eps', 'psc2')

Since zero coupon instruments typically do not trade trade throughout the maturity spectrum, it is necessary to construct the zero coupon curve from underlying (traded) coupon bonds. Bloomberg and other providers of yield curve data perform such calculations by default; but it is naturally relevant to be able to reproduce such calculations. Firstly, as a means to double check and to gain an understanding of how such processes work, and secondly, and more importantly, by mastering the curve extraction methodology it is possible to create proprietary curves for yield curve segments that are not available from the licensed curve contributors, or as a means to create curves for combinations of curve segments that may be needed for ad-hoc analyses.

In principle, the required calculations are very simple, and of course, relies on the discounting of future cashflows. Define a universe of traded bonds that forms the foundation for the zero coupon yield curve calculations. These bonds should be activily traded and representative for the market segments under investigation. Let there by $j$ bonds in the universe, and denote by $p$ the vector that collects the $j$ observed prices for these bonds, and by $C$ the matrix that collects their cashflows. Since the price of a bond is equal to the sum of the discounted future payments, the following holds: $p=C \cdot d$, where $d$ is the vector of discount factors. Zero coupon yields $y$ are related to discount factors as $y_j=d^{-1/j}-1$. Under the assumption that $C$ is invertible, the discount factors can then be found as: $d=C^{-1} \cdot p$, which then subsequiently can be converted in to zero coupn yield. The following examples illustrates this principle for the first four maturities, using constructed/hypothetical bond price data:

j  = (1:1:4)';
p  = [ 99.87; 105.82; 104.35; 112.75 ];

C  = [ 100.00    0.00    0.00    0.00 ;
         3.25  103.25    0.00    0.00 ;
         2.10    2.10  102.10    0.00 ;
         4.10    4.10    4.10  104.10 ];

d  = C^-1*p;

yz = ( d.^(-1./j)-1 ).*100;

plot(tau_,y_ES,'-ob', j, yz,'-*r'), ... % title('Figure 1: Spain - Gov Curve on 11 AUG 2014'), ...
                     xlabel('Maturity in years'), ylabel('Yield in pct')

%saveas(gcf,'ZeroYC_stylised.eps', 'psc2')

Figure 2 plots the calculated (example) zero rates together with the observed curve. The data used for this illustration were constructed on the basis of the characteristics of the bonds that actually underpin Bloombergs zero curve calculations. Using the coupons of these bonds, the prices were adjusted slightly to provide zero yields that are close to those seen in Figure 1.

While the above example seems easy to implement and provides results that are close to the observed Bloomberg zero yields, it is far from being implementable to real-life data. Bonds that trade in the markets are not necessary well spaced through out the maturity spectrum. Government bonds are used by the sovereign treasury department to fund the government debt. In doing so, issuances are timed and maturities are structured such that the debt is best serviced. Therefore, it cannot be guaranteed that liquidly traded bonds span all maturities as desired by the financial econometrician. In case that gaps emerge in the maturity spectrum, it is necessary to apply some form of interpolation methodology, to close the gaps. In addition, the maturity profile of the outstanding bonds change as time passes: a bond with a maturity of one year today, becomes a bond with a residual maturity of eleven months, one month from now. Another complicating issue relates to the 'book-keeping' of nitty-gritty bond specific details and how these affect the discount factors. Bonds can be issued using different 'day-count' conventions, i.e. different ways by which interests accrue according to agreed calendar-function schemes, e.g. actual/actual, 30/360, actual/360 and so on. Furthermore, coupons can be paid annually, semi-annually, quarterly, etc.

Matlab's Financial Instruments Toolbox provides a means to estimate zero coupon curves while ensuring correct 'book-keeping' of the bond specific details, as well as three readily available curve interpolation techniques. The R2012b version of this toolbox allows for interpolation to be performed using the Nelson-Siegel, Svensson (-Soderlind), and Smoothing Splines. As shown below, user specific functions can be implemented as well. Once the parameters of the chosen model specification has been estimated, Matlab will provide the zero coupon curve, the instantaneous forward curve, the par-yield curve, and discount factors.

In order to perform the calculations, the 'IRFunctionCurve' object has to be populated with data and parameters. The parameters are entered in several rounds. First, the ones that characterise the precision of the estimated parameters. These are entered via the Matlab structure 'optimset':

optOptions_  = optimset('TolFun', 1e-12, 'TolX', 1e-12, 'MaxFunEvals', 1e12, ...
                        'MaxIter', 1e12, 'Display', 'final' );

Here it is specified that convergence of the objective function is achieved when the it changes by less than 1e-12 from one iteration to the next ('TolFun'), or when the parameters of the objective function changes by less than 1e-12 from one iteration to the next ('TolX'). The number of calculations that are allowed before the optimisation algorithm terminates without finding am optimal solution is specified by 'MaxFunEvals' and 'MaxIter'. In other words, the algorithm will terminate after 1e12 calculation atempts unless the an optimal solution is found before, where the optimality criterion is specified by the values entered for the variables 'TolFun' and 'Tolx'.

Options that govern the interest rate curve fit are provided through the 'IRFitOptions' structure. The choices made for the optimisation algorithm via 'optimset' must also be passed to Matlab via this structure.

b_0  = [ 4.00  -3.00  0.40  -0.50  3.50  1.50 ];
lb_  = [ 0.00  -inf  -inf  -inf   0.00  0.00 ];
ub_  = [ 25.0   inf   inf   inf    inf   inf ];

fitOptions_ss  = IRFitOptions(b_0, 'FitType', 'durationweightedprice', ...
                                   'LowerBound', lb_, 'UpperBound', ub_, ...
                                   'OptOptions', optOptions_ );

Starting values for the optimisation process is provided via the variable 'b_0' and lower and upper bounds for the parameters are fed to Matlab via the variables lb_ and ub_, respectively. It is important to note that different targets can be specified for the employed optimisation algorithm. This is done via variable 'FitType', which can take on the values 'yield', 'price', 'DurationWeightedPrice', depending on which criterion the user wants to employ. The parameters of the objective function is found by minimising the sum of squared residuals between an observed quantity and the corresponding model quantity. 'FitType' allows the user to decide whether it is the residuals calculated from yields, prices or duration weighted prices that are minimised.

Lastly, the model must be estimated. This is done through an 'IRFunctionCurve' call to a specific fitting function: '.fitSvensson', '.fitNelsonSiegel', or 'fitSmoothingSpline'.

The data used in the example below are again taken from Bloomberg. In fact, the bonds included are the ones that underpin Bloomberg's calculations of the Spanish zero coupon curve shown in Figure 1.

p_   = [99.989; 99.974; 99.937; 99.897; 105.128; 104.15; ...
        112.545; 107.49; 114.015; 123.958; 126.74; 116.093; ...
        100.863; 121.748; 131.85; 109.628; 122.798 ];

cpn_ = [ 0.00; 0.00; 0.00; 0.00; 3.25; 2.10; 4.10; 2.75; ...
         4.00; 5.50; 5.85; 4.40; 2.75; 5.15; 5.75; 4.20; 5.15]./100;

maturity_ = datenum([ '17/10/2014'; '23/01/2015'; '10/04/2015'; '17/07/2015'; ...
             '30/04/2016'; '30/04/2017'; '30/07/2018'; '30/04/2019'; ...
             '30/04/2020'; '30/04/2021'; '31/01/2022'; '31/10/2023'; ...
             '31/10/2024'; '31/10/2028'; '30/07/2032'; '31/01/2037'; ...
             '31/10/2044' ], 'dd/mm/yyyy');

settle_      = repmat(datenum(['11/08/2014'],'dd/mm/yyyy'),length(p_),1);

Instruments_ = [ settle_ maturity_ p_ cpn_ ];

cpn_freq     = [ 0; 0; 0; 0; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1; 1 ];

day_cnt      = [ 0; 0; 0; 0; 8; 8; 8; 8; 8; 8; 8; 8; 8; 8; 8; 8; 8 ];

The variable 'p_' collects the observed prices, 'cpn_' contains the coupon rates paid by the bonds, 'maturity_' holds the dates at which the bonds mature and 'settle_' is the settlement date.

est_Svensson = IRFunctionCurve.fitSvensson('Zero', settle_(1,1), Instruments_, ...
                                                 'Compounding', -1, ...
                                                 'Basis', 8, ...
                                                 'InstrumentPeriod', cpn_freq, ...
                                                 'InstrumentBasis', day_cnt, ...
                                                 'IRFitOptions', fitOptions_ss);
Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the selected value of the function tolerance.

In addition to the data, it is also necessary to specify that the 'Zero' curve should be estimated, which compounding and day count convention the estimated curve should follow: ('Compounding', -1) and ('Basis', 8), respectively. This means that continous compounding is used and the day count convention 'actual/actual'.\footnote{See the Matlab help documentation for a list of day count conventions and their numerical equivalent.}. Similarly, for the bond universe that provides the foundation for the calculations, the coupon frequency should be provided for each bond ('InstrumentPeriod', cpn_freq,), and the day count convention for each bond ('InstrumentBasis', day_cnt). For the bond universe used in the current example, the first four bonds are qouted using the convention of 'actual/360', which in Matlab terminology is equivalent to '0'; and the remaining bonds are quoted using 'actual/actual', which is equivalent to a nummerical value in Matlab of of 8.

After the model parameters have been estimated the Svensson zero curve is found by extracting the vector from the structured variable 'est_Svensson'. This is done here:

tau_m     = settle_(1,1) + 365*tau_;
y_zero_SS = est_Svensson.getZeroRates(tau_m).*100;

And, the resulting Svensson-Soderlind zero coupon curve is plotted along side the zero curve reported by Bloomberg.

plot(tau_,y_ES,'-ob'), ...% title('Figure 3: Spain - Gov Curve on 11 AUG 2014'), ...
                     xlabel('Maturity in years'), ylabel('Yield in pct')
hold on
plot(tau_,y_zero_SS,'-*r'), legend('Bloomberg','Svensson','Location','SouthEast')

%saveas(gcf,'SvenssonZC.eps', 'psc2')

There are two additional built-in curve fitting models available in Matlab's toolbox: Nelson-Siegel and Smoothing Splie. These models are estimated below following the principles outlined above and using the same data as above.

The Nelson-Siegel model requires less parameters than the Svensson model. While Svensson relies on a parametrisation of the yield curve using a Level, Slope, Curvature1 and Curvature2 factors, and estimates time-decay parameters for the two curveature factors, the Nelson-Siegel model relies on a more parsimonious parametrisation comprising Level, Slope and just one Curvature factor. For this reason the input parameters need to be adjusted.

b_1   = [ 4.00  -3.00  0.40  1.00 ];
lb_1  = [ 0.00  -inf   -inf  0.00 ];
ub_1  = [ 25.0   inf    inf  inf  ];

fitOptions_ns  = IRFitOptions(b_1, 'FitType', 'durationweightedprice', ...
                                   'LowerBound', lb_1, 'UpperBound', ub_1, ...
                                   'OptOptions', optOptions_ );

est_NS = IRFunctionCurve.fitNelsonSiegel('Zero', settle_(1,1), Instruments_, ...
                                                 'Compounding', -1, ...
                                                 'Basis', 8, ...
                                                 'InstrumentPeriod', cpn_freq, ...
                                                 'InstrumentBasis', day_cnt, ...
                                                 'IRFitOptions', fitOptions_ns);
y_zero_NS = est_NS.getZeroRates(tau_m).*100;
Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the selected value of the function tolerance.

To estimate the Smoothing Spline curve, the input parameters have to be adjusted accordingly. The 'penalty_fx' is a penalty funcition that is used to ensure a certain smoothness of the forward curve. The parameter values for L, S and mu are based on the example shown in the Matlab help documentation.

L = 9.2; S = -1; mu = 1;
penalty_fx = @(t) exp(L - (L-S)*exp(-t/mu));

est_SmSpline = IRFunctionCurve.fitSmoothingSpline('Zero', settle_(1,1), Instruments_, ...
                                                 penalty_fx, ...
                                                 'Compounding', -1, ...
                                                 'Basis', 8, ...
                                                 'InstrumentPeriod', cpn_freq, ...
                                                 'InstrumentBasis', day_cnt );

y_zero_SmSpline = est_SmSpline.getZeroRates(tau_m).*100;

For comparison, Figure 4 plots the estimated Nelson-Siegel and Smoothing Spline curves together with the Svensson and Bloomberg curves.

subplot(2,2,1), plot(tau_,y_ES,':ob','LineWidth',2), title('Gov Curve on 11 AUG 2014 - Bloomberg'), ...
                     xlabel('Maturity in years'), ylabel('Yield in pct')

subplot(2,2,2), plot(tau_,y_ES,':ob','LineWidth',2), title('And, Svensson'), ...
hold on
plot(tau_,y_zero_SS,':*r','LineWidth',2), ...
                      xlabel('Maturity in years'), ylabel('Yield in pct')

subplot(2,2,3), plot(tau_,y_ES,':ob','LineWidth',2), title('And, Nelson-Siegel')
hold on
plot(tau_,y_zero_NS,':*r','LineWidth',2), ...
                      xlabel('Maturity in years'), ylabel('Yield in pct')

subplot(2,2,4), plot(tau_,y_ES,':ob','LineWidth',2), title('And, Smoothing Spline')
hold on
plot(tau_,y_zero_SmSpline,':*r','LineWidth',2), ...
                      xlabel('Maturity in years'), ylabel('Yield in pct')

%saveas(gcf,'AllZC.eps', 'psc2')

Now that all three models are estimated, it is easy to also compare their instantaneous forward rate curves. First, the forward rates are extracted from the structured variables that hold the model results, and then the curves are plotted in Figure 5.

forward_SS       = est_Svensson.getForwardRates(tau_m).*100;
forward_NS       = est_NS.getForwardRates(tau_m).*100;
forward_SmSpline = est_SmSpline.getForwardRates(tau_m).*100;

plot(tau_,[forward_SS forward_NS forward_SmSpline]), ...% title('Instantenous Forward Rate Curves'), ...
                  legend('Svensson','Nelson-Siegel','Smoothing Spline','Location','SouthEast'), ...
                  xlabel('Maturity in years'), ylabel('Yield in pct')

%saveas(gcf,'AllForwards.eps', 'psc2')

Fitting curves to observed yields

Yield curve data is often directly available, and it may therefore not be necessary to apply estimation techniques to the bond data (prices, coupons, and maturity dates) to uncover the curves, as it was done in the above section. Instead, it is possible to work directly with the observed yields to fit models that allow for the calculation of forward rate curves, and to facilitate interpolation and extrapolation to maturities that are not covered by the observed data.

Two types of yield curve data are used in the example. One is the zero coupon data for Spanish Government curve observed on 11 August 2014, which is also used in the above example. The second data set comprises observations on the EONIA 3 months swap curve (Bloomberg code: EUSWE[X] CMPN Curncy). These data are shown in Figure 6.

y_ES     = [0.075; 0.133; 0.334; 0.628; 0.884; 1.216; 1.532; ...
            1.841; 2.178; 2.474; 2.711; 3.499; 3.716; 4.341];

y_Swap   = [ 0.060; 0.050; 0.058; 0.117; 0.209; 0.340; 0.496; ...
             0.663; 0.826; 0.977; 1.112; 1.603; 1.828; 1.962 ];

tau_     = [ 3; 12; 24; 36; 48; 60; 72; 84; 96; 108; 120; 180; 240; 360]./12;
nTau     = length(tau_);

hold on
plot(tau_,y_Swap,'-or'), ... % title('Figure 6: Yield Curves Observed on 11 AUG 2014'), ...
                     xlabel('Maturity in years'), ylabel('Yield in pct'), ...
                     legend('Spanish Gov zc', 'EONIA Swap', 'Location', 'SouthEast')
%saveas(gcf,'SpainEonia.eps', 'psc2')

To fit the available models to the Spanish zero coupon curve data it is first necessary to construct the needed input. The settlement date is set equal to the date at which the curve is observed, and the maturities at which the curve points are observed, contained in the vector 'tau_', they need to be converted into corresponding 'dates'. This is done by defining the vector 'maturity_m', which takes the business day closest to the settlement date plus the appropriate numbers of days that correspond to the maturities contained in 'tau_'. It is assumed that there are 365 days in a year. Since zero coupon data is used, the coupon input variable 'cpn_m' is set equal to a vector of zeros, of appropriate length. Finally, the prices of the zero coupon 'instruments' are found as the discount factors for the observed yields. This is done using the built-in Matlab function 'zero2disc'.\footnote{The function 'zero2disc' requires two additional inputs: the type of compounding used (-1 corresponds to continous compounding) and the day count convention (8 corresponds to actual/actual).}

settle_m     = repmat(datenum(['11/08/2014'],'dd/mm/yyyy'), nTau,1);
maturity_m   = busdate( settle_m(1,1) + ceil(tau_.*365) );
cpn_m        = zeros(nTau,1);
price_m      = zero2disc(y_ES./100, maturity_m, settle_m, -1, 8).*100;

Instruments_ = [settle_m maturity_m price_m cpn_m ];

The obtained prices are plotted in Figure 7.

plot(tau_, price_m,'-ob'), ... % title('Figure 7: Zero Coupon Prices'), ...
                     xlabel('Maturity in years'), ylabel('Price')
%saveas(gcf,'SpainPrices.eps', 'psc2')

Once starting values ('b_S') and possible lower and upper bound parameter constraints (lb_S and ub_S) have been decided on, the chosen model can be estimated. Here, the Svensson model is used to fit the observed Spanish zero coupon curve. Results are shown in Figure 8.

b_S   = [ 4.00  -3.00  0.40  -0.50  3.50  1.50 ];
lb_S  = [ 0.00  -inf  -inf  -inf   0.00  0.00 ];
ub_S  = [ 25.0   inf   inf   inf    inf   inf ];

optOptions_  = optimset('TolFun', 1e-8, 'TolX', 1e-8, 'MaxFunEvals', 1e12, ...
                        'MaxIter', 1e12, 'Display', 'none' );

fitOptions_ss  = IRFitOptions(b_S, 'FitType', 'durationweightedprice', ...
                                   'LowerBound', lb_S, 'UpperBound', ub_S, ...
                                   'OptOptions', optOptions_ );

est_SvenssonES = IRFunctionCurve.fitSvensson('Zero', settle_(1,1), Instruments_, ...
                                                 'Compounding', -1, ...
                                                 'Basis', 8, ...
                                                 'IRFitOptions', fitOptions_ss);

y_ES_SS   = est_SvenssonES.getZeroRates(maturity_m).*100;
fwd_ES_SS = est_SvenssonES.getForwardRates(maturity_m).*100;

hold on
hold on
plot(tau_,fwd_ES_SS,':*r'), ...% title('Figure 8: Observed and Fitted Yield and Forward Curves'), ...
                     xlabel('Maturity in years'), ylabel('Yield in pct'), ...
                     legend('Spanish Gov zc', 'Svensson Curve Fit','Svensson Forward', 'Location', 'SouthEast')
%saveas(gcf,'SpainSvensson.eps', 'psc2')

The second data set used is the EONIA swap curve, observed on 11 August 2014. These data are used to illustrate how a zero coupon curve can be extracted from swap data. It is important to emphasise that the exercise does not aim to fit the observed data, rather the purpose is to extract the zero coupon curve from the data. Swap rates are not zero coupon rates, but more similar to coupon bonds. In particular, a swap is an arrangement between two parties to exchange a stream of future floating payment for a stream of future fixed payments. In the case of the EONIA swap the floating payment is constituted by the future EONIA rate, and the fixed payment is equal to the swap rate. In general, the value of a swap contract is equal to zero at the time it is initiated. This means that the observed swap curve, at a given point in time, in some sense represents market equilibrium rates, at that point in time, at which traders are willing to exchange fixed payments for future floating (EONIA) payments/rates. The market equilibrium consideration implies that the swap rate can be treated as a coupon payment on a newly issued par bond, i.e. a bond which price is 100. For this reason, the inputs to the swap rate estimation are as follows:

settle_swap   = repmat(datenum(['11/08/2014'],'dd/mm/yyyy'), nTau,1);
maturity_swap = busdate( settle_m(1,1) + ceil(tau_.*365) );
cpn_swap      = y_Swap./100;
price_swap    = 100.*ones(nTau,1);

Instruments_swap  = [settle_swap maturity_swap price_swap cpn_swap ];

est_Svensson_swap = IRFunctionCurve.fitSvensson('Zero', settle_(1,1), Instruments_swap, ...
                                                 'Compounding', -1, ...
                                                 'Basis', 8, ...
                                                 'IRFitOptions', fitOptions_ss);

y_swap_SS      = est_Svensson_swap.getZeroRates(maturity_m).*100;
fwd_swap_SS    = est_Svensson_swap.getForwardRates(maturity_m).*100;

hold on
hold on
plot(tau_,fwd_swap_SS,':*r'), ... %title('Figure 8: Observed and Fitted Yield and Forward Curves'), ...
                     xlabel('Maturity in years'), ylabel('Yield in pct'), ...
                     legend('EONIA Swap', 'Svensson Zero Curve','Svensson Forward', 'Location', 'SouthEast')
%saveas(gcf,'fitEONIA.eps', 'psc2')


The objective of this note is to provide a brief and hands-on introduction to fitting of single dated (as opposed to panel data) yield curves using Matlab's Financial Instruments Toolbox. The toolbox contains three pre-programmed models, namely (i) Nelson-Siegel, (ii) Svensson (-Sodelind), and (iii) Smoothing Spline, and it facilitates the integration of user defined models, althogh this is not covered by the note. Using the three standard models, it is demonstrated how zero and forward curves can be extracted, on the basis of traded bonds, i.e. how information on bonds prices, coupons, coupon frequencies, maturities, day count conventions and compounding type, can be translated into zeros and forwards. It is also shown how to perform estimations using swap and zero coupon data.