Yield Curves with Mathematica 6.0

Jan Hurt

Charles University of Prague, Czech Republic

Abstract

This contribution is concerned with analysis of the profile of interest rates (IR's) over time, also known as term structure of IR's. The graphical representation of the dependence of IR's on time to maturity is the yield curve (YC). A common problem for the financial analysts is either to construct or estimate YC from market data. There is a variety of numerical and statistical techniques for doing it. Some parametric and nonparametric methods of fitting the YC will be discussed and a comparison with results obtained by interpolation will be given. Among other criteria, the integrated mean squared error will be used and an economic interpretation of it will be given. An illustration on real and simulated data will be presented.

Introduction

A criticism should be taken about the approach to yield curves by financial analysts. They usually insist on exact values in quoted data, which may, after an interpolation, be dangerous in the intervals between  quoted data. The author of this contribution is convinced that smoothing would eventually give better results than interpolation in practise. Some examples will give the evidence for this argument.

Preliminaries

Term structure of interest rates

Discrete case

Suppose that at time instances t=1,2,...,T we have spot interest rates (Subscript[R, 1]),...,(Subscript[R, T]). The forward rates (Subscript[f, 1]),...,(Subscript[f, T]) are defined by relation

(((1+(Subscript[R, t])))^t)=(∏_(j=1)^t)(1+(Subscript[f, j]))      or         (Subscript[f, t])=((((1+(Subscript[R, t])))^t)/(((1+(Subscript[R, t-1])))^(t-1)))-1 .

Finally, the discount factors (Subscript[v, 1]),...,(Subscript[v, T]) are defined as

(Subscript[v, t])=(1/(((1+(Subscript[R, t])))^t)) .

Any of these sequences uniquely defines the term structure of interest rates. Note that also

log(1+(Subscript[R, t]))=(1/t)(∑_(j=1)^t)log(1+(Subscript[f, j])) .

Continuous case

Here we start with forward or instant rate of interest f(t) which is a function defined on [0,T]. The quantity f(t)Δt roughly means the accrued interest during the time interval [t,t+Δt]. The spot rate of interest is defined as the integral mean of the forward rate and resembles (1) in discrete case:

R(t)=(1/t)(∫_0^t)f(τ)dτ .

The discount function v(t) is

v(t)= exp{-t R(t)}=exp{-(∫_0^t)f(τ)dτ} .

The data

Needs["NonlinearRegression`"]

This set of data come from real world. They are in the form {maturity (number of days), spot interest rate}.

plot3=ListPlot[(ir3={{7,4.74},{14,4.72833333333333},{30,4.705},{60,4.595},{90,4.53},{180,4.42},{270,4.415},{360,4.41},{730,4.51},{1095,4.68},{1460,4.85},{1825,4.99},{2190,5.09},{2555,5.19},{2920,5.27},{3285,5.33},{3650,5.38`},{4015,5.42},{4380,5.46},{4745,5.48666666666667},{5110,5.51333333333333},{5475,5.54}}),PlotStyle->{Red,PointSize[0.01]}]

int3=Interpolation[ir3];
int31=Interpolation[ir3,InterpolationOrder->1];
plotint3=Plot[int3[t],{t,7,10000}];
Show[plot3,plotint3]

Parametric models

There is a number of parametric models used for the either fitting or interpolation of empirical data coming from markets. We will deal with just two of them which, under the author's opinion are most vital.

Nelson-Siegel

Clear[nelsonsiegelspot]; (* Spot curve or zero-coupon*)
nelsonsiegelspot[b0_,b1_,b2_,c_,t_]:=b0+(((b1+b2) (1-(E^(-(t/c)))))/(t/c))-b2 (E^(-(t/c)))

The parameters have the following meaning: the function starts at level b0+b1 at time 0 and has asymptote at level b0; a local extreme is determined by parameter c.

Limit[nelsonsiegelspot[b0,b1,b2,c,t],t->0]
Limit[nelsonsiegelspot[b0,b1,b2,c,t],t->∞,Assumptions->{b0>0∧b1>0∧b2>0∧c>0}]

Manipulating spot:

Manipulate[Plot[nelsonsiegelspot[b0,b1,b2,c,t],{t,0,10000},PlotRange->{{0,10000},{3,6}},AxesOrigin->{0,3}],
{{b0,5.8,"b0"},5,6,0.1},
{{b1,-1,"b1"},2,-0.5,0.1},
{{b2,-3,"b2"},-5,-1,0.2},
{{c,400,"c"},100,1000,10}
]

(* Forward curve *)
Clear[nelsonsiegelforward];
nelsonsiegelforward[b0_,b1_,b2_,c_,t_]:=b0+b1 (E^(-t/c))+b2 (t/c) (E^(-t/c))

Manipulating forward:

Manipulate[Plot[nelsonsiegelforward[b0,b1,b2,c,t],{t,0,10000},PlotRange->{{0,10000},{3,6}},AxesOrigin->{0,3},PlotStyle->Orange],
{{b0,5.8,"b0"},5,6,0.1},
{{b1,-1,"b1"},2,-0.5,0.1},
{{b2,-3,"b2"},-5,-1,0.2},
{{c,400,"c"},100,1000,10}
]

Fitting Nelson-Siegel: on the graph, estimated spot is blue and estimated forward is orange.

fitnelsonsiegel3= FindFit[ir3,nelsonsiegelspot[b0,b1,b2,c,t],{b0,b1,b2,c},t]
curvenel3=nelsonsiegelspot[b0,b1,b2,c,t]/.fitnelsonsiegel3;
curvenelfc3[t_]:=nelsonsiegelspot[b0,b1,b2,c,t]/.fitnelsonsiegel3;
curvenel3forward=nelsonsiegelforward[b0,b1,b2,c,t]/.fitnelsonsiegel3;
Show[plot3,Plot[curvenel3forward,{t,0,10000},PlotRange->All,PlotStyle->Orange],Plot[curvenel3,{t,0,10000},PlotRange->All],PlotRange->All]

Swensson

By adding the 5th parameter in Nelson-Siegel model we get Swensson model. Note that possibly two extremes are determined by parameters c1 and c2.

Clear[swenssonforward];
swenssonforward[b0_,b1_,b2_,b3_,c1_,c2_,t_]:=b0+b1 (E^(-t/c1))+b2 (t/c1) (E^(-t/c1))+b3 (t/c2)(E^(-t/c2))
Clear[swenssonspot];
swenssonspot[b0_,b1_,b2_,b3_,c1_,c2_,t_]=(1/t)(∫_0^t)swenssonforward[b0,b1,b2,b3,c1,c2,τ]\[DifferentialD]τ//Simplify

Limit[swenssonspot[b0,b1,b2,b3,c1,c2,t],t->0,Assumptions->{b0>0&&b1>0&&b2>0&&b3>0&&c1>0&&c2>0}]
Limit[swenssonspot[b0,b1,b2,b3,c1,c2,t],t->∞,Assumptions->{b0>0&&b1>0&&b2>0&&b3>0&&c1>0&&c2>0}]

Manipulate[Plot[swenssonspot[b0,b1,b2,b3,c1,c2,t],{t,0,5000},PlotRange->{{0,5000},{3,6}},AxesOrigin->{0,3}],
{{b0,5.8,"b0"},5,6,0.1},
{{b1,-1,"b1"},2,-0.5,0.1},
{{b2,-3,"b2"},-5,-1,0.2},
{{b3,-3,"b3"},-5,-1,0.2},
{{c1,200,"c1"},100,1000,10},
{{c2,5000,"c2"},100,10000,10}
]

Manipulate[Plot[swenssonforward[b0,b1,b2,b3,c1,c2,t],{t,0,5000},PlotRange->{{0,5000},{3,6}},AxesOrigin->{0,3},PlotStyle->Orange],
{{b0,5.8,"b0"},5,6,0.1},
{{b1,-1,"b1"},2,-0.5,0.1},
{{b2,-3,"b2"},-5,-1,0.2},
{{b3,-3,"b3"},-5,-1,0.2},
{{c1,200,"c1"},100,1000,10},
{{c2,5000,"c2"},100,10000,10}
]

fitswensson3= FindFit[ir3,swenssonspot[b0,b1,b2,b3,c1,c2,t],{b0,b1,b2,b3,c1,c2},t]
(curveswe3=swenssonspot[b0,b1,b2,b3,c1,c2,t]/.fitswensson3)//Simplify;
curveswefc3[t_]:=swenssonspot[b0,b1,b2,b3,c1,c2,t]/.fitswensson3
(curveswe3forward=swenssonforward[b0,b1,b2,b3,c1,c2,t]/.fitswensson3)//Simplify;
Show[plot3,Plot[curveswe3,{t,0,10000},PlotRange->All],Plot[curveswe3forward,{t,0,10000},PlotRange->All,PlotStyle->Orange],PlotRange->All]

Nonparametric modelling

Interpolation by polynomials

Polynomials of  the degree n-1

Stupid but sometimes observable method: for data size n, use polynomial of the degree n-1. Since the maximum allowed order for interpolation is 21 we use it for the demonstration of the catastrophic behaviour of the interpolated function

intpol3=Interpolation[ir3,InterpolationOrder->Min[Length[ir3],21]];
plotpol3=Plot[intpol3[t],{t,Min[First[ir3][[1]]],Max[Last[ir3][[1]]]},PlotRange->All];
Show[plot3,plotpol3,PlotRange->All]

Manipulate[intpol3=Interpolation[ir3,InterpolationOrder->order];
Show[plotpol3=Plot[intpol3[t],{t,Min[First[ir3][[1]]],Max[Last[ir3][[1]]]},PlotRange->All],plot3,PlotRange->All],
{{order,4,"order"},1,21,1}]

Splines

Procedure for fitting data by cubic splines

SplinesJanHurt.nb
Reference 1: Späth, Helmuth: Algorithmen für elementare Ausgleichsmodelle. Oldenbourg Verlag, München Wien 1973 (pp. 59-61). Beware of a lot of errors in the formulas!!!

Reference 2 (added in proofs) : Fisher, M. and Zervos, D.: YieldCurve. In: Varian, H. R. (ed.) Computational Economics and Finance. Springer. New York 1996. (pp. 269-302).

Notice that with size n of the input data the system of 3n linear equations must be solved. Generally, this is simple task for Mathematica since the matrix of the system is five-diagonal, but rendering graphics of the spline functions usually takes a longer time. 

It is also useful to compare this method with the method of running medians (Tukey).

splinesJH::usage = "splinesJH[x, u, p, z] gives a list of cubic splines in the respective intervals given by the list x of independent observations, u being the list of dependent observations at points x, weights list p, and z as the argument of the spline functions. If p->∞, the method leads to usual numerical interpolation while for p->0 it leads to the least squares method."

plotsplines::usage = "plotsplines[x, u, p]  produces the plot of the spline function given the list x of independent variables, the list u of dependent variables with weights list p."

Clear[splinesJH];splinesJH[x_,u_,p_,z_]:=Module[{pp=p,dz,dyy,yy2,eq1,eq2,eq3,eq4,eq5,variables,sol,yyy,yyy2,dyyy,dyyy2,aa,bb,cc,dd,xk,n},
n=Length[x];
dx=Take[RotateLeft[x]-x,n-1];
yy=Array[y,n];
dyy=Take[RotateLeft[yy]-yy,n-1];
yy2=Array[y2,n];tt=Array[t,n];
eq1=Table[dx[[k-1]]*yy2[[k-1]]+dx[[k]]*yy2[[k+1]]+2*(dx[[k-1]]+dx[[k]])*yy2[[k]]==6*(dyy[[k]]/dx[[k]]-dyy[[k-1]]/dx[[k-1]]),{k,2,n-1}];
eq2={yy2[[1]]==0,yy2[[n]]==0};
eq3=Table[pp[[k]]*(u[[k]]-yy[[k]])==tt[[k]],{k,n}];
eq4={tt[[1]]==(yy2[[2]]-yy2[[1]])/dx[[1]],tt[[n]]==-((yy2[[n]]-yy2[[n-1]])/dx[[n-1]])};
eq5=Table[tt[[k]]==(yy2[[k+1]]-yy2[[k]])/dx[[k]]-(yy2[[k]]-yy2[[k-1]])/dx[[k-1]]==yy2[[k-1]]/dx[[k-1]]+yy2[[k+1]]/dx[[k]]-(1/dx[[k]]+1/dx[[k-1]])*yy2[[k]],{k,2,n-1}];
variables=Join[yy,yy2,tt];
sol=Solve[Join[eq1,eq2,eq3,eq4,eq5],variables];
yyy=Flatten[Partition[yy/.sol,1]];
yyy2=Flatten[Partition[yy2/.sol,1]];
dyyy=Take[RotateLeft[yyy]-yyy,n-1];
dyyy2=Take[RotateLeft[yyy2]-yyy2,n-1];
aa=dyyy2/(6*dx);
bb=1/2*Take[yyy2,n-1];
cc=dyyy/dx-1/6*dx*(4*bb+Take[RotateLeft[yyy2,1],n-1]);
dd=Take[yyy,n-1];xk=Take[x,n-1];
aa*(z-xk)^3+bb*(z-xk)^2+cc*(z-xk)+dd]

Clear[plotsplines];plotsplines[x_,u_,p_]:=Module[{s1,t,x1,tt,intervals,tab1,k,n,plot1,plot2},
(* Needs["PlotLegends`"];*)
s1=splinesJH[x,u,p,t];n=Length[x];x1=RotateLeft[x];tt=Table[t,{n}];intervals=Take[Transpose[{tt,x,x1}],n-1];tab1=Table[Plot[s1[[k]],Evaluate[intervals[[k]]],DisplayFunction->Identity],{k,1,n-1}];(*plot1=Show[tab1,DisplayFunction->$DisplayFunction];*)plot2=ListPlot[Transpose[{x,u}],Joined->True,PlotStyle->Red];
Show[tab1,plot2,PlotRange->Automatic]]

Sorry, the above procedure has been created when the author was junior beginner with Mathematica. The system of the equations may be created more effectively by direct symbolic calculation.

Example 1

Since the interest rates data are too smooth we will illustrate the procedure on somewhat wilder data. Red are the data, blue the estimated spline function.

n=25;
SeedRandom[13];
uu=Range[n]+20 Sin[0.5 RandomReal[{0,1},n]-0.5];
xx=Range[n];
ppp=Table[10.500,{n}];
(* plotsplines[xx,uu,ppp]  *)
Manipulate[ppp=Table[p,{Length[xx]}];plotsplines[xx,uu,ppp],{{p,0.1,"weights"},0.0001,40.0001,2}]

Example 2

xa=Range[40];ua=((Range[40])^3)/1000.-15 ((Range[40])^2)/1000+10 Sin[ ((Range[40])^2)];
pa=0.0001;
ppa=Table[pa,{40}];
(* plotsplines[xa,ua,ppa] *)
Manipulate[ppa=Table[pa,{Length[xa]}];plotsplines[xa,ua,ppa],{{pa,0.1,"weights"},0.0001,100.0001,5}]

Interest rates data

Since there is almost impossible to see the differences between the fitted curve and the original data if we use all of them, we take just first 10 pairs. The following picture demonstrates that too small weight leads to numerical problems. But we usually do not fit the data by a straight line using cubic splines function.

x3=Take[Transpose[ir3][[1]],10];
u3=Take[Transpose[ir3][[2]],10];
p3=0.000000010;
pp3=Table[p3,{Length[x3]}];plotsplines[x3,u3,pp3]

Next we see that sensitivity of the shape of the fitted curve on the weighting parameter.

Manipulate[x3=Take[Transpose[ir3][[1]],10];u3=Take[Transpose[ir3][[2]],10];pp3=Table[p3,{Length[x3]}];plotsplines[x3,u3,pp3] ,{{p3,0.1,"weights"},0.0000001,0.050000001,0.001}]

A comparison

Mean squared error (MSE)

Mean squared error is defined in a standard way. Here x are independent, y dependent observations at points x, fitfc is the fitted function.

Clear[mse];mse[x_,y_,fitfc_]:=If[Length[x]==Length[y],(1/(Length[x]))(∑_(i=1)^(Length[x]))(((y[[i]]-fitfc[x[[i]]]))^2),Print["Incorrect data"]]

Nelson-Siegel

msenel=mse[Transpose[ir3][[1]],Transpose[ir3][[2]],curvenelfc3]

Swensson

msesve=mse[Transpose[ir3][[1]],Transpose[ir3][[2]],curveswefc3]

As expected, the more parameters, the better fitting.

Integrated mean squared error

From data the interpolation function is calculated and used as a numaraire for the comparison. Implicitly, the order of interpolation is 1.

imse[x_,y_,fitfc_,order_:1]:=Module[{interpoldata,xmin,xmax},
interpoldata=Interpolation[Transpose[{x,y}],InterpolationOrder->order];
xmin=Min[x];xmax=Max[x];
(1/(xmax-xmin))NIntegrate[(((interpoldata[τ]-fitfc[τ]))^2),{τ,xmin,xmax}]
]

Nelson-Siegel

imsenel=imse[Transpose[ir3][[1]],Transpose[ir3][[2]],curvenelfc3]

Swenson

imsesve=imse[Transpose[ir3][[1]],Transpose[ir3][[2]],curveswefc3]

Interpolating polynomials

Table[{order,imse[Transpose[ir3][[1]],Transpose[ir3][[2]],Interpolation[ir3,InterpolationOrder->order]]},{order,0,9}]//TableForm

Acknowledgement

The research has been supported by Grant of Czech Ministry of Education MSM 0021620839. The author is very grateful to Mark Fisher for appointing his article on splines [9].

References

[1] Shaw, W.: Modelling Financial Derivatives with Mathematica. Cambridge University Press. Cambridge, 1998.
[2] UnRisk PRICING ENGINE. MathConsult GmbH. Linz, 2006.
[3] Maclachlan, F.: Price-Yield Curve. Wolfram Demonstrations Project.  http://demonstrations.wolfram.com/PriceYieldCurve/
[4] Hurt, J. and Wertz, W.: Asymptotic properties of the invariantly optimal estimator of reliability in the exponential case. Statistics & Decision 1 (1983). 197-204.
[5] Nelson, C. R. and Siegel, A. F.: Parsimonious modeling of yield curves. J. of Business 60 (1987). 473-489.
[6] Hawkins, R. J. et al.: Ab initio yield curve dynamics. Physics Letters A 344 (2005). 317-323.
[7] Linton, O. et al.: Yield curves. J. of Econometrics 105 (2001). 185-223.
[8] Ioannides, M.: A comparison of yield curve estimation techniques using UK data. J. of Banking & Finance 27 (2003). 1-26.
[9] Fisher, M. and Zervos, D.: YieldCurve. In: Varian, H. R. (ed.) Computational Economics and Finance. Springer. New York (1996). pp. 269-302.