reading the paper from Joshi, Yang (2009) I decided to implement the algorithm suggested for the Euler drift freezing scheme. But it won't work when I compare my results with the ones provided in table 3.1.
Here is the Matlab code: it retrieves all the parameters we need to implement the step. I use initial suggested correlation matrix and displacement coefficients.
- Code: Select all
function [barA, a, mu]=fun(SR)
% From a state of the vector SR_1,... SR_10 we find the parameters
% barA = discounted annuities, a=pseudo-square root of Cov, mu=drifts
rho=zeros(10);
mu=zeros(1,10);
for i=1:10
for j=1:10
rho(i,j)=exp(-0.05*abs(i-j));
end
end
%IV=csvread('IV.csv');
%alpha=csvread('alpha.csv');
IV=.2; %
alpha=.02*ones(1,10); %displacements all set at 2%
sigma=zeros(1,10); %this is sigma^alpha
for i=1:10
sigma(i)=IV*SR(i)/(SR(i)+alpha(i)); %Rebonato approximation
end
covmat=diag(sigma)*rho*diag(sigma); %this is where I could be wrong
a=chol(covmat)'; %pseudo square root..
P=zeros(1,10); %bond prices deflated by P0
barA=zeros(1,10);
for i=1:10
if i==1
aaa=0;
else
aaa=barA(i-1);
end
P(i)=(1-SR(i)*aaa)/(1+SR(i)); %formula (2.10)
barA(i)=aaa+P(i); %(2.8)
end
cv=zeros(10); %cv(i,j)= cross variation Z_k, A_j
for k=1:10
for i=1:10
if i==1
aaa=0;
else
aaa=cv(k,i-1);
end
temp=aaa-a(i,k)*(SR(i)+alpha(i))*barA(i); %(2.6)
cv(k,i)=temp/(1+SR(i));
end
end
the step itself is done like
- Code: Select all
function SR=step(S)
%we implement the formula (2.11)
SR=S;
alpha=0.02*ones(1,10);
[barA,a,mu]=fun(S);
Z=randn(1,10);
for i=1:10
temp=0;
for k=1:10
temp=temp+a(i,k)*Z(k)-a(i,k)^2/2;
end
SR(i)=(S(i)+alpha(i))*exp(mu(i)+temp)-alpha(i);
end
