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