% This is va2fgrin, a 2-factor Vasicek model for Germany spot curve estimated with identified factors
% 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 below
function lf=va2fgrin(para);
% loading data
load 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 r77 r88 r99 % 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 r77 r88 r99 r1010 r1111 r1212 r1313 r1414 r1515 r1616 r1717 r1818 r1919 s10 s20;
global AAr HH1r HH2r A H F R Q S10 Q0 p10 vp10 z;
% transposing the data column vector
Y=datgerin';
% computing the number of observations
obs=length(Y(1,:));
% creating a row vector of 1's with No. of columns equal to the number of observations
X=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'1 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 chapter
c1=para(5); c2=para(6);
% variances of the residuals in observation equation errors
d1=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 less than one;
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
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 parameters
P=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 (48)
AAr=zeros(1,NN);
HH1r=zeros(1,NN);
HH2r=zeros(1,NN);
%============================Recursive Restrictions on the parameters====================================
for n=2:NN
% Equation (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 (25)
HH1r(n)=1+phi1*HH1r(n-1);
HH2r(n)=1+phi2*HH2r(n-1);
end
%=============================Defining the matrices for estimation====================================
% Estimated parameters for the maturities considered in the sample
% 1m
A1=AAr(2);
H11=HH1r(2);
H12=HH2r(2);
% 3m
A2=AAr(4);
H21=HH1r(4);
H22=HH2r(4);
% 1y
A3=AAr(13);
H31=HH1r(13);
H32=HH2r(13);
% 2y
A4=AAr(25);
H41=HH1r(25);
H42=HH2r(25);
% 3y
A5=AAr(37);
H51=HH1r(37);
H52=HH2r(37);
% 4y
A6=AAr(49);
H61=HH1r(49);
H62=HH2r(49);
% 5y
A7=AAr(61);
H71=HH1r(61);
H72=HH2r(61);
% 7y
A8=AAr(85);
H81=HH1r(85);
H82=HH2r(85);
% 10y
A9=AAr(121);
H91=HH1r(121);
H92=HH2r(121);
% creating a column vector (18 rows) of A's (independent term of the observation equation in (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 (48))
H=zeros(19,2);
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 (49))
F=zeros(2,2);
F(1,1)=phi1;
F(2,2)=phi2;
% creating the variance-covariance matrix of the observation equation errors in (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;
R(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 errors in (49)
Q=zeros(2,2);
Q(1,1)=sgm1^2;
Q(2,2)=sgm2^2;
% creating the matrix of the volatility curve in equation (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 (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;
%=========================================== Calculate 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 factors
S10(1,1)=s10;
S10(2,1)=s20;
% 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 (A7)
vp10=inv(eye(rr^2)-kron(F,F))*reshape(Q,rr*rr,1);
% Unconditional variance of the state vector (in matrix form)
p10=reshape(vp10,rr,rr)';
% Initial value for the likelihood function
lf=0;
%====================== Derivation of the KF ===================================
for t=1:obs
% Forecasting Yt|t-1 (equation (A10))
y10=(A*X(:,t)+H*S10(:,t))*1200;
% Creating matrix of estimated Y
YE(:,t:t)=y10;
% covariance of Yt-Yt|t-1 (equation (A12))
O10=H*p10*H'+R;
kk=length(R(:,1));
%====================== Calculate the likelihood function ======================
% Forecast error (equation (A11))
res=Y(:,t)-y10;
% log likelihood function (equation (A23))
ff=-0.5*log(2*pi*det(O10))-0.5*res'*inv(O10)*res;
% Sum of the negative of the log likelihood function
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 (A18))
S=S10(:,t)+p10*H'*inv(O10)*(Y(:,t)-y10);
% Covariance of St|t (equation (A19))
p=p10-p10*H'*inv(O10)*H*p10;
% Forecast St+1|t (equation (A20))
S10(:,t+1)=F*S;
% Update Pt+1|t (equation (A21))
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;
T(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 matrix
YE1=YE(1:9,:);
% Volatility matrix
YE2=YE(10:18,:);
YE3=YE(19:19,:);
YET=YE1';
VET=YE2;
IET=YE3';
% Observation equation matrices
AArT=AAr';
HH1rT=HH1r';
HH2rT=HH2r';
L=[AArT,HH1rT,HH2rT];
%VT=V';
%TT=T';
S10T=S10';
%============================================== Plots =========================================================================
% Yields
figure(1)
plot(YET)
% Volatilities
figure(2)
plot(VET)
% Factors
figure(3)
plot(S10T)
% Estimated inflation
figure(4)
plot(IET)
%============================================== Saving files ==================================================================
% Estimated Yields
save c:\matlab\V2YETgri YET -ascii;
% Estimated Inflation
save c:\matlab\V2IETgri IET -ascii;
% Factors
save c:\matlab\V2S10gri S10T -ascii;
% Factor loadings
save c:\matlab\V2loadgi L -ascii;
% Volatilities
save c:\matlab\V2Volgri VET -ascii;
% Term premium
save c:\matlab\V2tprgri TT -ascii;