t isalso independent from St and Ŝt|t−1 , from equations (3.48) and (3.61), the forecast errorcorresponds to:yt − E(yt |Xt , ζt−1 ) = (AXt + H St + wt ) − (AXt + H Ŝt|t−1 ) = H (St − Ŝt|t−1 ) + wt(3.62)From (3.62) the conditional variance–covariance matrix of the estimation error of theobservation vector will be:34E{[yt − E(yt |Xt , ζt−1 )][yt − E(yt |Xt , ζt−1 )] }= E{[H (St − Ŝt|t−1 ) + wt ][H (St − Ŝt|t−1 ) + wt ] }= HE [(St − Ŝt|t−1 )(St − Ŝt|t−1 ) ]H  + E(wt wt )= HP t|t−1 H  + R
(3.63)
After the updates of the mean and variance–covariance matrices of the dependentvariables, the log-likelihood function is computed. In our Kalman ﬁlter exercise, thelikelihood function is deﬁned as follows:35f (yt |It−1 ) = (2π )−1/2 |HP t|t−1 H  + R|−1/2× exp − 12 (yt − A − H Ŝt|t−1 ) (HP t|t−1 H  + R)−1 (yt − A − H Ŝt|t−1 )(3.64)34
35
Notice that from equations (3.48) and (3.61), the forecast error corresponds to:yt − E(yt |Xt , ζt−1 ) = (AXt + H St + wt ) − (AXt + H Ŝt|t−1 ) = H (St − Ŝt|t−1 ) + wt
As the latter depends on the residuals of the observation equations, these are previously deﬁned in the Matlabroutines.
92
Applied Quantitative Methods for Trading and Investment
As usual, the optimisation procedure will be performed as a minimisation problem.Thus, the negative of the log-likelihood function is computed:− log L(YT ) = −
T
log f (yt |It−1 )
(3.65)
t=1
It can be shown that the distribution of St given yt , Xt and ζt−1 is N (Ŝt|t , Pt|t ), whereŜt|t and Pt|t are respectively the optimal forecast of St given Pt|t and the mean square errorof this forecast, corresponding to the following updating equations of the Kalman ﬁlter.After the initial estimate of the dependent variables, the mean and the variance–covariancematrices of the Kalman ﬁlter may be updated as follows:Ŝt|t = Ŝt|t−1 + Pt|t−1 H  (HP t|t−1 H  + R)−1 [yt − (AXt + Ht|t−1 )]
(3.66)
Pt|t = Pt|t−1 − Pt|t−1 H  (HP t|t−1 H  + R)−1 HP t|t−1
(3.67)
Following this update, a new estimate for these matrices can be obtained, generalising(3.58) and (3.59):36Ŝt+1|t = C + F Ŝt|t= C + F {Ŝt|t−1 + Pt|t−1 H  (HP t|t−1 H  + R)−1 [yt − (AXt + Ht|t−1 )]}= C + F Ŝ t|t−1 + FP t|t−1 H  (HP t|t−1 H  + R)−1 [yt − (AXt + Ht|t−1 )]
(3.68)
Pt+1|t = FP t|t F  + GQ t G= F [Pt|t−1 − Pt|t−1 H  (HP t|t−1 H  + R)−1 HP t|t−1 ]F  + GQ t G= FP t|t−1 F  − FP t|t−1 H  (HP t|t−1 H  + R)−1 HP t|t−1 F  + GQ t G
(3.69)
At the end of each update of the Kalman ﬁlter, the volatility and the term premiumcurves are computed. Considering that only Gaussian models were ﬁtted, these curvesmove only when different values for the parameters are considered, as they do not dependon the factor values. Finally, concluding the routines, there are two blocks where theresults are stored in matrices and ﬁles and charts for the yields, volatility curves, factorsand inﬂation (in the identiﬁed models) are produced. The Matlab routine written to runthe Kalman ﬁlter for the two-factor identiﬁed model is presented next.% This is va2fgrin, a 2-factor Vasicek model, with an identified% factor, for Germany spot curve estimated using the Kalman Filter% By Jorge Barros Luis% Estimation on 1m, 1y, 2y, 3y, 4y, 5y, 7y, 10y% ==============================GENERAL PROCEDURES======================% computing the log-likelihood function defined belowfunction lf=va2fgrin(para);The matrix FP t|t−1 H  (HP t|t−1 H  + R)−1 is usually known as the gain matrix, since it determines the updatein Ŝt+1|t due to the estimation error of yt . Equation (3.68) is known as a Ricatti equation.36
Modelling the Term Structure of Interest Rates
93
% loading dataload c:\matlab\datgerin.txt;% creating the original parameters to be estimated% r11 to r1919 are the variances of the residuals in observation% equations% s10 and s20 are the starting values of the factors (subtracted by the% respective means)global delta phi1 phi2 sgm1 sgm2 lmda1 lmda2 r11 r22 r33 r44 r55 r66 r77r88 r99 r1010 r1111 r1212 r1313 r1414 r1515 r1616 r1717 r1818 r1919 s10s20;global AAr HH1r HH2r A H F R Q S10 Q0 p10 vp10 z;% transposing the data column vectorY=datgerin’;% computing the number of observationsobs=length(Y(1,:));% creating a row vector of 1’s with No. of columns equal to the number% of observationsX=ones(1,obs);% computing the longest term for the restrictions on the parameters (in% months)NN=10*12+1;%=========DEFINING THE TRANSFORMED OR UNRESTRICTED PARAMETERS==========% these parameters are defined as to impose restrictions on the% original parameters (e.g. sigma’s have to be positive, phi’s have to% be between 0 and 1,% one of the lambda’s has to be negative in order to have positive term% premia and the variances of the residuals have to be positive)% persistence parameters (phi1 and phi2)a1=para(1); a2=para(2);% volatilities (sigma 1 and sigma2)b1=para(3); b2=para(4);% lambda’s (the market price of risks), defined as the product between% the lambda’s and the sigma’s in the chapterc1=para(5); c2=para(6);% variances of the residuals in observation equation errorsd1=para(7); d2=para(8); d3=para(9); d4=para(10); d5=para(11);d6=para(12); d7=para(13); d8=para(14); d9=para(15);d10=para(16); d11=para(17); d12=para(18); d13=para(19); d14=para(20);d15=para(21); d16=para(22); d17=para(23); d18=para(24); d19=para(25);% conditional mean of short-term interest rate;e1=para(26);%==================DEFINING THE ORIGINAL PARAMETERS==================% conditional mean of short-term interest rate;delta=exp(e1);% constrain persistence parameters (phi’s) to be positive and <1;phi1=0.999/(1+exp(-a1));phi2=0.999/(1+exp(-a2));% constrain volatilities (sigma’s) to be positive;sgm1=exp(b1);sgm2=exp(b2);% constrain the 2st market price of risk (lamda2) to be negative;lmda1=c1;lmda2=-exp(c2);% constrain variances of the observation equation errors
94
Applied Quantitative Methods for Trading and Investment
r11=exp(d1);r22=exp(d2);r33=exp(d3);r44=exp(d4);r55=exp(d5);r66=exp(d6);r77=exp(d7);r88=exp(d8);r99=exp(d9);r1010=exp(d10);r1111=exp(d11);r1212=exp(d12);r1313=exp(d13);r1414=exp(d14);r1515=exp(d15);r1616=exp(d16);r1717=exp(d17);r1818=exp(d18);r1919=exp(d19);% initial values for the state variables;s10=0;s20=0;%======================================================================% creating a row vector of the original parametersP=zeros(1,6);P(1,1)=phi1;P(1,2)=phi2;P(1,3)=sgm1;P(1,4)=sgm2;P(1,5)=lmda1;P(1,6)=lmda2;%===========DEFINING ZERO ROW VECTORS WITH DIMENSION=NN================%Matrices of observation equations in (3.48)AAr=zeros(1,NN);HH1r=zeros(1,NN);HH2r=zeros(1,NN);%=============Recursive Restrictions on the Factor Loadings============for n=2:NN% Equation (3.24)AAr(n)=delta+AAr(n-1)+0.5*(lmda1^2+lmda2^2-(lmda1+HH1r(n-1)*sgm1)^2(lmda2+HH2r(n-1)*sgm2)^2);% Equation (3.25)HH1r(n)=1+phi1*HH1r(n-1);HH2r(n)=1+phi2*HH2r(n-1);end%==============Defining the matrices for estimation====================% Estimated factor loadings for the maturities considered in the sample% 1mA1=AAr(2);H11=HH1r(2);H12=HH2r(2);% 3mA2=AAr(4);H21=HH1r(4);
Modelling the Term Structure of Interest RatesH22=HH2r(4);% 1yA3=AAr(13);H31=HH1r(13);H32=HH2r(13);% 2yA4=AAr(25);H41=HH1r(25);H42=HH2r(25);% 3yA5=AAr(37);H51=HH1r(37);H52=HH2r(37);% 4yA6=AAr(49);H61=HH1r(49);H62=HH2r(49);% 5yA7=AAr(61);H71=HH1r(61);H72=HH2r(61);% 7yA8=AAr(85);H81=HH1r(85);H82=HH2r(85);% 10yA9=AAr(121);H91=HH1r(121);H92=HH2r(121);% creating a column vector (18 rows) of A’s (independent term of the% observation equation in (3.48))A=zeros(19,1);A(1,1)=A1/1;A(2,1)=A2/3;A(3,1)=A3/12;A(4,1)=A4/24;A(5,1)=A5/36;A(6,1)=A6/48;A(7,1)=A7/60;A(8,1)=A8/84;A(9,1)=A9/120;A(10,1)=1200/(1^2)*(H11^2*sgm1^2+H12^2*sgm2^2);A(11,1)=1200/(3^2)*(H21^2*sgm1^2+H22^2*sgm2^2);A(12,1)=1200/(12^2)*(H31^2*sgm1^2+H32^2*sgm2^2);A(13,1)=1200/(24^2)*(H41^2*sgm1^2+H42^2*sgm2^2);A(14,1)=1200/(36^2)*(H51^2*sgm1^2+H52^2*sgm2^2);A(15,1)=1200/(48^2)*(H61^2*sgm1^2+H62^2*sgm2^2);A(16,1)=1200/(60^2)*(H71^2*sgm1^2+H72^2*sgm2^2);A(17,1)=1200/(84^2)*(H81^2*sgm1^2+H82^2*sgm2^2);A(18,1)=1200/(120^2)*(H91^2*sgm1^2+H92^2*sgm2^2);A(19,1)=0;% creating a 18*2 matrix of H’s (coefficient matrix of observation% equation in (3.48))H=zeros(19,2);
95
96
Applied Quantitative Methods for Trading and Investment
H(1,1)=H11/1;H(1,2)=H12/1;H(2,1)=H21/3;H(2,2)=H22/3;H(3,1)=H31/12;H(3,2)=H32/12;H(4,1)=H41/24;H(4,2)=H42/24;H(5,1)=H51/36;H(5,2)=H52/36;H(6,1)=H61/48;H(6,2)=H62/48;H(7,1)=H71/60;H(7,2)=H72/60;H(8,1)=H81/84;H(8,2)=H82/84;H(9,1)=H91/120;H(9,2)=H92/120;H(10,1)=0;H(10,2)=0;H(11,1)=0;H(11,2)=0;H(12,1)=0;H(12,2)=0;H(13,1)=0;H(13,2)=0;H(14,1)=0;H(14,2)=0;H(15,1)=0;H(15,2)=0;H(16,1)=0;H(16,2)=0;H(17,1)=0;H(17,2)=0;H(18,1)=0;H(18,2)=0;H(19,1)=0;H(19,2)=1/phi2;% creating a 2*2 matrix of F’s (coefficient matrix of state equations% in (3.49))F=zeros(2,2);F(1,1)=phi1;F(2,2)=phi2;% creating the variance-covariance matrix of the observation equation% errors in (3.48)R=zeros(19,19);R(1,1)=r11;R(2,2)=r22;R(3,3)=r33;R(4,4)=r44;R(5,5)=r55;R(6,6)=r66;R(7,7)=r77;R(8,8)=r88;
Modelling the Term Structure of Interest RatesR(9,9)=r99;R(10,10)=r1010;R(11,11)=r1111;R(12,12)=r1212;R(13,13)=r1313;R(14,14)=r1414;R(15,15)=r1515;R(16,16)=r1616;R(17,17)=r1717;R(18,18)=r1818;R(19,19)=r1919;% creating the variance-covariance matrix of the state equation errorsin (3.49)Q=zeros(2,2);Q(1,1)=sgm1^2;Q(2,2)=sgm2^2;% creating the matrix of the volatility curve in Equation (3.29), in %V=zeros(9,obs+1);V(1,1)=1/(1^2)*(H11^2*sgm1^2+H12^2*sgm2^2)*1200^2;V(2,1)=1/(3^2)*(H21^2*sgm1^2+H22^2*sgm2^2)*1200^2;V(3,1)=1/(12^2)*(H31^2*sgm1^2+H32^2*sgm2^2)*1200^2;V(4,1)=1/(24^2)*(H41^2*sgm1^2+H42^2*sgm2^2)*1200^2;V(5,1)=1/(36^2)*(H51^2*sgm1^2+H52^2*sgm2^2)*1200^2;V(6,1)=1/(48^2)*(H61^2*sgm1^2+H62^2*sgm2^2)*1200^2;V(7,1)=1/(60^2)*(H71^2*sgm1^2+H72^2*sgm2^2)*1200^2;V(8,1)=1/(84^2)*(H81^2*sgm1^2+H82^2*sgm2^2)*1200^2;V(9,1)=1/(120^2)*(H91^2*sgm1^2+H92^2*sgm2^2)*1200^2;% creating the matrix of the term premium in Equation (3.30), in %T=zeros(9,obs+1);T(1,1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H11*sgm1)^2+(lmda2+H12*sgm2)^2))*1200;T(2,1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H21*sgm1)^2+(lmda2+H22*sgm2)^2))*1200;T(3,1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H31*sgm1)^2+(lmda2+H32*sgm2)^2))*1200;T(4,1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H41*sgm1)^2+(lmda2+H42*sgm2)^2))*1200;T(5,1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H51*sgm1)^2+(lmda2+H52*sgm2)^2))*1200;T(6,1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H61*sgm1)^2+(lmda2+H62*sgm2)^2))*1200;T(7,1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H71*sgm1)^2+(lmda2+H72*sgm2)^2))*1200;T(8,1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H81*sgm1)^2+(lmda2+H82*sgm2)^2))*1200;T(9,1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H91*sgm1)^2+(lmda2+H92*sgm2)^2))*1200;%===================== Starting the Kalman Filter ====================% creating the matrix of the unobservable factors (2*obs+1, as the% starting values are included)S10=zeros(2,obs+1);% assuming starting values for the factorsS10(1,1)=s10;S10(2,1)=s20;
97
98
Applied Quantitative Methods for Trading and Investment
% length of the vector of the initial values for the factors;rr=length(S10(:,1));% Unconditional variance of the state vector (in vector form)% Equation (3.60)vp10=inv(eye(rr^2)-kron(F,F))*reshape(Q,rr*rr,1);% Unconditional variance of the state vector (in matrix form)% Equation (3.59)p10=reshape(vp10,rr,rr)’;% Initial value for the likelihood functionlf=0;%====================== Derivation of the KF ==========================for t=1:obs% Forecasting Yt|t-1 (equation (3.61))y10=(A*X(:,t)+H*S10(:,t))*1200;% Creating matrix of estimated YYE(:,t:t)=y10;% covariance of Yt-Yt|t-1 (equation (3.63))O10=H*p10*H’+R;kk=length(R(:,1));%================== Calculate the likelihood function =================% Forecast error (equation (3.62))res=Y(:,t)-y10;% log likelihood function (equation (3.64))ff=-0.5*log(2*pi*det(O10))-0.5*res’*inv(O10)*res;% Sum of the negative of the log likelihood function (equation (3.65))lf=lf-ff;%======== Continuing updating the state vector and covariances ========% Update St|t-1, in order to get the expected value of the state% vector, given the information available up to t (equation (3.66))S=S10(:,t)+p10*H’*inv(O10)*(Y(:,t)-y10);% Covariance of St|t (equation (3.67))p=p10-p10*H’*inv(O10)*H*p10;% Forecast St+1|t (equation (3.68))S10(:,t+1)=F*S;% Update Pt+1|t (equation (3.69))p10=F*p*F’+Q;%===================Calculate the volatility curve =====================V(1,t+1)=1/(1^2)*(H11^2*sgm1^2+H12^2*sgm2^2)*1200^2;V(2,t+1)=1/(3^2)*(H21^2*sgm1^2+H22^2*sgm2^2)*1200^2;V(3,t+1)=1/(12^2)*(H31^2*sgm1^2+H32^2*sgm2^2)*1200^2;V(4,t+1)=1/(24^2)*(H41^2*sgm1^2+H42^2*sgm2^2)*1200^2;V(5,t+1)=1/(36^2)*(H51^2*sgm1^2+H52^2*sgm2^2)*1200^2;V(6,t+1)=1/(48^2)*(H61^2*sgm1^2+H62^2*sgm2^2)*1200^2;V(7,t+1)=1/(60^2)*(H71^2*sgm1^2+H72^2*sgm2^2)*1200^2;V(8,t+1)=1/(84^2)*(H81^2*sgm1^2+H82^2*sgm2^2)*1200^2;V(9,t+1)=1/(120^2)*(H91^2*sgm1^2+H92^2*sgm2^2)*1200^2;%================== Calculate the term premium curve ==================T(1,t+1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H11*sgm1)^2+(lmda2+H12*sgm2)^2))*1200;T(2,t+1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H21*sgm1)^2+(lmda2+H22*sgm2)^2))*1200;T(3,t+1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H31*sgm1)^2+(lmda2+H32*sgm2)^2))*1200;
Modelling the Term Structure of Interest RatesT(4,t+1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H41*sgm1)^2+(lmda2+H42*sgm2)^2))*1200;T(5,t+1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H51*sgm1)^2+(lmda2+H52*sgm2)^2))*1200;T(6,t+1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H61*sgm1)^2+(lmda2+H62*sgm2)^2))*1200;T(7,t+1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H71*sgm1)^2+(lmda2+H72*sgm2)^2))*1200;T(8,t+1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H81*sgm1)^2+(lmda2+H82*sgm2)^2))*1200;T(9,t+1)=(0.5*(lmda1^2+lmda2^2)0.5*((lmda1+H91*sgm1)^2+(lmda2+H92*sgm2)^2))*1200;end%======================= Transposing matrices ========================YT=Y’;% Yield matrixYE1=YE(1:9,:);% Volatility matrixYE2=YE(10:18,:);YE3=YE(19:19,:);YET=YE1’;VET=YE2;IET=YE3’;% Observation equation matricesAArT=AAr’;HH1rT=HH1r’;HH2rT=HH2r’;L=[AArT,HH1rT,HH2rT];%VT=V’;%TT=T’;S10T=S10’;%============================== Plots ================================% YieldsFigure(1)plot(YET)% VolatilitiesFigure(2)plot(VET)% FactorsFigure(3)plot(S10T)% Estimated inflationFigure(4)plot(IET)%============================ Saving files =========================% Estimated Yieldssave c:\matlab\V2YETgri YET -ascii;% Estimated Inflationsave c:\matlab\V2IETgri IET -ascii;% Factorssave c:\matlab\V2S10gri S10T -ascii;% Factor loadingssave c:\matlab\V2loadgi L -ascii;% Volatilities
99
100
Applied Quantitative Methods for Trading and Investment
save c:\matlab\V2Volgri VET -ascii;% Term premiumsave c:\matlab\V2tprgri TT -ascii;
This routine runs the Kalman ﬁlter for any given set of parameters. In order to use itin Matlab, an instruction must be written in the software editor comprising the name ofthe ﬁle and the vector of the restricted parameters between brackets, as follows:va2fgrin([4.035,4.066,-6.465,-6.822,3.098,1.492,-0.690,-0.170,-1.780,-1.236,-0.681,0.610,-0.083,-0.147, -1.068,0.415,-0.317,-0.079,-0.378,0.735,0.805,0.511,0.392,-0.561,-2.794,-5.302])Though running the Kalman ﬁlter for given sets of parameters may be interesting toassess the behaviour of the model, a more relevant exercise is the estimation of theoptimal set of