I have completely re-worked the structure of the program. I have removed looping from monte carlo and written all terms in a serial fashion both for 2nd order and fourth order monte carlo. I have also re-worked the part that calculates the drift, its derivatives and volatility in terms of hermite polynomials. With the new structure, you can easily know about the contribution to drift and volatility from a certain derivative as the indexing is based on derivatives of drift and volatility. Here is the new program
function [TrueMean,ItoHermiteMean,MCMean] = ItoHermiteMethodWilmott08Q40(x0,theta,kappa,gamma,sigma0,T)
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=mu1 x(t)^beta1 dt + mu2 x(t)^beta2 dt +sigma x(t)^gamma dz(t)
%I have not directly simulated the SDE but simulated the transformed
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.
%In this version, the mean correction has been disabled. If you are
%simulating standard mean-reverting SV type SDEs, please enable the mean
%correction by uncommenting the appropriate line in the body of code below.
dt=.125/16; % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
%T=5;
Tt=128*T;%16*16; % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
OrderA=4; %
OrderM=4; %
%dtM=.125/8;%Monte carlo time interval size dtM.
%TtM=8*8;%Monte carlo number of simulation intervals.
dtM=dt*2;
TtM=Tt/2;
dNn=.2/1; % Normal density subdivisions width. would change with number of subdivisions
Nn=48; % No of normal density subdivisions
NnMidl=25;%One half density Subdivision left from mid of normal density(low)
NnMidh=26;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;
%theta=mu1/(-mu2);
%kappa=-mu2;
%x0=.25; % starting value of SDE
beta1=0.0;
beta2=1.0; % Second drift term power.
%gamma=.95;%50; % volatility power.
%kappa=.5; %mean reversion parameter.
%theta=.075;%mean reversion target
%you can specify any general mu1 and mu2 and beta1 and beta2.
%The present program is not perfect but quite reasonably good.
%if there is only one term in drift with beta=0, there will be minor
%errors. This case is yet not covered.
mu1=+1*theta*kappa; %first drift coefficient.
mu2=-1*kappa; % Second drift coefficient.
%sigma0=1.0;%Volatility value
alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;
w(1:Nn)=x0^(1-gamma)/(1-gamma);
Z(1:Nn)=(((1:Nn)-5.5)*dNn-NnMid);
Z
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);
for nn=2:Nn-1
ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1)-normcdf(.5*Z(nn)+.5*Z(nn-1),0,1);
%Above calculate probability mass in each probability subdivision.
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigma11(1:OrderA+1)=0;
mu11(1:OrderA+1)=0;
mu22(1:OrderA+1)=0;
sigma22(1:OrderA+1)=0;
% index 1 correponds to zero level since matlab indexing starts at one.
sigma11(1)=1;
mu11(1)=1;
mu22(1)=1;
sigma22(1)=1;
for k=1:(OrderA+1)
if sigma0~=0
sigma11(k)=sigma0^(k-1);
end
if mu1 ~= 0
mu11(k)=mu1^(k-1);
end
if mu2 ~= 0
mu22(k)=mu2^(k-1);
end
if sigma0~=0
sigma22(k)=sigma0^(2*(k-1));
end
end
%Ft(1:TtM+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0; %General time powers on hermite polynomials
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers on coefficients of hermite polynomials.
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers for bessel transformed coordinates.
%YCoeff0 and YCoeff are coefficents for original coordinates monte carlo.
%YqCoeff0 and YqCoeff are bessel/lamperti version monte carlo.
YCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
YqCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
%Pre-compute the time and power exponent values in small multi-dimensional arrays
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma); %expand y^alpha where alpha=1;
YqCoeff = ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);%expand y^alpha1 where alpha1=(1-gamma)
YqCoeff=YqCoeff/(1-gamma); %Transformed coordinates coefficients have to be
%further divided by (1-gamma)
for k = 0 : (OrderA)
for m = 0:k
l4 = k - m + 1;
for n = 0 : m
l3 = m - n + 1;
for j = 0:n
l2 = n - j + 1;
l1 = j + 1;
%Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
- (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
- (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
YCoeff0(l1,l2,l3,l4) =YCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
YqCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wnStart=1;%
d2wdZ2(1:Nn)=0;
dwdZ(1:Nn)=0;
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
tic
for tt=1:Tt
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
%[wMu0dt,dwMu0dtdw,d2wMu0dtdw2,c1,c2] = CalculateDriftAndDerivatives(w,wnStart,Nn,mu1,mu2,beta1,beta2,sigma0,gamma,dt);
[wMu0dt,dwMu0dtdw,d2wMu0dtdw2,c1,c2,c3,c4] = CalculateDriftAndVolOrder4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt)
dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn) ;% ...
dw2(wnStart:Nn)=dw(wnStart:Nn).^2;
d2w(wnStart:Nn)=c2(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
d3w(wnStart:Nn)=c3(wnStart:Nn).*(Z(wnStart:Nn).^3-3*Z(wnStart:Nn));
w(isnan(w)==1)=0;
wMu0dt(isnan(wMu0dt)==1)=0;
wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[w0,Sigma00,Sigma11,Sigma22,Sigma33,Sigma1,Sigma2,Sigma3,Sigma0] = CalculateHermiteRepresentationOfDensityOrderThreeNew(w,Z,wnStart,Nn,dNn,NnMidh,NnMidl);
d2wdZ2(wnStart:Nn)=d2wdZ2(wnStart:Nn)+ d2wMu0dtdw2(wnStart:Nn).*dwdZ(wnStart:Nn).^2+dwMu0dtdw(wnStart:Nn).*d2wdZ2(wnStart:Nn);%+ ...
dwdZ(wnStart:Nn)=(Sigma00)+Sigma1(wnStart:Nn).*2.*Z(wnStart:Nn);
%dwdZN(wnStart+1:Nn-1)=(w(wnStart+2:Nn)-w(wnStart:Nn-2))/(2*dNn);
%dwdZN(wnStart)= InterpolateOrderN(4,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),0,dwdZN(wnStart+1),dwdZN(wnStart+2),dwdZN(wnStart+3),dwdZN(wnStart+4),0);
H2(wnStart:Nn)=(Z(wnStart:Nn).^2-1);
CTerm(wnStart:Nn)=(d2wMu0dtdw2(wnStart:Nn).*dwdZ(wnStart:Nn).^2 - ...
1*dwMu0dtdw(wnStart:Nn).*d2wdZ2(wnStart:Nn));
Correction1(wnStart:Nn)=-sqrt(3)/2/(22/7)^2.*CTerm(wnStart:Nn).*H2(wnStart:Nn);
if(tt==1)
w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+ ...
sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)));
end
if(tt>1)
[C0,C1,C2,C3,w0,C20,wMean,ZMean,nnMean,Cm1,Cm2,Cm3] = CalculateHermitesOrderThree(w,Z,ZProb,wnStart,Nn,dNn,NnMidh,NnMidl);
%Below is hermite decomposition of w
%w(wnStart:Nn)=C0+C1(wnStart:Nn).*Z(wnStart:Nn)+C2(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
%Below is hermite decomposition of wMean with associated ZMean.
%wMean=C0+Cm1.*ZMean+Cm2.*(ZMean^2-1);
%[mC0,mC1,mC2,mw0,mC20] = CalculateHermitesOrderTwoWithoutMean(wMu0dt,Z,ZProb,wnStart,Nn,dNn,NnMidh,NnMidl);
[mC0,mC1,mC2,mC3,Mu0,mC20,MuMean] = CalculateHermitesOrderThreeWithoutMean(wMu0dt,Z,ZProb,wnStart,Nn,dNn,NnMidh,NnMidl)
%Below is hermite decomposition of drift
%wMu0dt(wnStart:Nn)=mC0+mC1(wnStart:Nn).*Z(wnStart:Nn)+mC2(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
%Correction1(wnStart:Nn)=-1/sqrt(1)/(22/7)^1.*mC2(wnStart:Nn).*H2(wnStart:Nn);
w(wnStart:Nn)=wMeanPrev+mC0+mC1(wnStart:Nn).*Z(wnStart:Nn)+mC2(wnStart:Nn).*(Z(wnStart:Nn).^2-1)+mC3(wnStart:Nn).*(Z(wnStart:Nn).^3-3*Z(wnStart:Nn))+ ...
sign(C1(wnStart:Nn).*Z(wnStart:Nn)-Cm1.*ZMean+dw(wnStart:Nn)).* ...
sqrt(abs(sign(C1(wnStart:Nn).*Z(wnStart:Nn)-Cm1.*ZMean).*(C1(wnStart:Nn).*Z(wnStart:Nn)-Cm1.*ZMean).^2+ ...
sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)))+ ...
sign(C2(wnStart:Nn).*(Z(wnStart:Nn).^2-1)-Cm2.*(ZMean.^2-1)+d2w(wnStart:Nn)).* ...
sqrt(abs(sign(C2(wnStart:Nn).*(Z(wnStart:Nn).^2-1)-Cm2.*(ZMean.^2-1)).*(C2(wnStart:Nn).*(Z(wnStart:Nn).^2-1)-Cm2.*(ZMean.^2-1)).^2+ ...
sign(d2w(wnStart:Nn)).*d2w(wnStart:Nn).^2)) + ...
sign(C3(wnStart:Nn).*(Z(wnStart:Nn).^3-3*Z(wnStart:Nn))-Cm3.*(ZMean.^3-3*ZMean)+d3w(wnStart:Nn)).* ...
sqrt(abs((C3(wnStart:Nn).*(Z(wnStart:Nn).^3-3*Z(wnStart:Nn))-Cm3.*(ZMean.^3-3*ZMean)).^2+ ...
sign(d3w(wnStart:Nn)).*d3w(wnStart:Nn).^2))+ ...
Correction1(wnStart:Nn);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Below are instability correction equations. There is instability in the
% SDE representation equation roughly between -.5 SD to +.5 SD in normal terms.
%Instability correction Algorithm takes two ends that are here given as
%lower end=w(NnMidl-2) and upper end=w(NnMidh+3).
%In the block below the interpolation starts from lower end.
%1st equation is
%wll(NnMidl-1)=w(NnMidl-2)+(D1wl+CTheta*D2wl*dNn/1+.5*DeltaD1w)*dNn;
%Second equation is: wll(NnMidl)=wll(NnMidl-1)+(D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
%And so on.
%There is a sixth equation that is not written below(since it takes us to
%the upper interpolation point which is already known) which is
%wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
%Sum of these six equations add up to the distance between upper and lower
%interpolation points.
%In the sum there are six instances of first derivative which is already
%calculated at the boundary, six instances of second derivative which is also
%already calculated at lower boundary and 18.0 instances of DeltaD1w which is
%unknown. This equation calculates the value of one unit of DeltaD1w addition
%DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
% I hope to find a better analytical method to tackle this issue. Friends
% are welcome to play with slightly changed related alternatives to the
%above interpolation.
CTheta=1.0;
D1wl=(11/6*w(NnMidl-2)-3*w(NnMidl-3)+3/2*w(NnMidl-4)-1/3*w(NnMidl-5))/dNn;
D2wl=(35/12*w(NnMidl-2) -26/3*w(NnMidl-3)+ 19/2*w(NnMidl-4)-14/3*w(NnMidl-5)+11/12*w(NnMidl-6))/dNn.^2;
DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
wll(NnMidl-1)=w(NnMidl-2)+ (D1wl+CTheta*D2wl*dNn+.5*DeltaD1w)*dNn;
wll(NnMidl)=wll(NnMidl-1)+ (D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
wll(NnMidh)=wll(NnMidl)+ (D1wl+CTheta*D2wl*dNn+2.5*DeltaD1w)*dNn;
wll(NnMidh+1)=wll(NnMidh)+ (D1wl+CTheta*D2wl*dNn+3.5*DeltaD1w)*dNn;
wll(NnMidh+2)=wll(NnMidh+1)+(D1wl+CTheta*D2wl*dNn+4.5*DeltaD1w)*dNn;
%wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
%The last equation is not used since point w(NnMidh+3)is given but it is used in
%calculation of distance between two interpolation points. This distance is
%used to determine the value of DeltaD1w.
D1wh=(-11/6*w(NnMidh+3)+3*w(NnMidh+4)-3/2*w(NnMidh+5)+1/3*w(NnMidh+6))/dNn;
D2wh=(35/12*w(NnMidh+3) -26/3*w(NnMidh+4)+ 19/2*w(NnMidh+5)-14/3*w(NnMidh+6)+11/12*w(NnMidh+7))/dNn.^2;
DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wh*dNn-6.0*CTheta*D2wh*dNn.^2.0/1.0)/18.0/dNn;
% DeltaD1w=(w(NnMidh+2)-w(NnMidl-1)-4*D1wh*dNn-16.0*D2wh*dNn.^2.0/1.0)/6.0/dNn;
whh(NnMidh+2)=w(NnMidh+3)-(D1wh+CTheta*D2wh*dNn+.5*DeltaD1w)*dNn;
whh(NnMidh+1)=w(NnMidh+2)-(D1wh+CTheta*D2wh*dNn+1.5* DeltaD1w)*dNn;
whh(NnMidh)=whh(NnMidh+1)-(D1wh+CTheta*D2wh*dNn+2.5* DeltaD1w)*dNn;
whh(NnMidl)=whh(NnMidh)-(D1wh+CTheta*D2wh*dNn+3.5* DeltaD1w)*dNn;
whh(NnMidl-1)=whh(NnMidl)-(D1wh+CTheta*D2wh*dNn+4.5* DeltaD1w)*dNn;
%whh(NnMidl-2)=whh(NnMidl-1)-(D1wh+CTheta*D2wh*dNn/1+5.5* DeltaD1w)*dNn;
%
w(NnMidh+2)=.5*whh(NnMidh+2)+.5*wll(NnMidh+2);
w(NnMidh+1)=.5*whh(NnMidh+1)+.5*wll(NnMidh+1);
w(NnMidh)=.5*whh(NnMidh)+.5*wll(NnMidh);
w(NnMidl)=.5*whh(NnMidl)+.5*wll(NnMidl);
w(NnMidl-1)=.5*whh(NnMidl-1)+.5*wll(NnMidl-1);
w1(1:Nn-1)=w(1:Nn-1);
w2(1:Nn-1)=w(2:Nn);
w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
w(w<0)=0.0;
for nn=1:Nn
if(w(nn)<=0)
wnStart=nn+1;
end
end
%%1st order mean correction. We know the mean in original
%%coordinates so I shift the density into original coordinates,
%%apply the mean correction and then transform again into Lamperti
%%coordinates. Algebra solves two equations in two unknowns.
%%Equations are Sum_0^N{(Y_w(wnStart:Nn)-Y0)*W0}=1 and the second
%%equation is Sum_0^N{Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0}=u0
%%Two unknows are Y0 and W0. u0 is known mean.
tt;
u0=theta+(x0-theta)*exp(-kappa*(tt*dt)); %analytic mean of the density
%If you are not using stochastic volatility, replace above with
%true mean otherwise results would become garbage.
Y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
Yn2Pn=sum(Y_w(wnStart:Nn).*Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
YnPn=sum(Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
Pn=1.0;%%Sum(ZProb(1:Nn))
Y0=(Yn2Pn-u0*YnPn)/(YnPn-u0*Pn);
Y0Pn=Y0*Pn;
W0=1/(YnPn-Y0Pn);
YCorrected(wnStart:Nn)=Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0;
wCorrected(wnStart:Nn)=(YCorrected(wnStart:Nn).^(1-gamma))./(1-gamma);
% w(wnStart:Nn)=wCorrected(wnStart:Nn);
%I have disabled the mean correction. The above mean correction is only valid for
%standard mean reverting stochastic volatility type SDEs. To enable mean
%correction please uncomment the above last line in the block.
end
y_w(1:Nn)=0;
y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
Dfy_w(wnStart:Nn)=0;
Dfw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
Dfy_w(nn) = (y_w(nn + 1) - y_w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
Dfw(nn) = (w(nn + 1) - w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
fy_w(1:Nn)=0;
fw(1:Nn)=0;
for nn = wnStart:Nn-1
fy_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
fw(nn) = (normpdf(Z(nn),0, 1))/abs(Dfw(nn));
end
toc
%str=input('Analytic Values calculated; Press a key to start monte carlo');
rng(29079137, 'twister')
paths=100000;
YY(1:paths)=x0; %Original process monte carlo.
Random1(1:paths)=0;
YYMean(1:TtM)=0;
for tt=1:TtM
Random1=randn(size(Random1));
HermiteP1(1,1:paths)=1;
HermiteP1(2,1:paths)=Random1(1:paths);
HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;
YY(1:paths)=YY(1:paths) + ...
(YCoeff0(1,1,2,1).*YY(1:paths).^Fp(1,1,2,1)+ ...
YCoeff0(1,2,1,1).*YY(1:paths).^Fp(1,2,1,1)+ ...
YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1))*dtM + ...
(YCoeff0(1,1,3,1).*YY(1:paths).^Fp(1,1,3,1)+ ...
YCoeff0(1,2,2,1).*YY(1:paths).^Fp(1,2,2,1)+ ...
YCoeff0(2,1,2,1).*YY(1:paths).^Fp(2,1,2,1)+ ...
YCoeff0(1,3,1,1).*YY(1:paths).^Fp(1,3,1,1)+ ...
YCoeff0(2,2,1,1).*YY(1:paths).^Fp(2,2,1,1)+ ...
YCoeff0(3,1,1,1).*YY(1:paths).^Fp(3,1,1,1))*dtM^2 + ...
((YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
(YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+ ...
YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5) .*HermiteP1(2,1:paths) + ...
((YCoeff0(1,1,1,3).*YY(1:paths).^Fp(1,1,1,3) *dtM) + ...
(YCoeff0(1,1,2,3).*YY(1:paths).^Fp(1,1,2,3)+ ...
YCoeff0(1,2,1,3).*YY(1:paths).^Fp(1,2,1,3)+ ...
YCoeff0(2,1,1,3).*YY(1:paths).^Fp(2,1,1,3)).*dtM^2).*HermiteP1(3,1:paths) + ...
((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )).*HermiteP1(4,1:paths) + ...
(YCoeff0(1,1,1,5).*YY(1:paths).^Fp(1,1,1,5)*dtM^2.0).*HermiteP1(5,1:paths);
%Uncomment for fourth order monte carlo
%
% YY(1:paths)=YY(1:paths) + ...
% (YCoeff0(1,1,2,1).*YY(1:paths).^Fp(1,1,2,1)+YCoeff0(1,2,1,1).*YY(1:paths).^Fp(1,2,1,1)+ ...
% YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1))*dtM + ...
% (YCoeff0(1,1,3,1).*YY(1:paths).^Fp(1,1,3,1)+YCoeff0(1,2,2,1).*YY(1:paths).^Fp(1,2,2,1)+ ...
% YCoeff0(2,1,2,1).*YY(1:paths).^Fp(2,1,2,1)+YCoeff0(1,3,1,1).*YY(1:paths).^Fp(1,3,1,1)+ ...
% YCoeff0(2,2,1,1).*YY(1:paths).^Fp(2,2,1,1)+YCoeff0(3,1,1,1).*YY(1:paths).^Fp(3,1,1,1))*dtM^2 + ...
% (YCoeff0(1,1,4,1).*YY(1:paths).^Fp(1,1,4,1)+YCoeff0(1,2,3,1).*YY(1:paths).^Fp(1,2,3,1)+ ...
% YCoeff0(2,1,3,1).*YY(1:paths).^Fp(2,1,3,1)+YCoeff0(1,3,2,1).*YY(1:paths).^Fp(1,3,2,1)+ ...
% YCoeff0(2,2,2,1).*YY(1:paths).^Fp(2,2,2,1)+YCoeff0(3,1,2,1).*YY(1:paths).^Fp(3,1,2,1)+ ...
% YCoeff0(1,4,1,1).*YY(1:paths).^Fp(1,4,1,1)+YCoeff0(2,3,1,1).*YY(1:paths).^Fp(2,3,1,1)+ ...
% YCoeff0(3,2,1,1).*YY(1:paths).^Fp(3,2,1,1)+YCoeff0(4,1,1,1).*YY(1:paths).^Fp(4,1,1,1))*dtM^3 + ...
% (YCoeff0(1,1,5,1).*YY(1:paths).^Fp(1,1,5,1)+YCoeff0(1,2,4,1).*yy(1:paths).^Fp(1,2,4,1)+ ...
% YCoeff0(2,1,4,1).*YY(1:paths).^Fp(2,1,4,1)+YCoeff0(2,1,1,1).*yy(1:paths).^Fp(2,1,1,1)+ ...
% YCoeff0(2,2,3,1).*YY(1:paths).^Fp(2,2,3,1)+ ...
% YCoeff0(3,1,3,1).*YY(1:paths).^Fp(3,1,3,1)+YCoeff0(1,4,2,1).*YY(1:paths).^Fp(1,4,2,1)+ ...
% YCoeff0(2,3,2,1).*YY(1:paths).^Fp(2,3,2,1)+YCoeff0(3,2,2,1).*YY(1:paths).^Fp(3,2,2,1)+ ...
% YCoeff0(4,1,2,1).*YY(1:paths).^Fp(4,1,2,1)+YCoeff0(1,5,1,1).*YY(1:paths).^Fp(1,5,1,1)+ ...
% YCoeff0(2,4,1,1).*YY(1:paths).^Fp(2,4,1,1)+ ...
% YCoeff0(3,3,1,1).*YY(1:paths).^Fp(3,3,1,1)+YCoeff0(4,2,1,1).*YY(1:paths).^Fp(4,2,1,1)+ ...
% YCoeff0(5,1,1,1).*YY(1:paths).^Fp(5,1,1,1))*dtM^4+ ...
% ((YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
% (YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
% YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5+ ...
% (YCoeff0(1,1,3,2).*YY(1:paths).^Fp(1,1,3,2)+YCoeff0(1,2,2,2).*YY(1:paths).^Fp(1,2,2,2)+ ...
% YCoeff0(2,1,2,2).*YY(1:paths).^Fp(2,1,2,2)+YCoeff0(1,3,1,2).*YY(1:paths).^Fp(1,3,1,2)+ ...
% YCoeff0(2,2,1,2).*YY(1:paths).^Fp(2,2,1,2)+YCoeff0(3,1,1,2).*YY(1:paths).^Fp(3,1,1,2)).*dtM^2.5+ ...
% (YCoeff0(1,1,4,2).*YY(1:paths).^Fp(1,1,4,2)+YCoeff0(1,2,3,2).*YY(1:paths).^Fp(1,2,3,2)+ ...
% YCoeff0(2,1,3,2).*YY(1:paths).^Fp(2,1,3,2)+YCoeff0(1,3,2,2).*YY(1:paths).^Fp(1,3,2,2)+ ...
% YCoeff0(2,2,2,2).*YY(1:paths).^Fp(2,2,2,2)+ YCoeff0(3,1,2,2).*YY(1:paths).^Fp(3,1,2,2)+ ...
% YCoeff0(1,4,1,2).*YY(1:paths).^Fp(1,4,1,2)+YCoeff0(2,3,1,2).*YY(1:paths).^Fp(2,3,1,2)+ ...
% YCoeff0(3,2,1,2).*YY(1:paths).^Fp(3,2,1,2)+YCoeff0(4,1,1,2).*YY(1:paths).^Fp(4,1,1,2)).*dtM^3.5) .*HermiteP1(2,1:paths) + ...
% ((YCoeff0(1,1,1,3).*YY(1:paths).^Fp(1,1,1,3) *dtM) + ...
% (YCoeff0(1,1,2,3).*YY(1:paths).^Fp(1,1,2,3)+YCoeff0(1,2,1,3).*YY(1:paths).^Fp(1,2,1,3)+ ...
% YCoeff0(2,1,1,3).*YY(1:paths).^Fp(2,1,1,3)).*dtM^2+ ...
% (YCoeff0(1,1,3,3).*YY(1:paths).^Fp(1,1,3,3)+YCoeff0(1,2,2,3).*YY(1:paths).^Fp(1,2,2,3)+ ...
% YCoeff0(2,1,2,3).*YY(1:paths).^Fp(2,1,2,3) + YCoeff0(1,3,1,3).*YY(1:paths).^Fp(1,3,1,3)+ ...
% YCoeff0(2,2,1,3).*YY(1:paths).^Fp(2,2,1,3)+YCoeff0(3,1,1,3).*YY(1:paths).^Fp(3,1,1,3)).*dtM^3).*HermiteP1(3,1:paths) + ...
% ((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )+ ...
% (YCoeff0(1,1,2,4).*YY(1:paths).^Fp(1,1,2,4)+YCoeff0(1,2,1,4).*YY(1:paths).^Fp(1,2,1,4)+ ...
% YCoeff0(2,1,1,4).*YY(1:paths).^Fp(2,1,1,4))*dtM^2.5).*HermiteP1(4,1:paths) + ...
% (YCoeff0(1,1,1,5).*YY(1:paths).^Fp(1,1,1,5)*dtM^2.0).*HermiteP1(5,1:paths);
%
YYMean(tt)=YYMean(tt)+sum(YY(1:paths))/paths;
end
YY(YY<0)=0;
disp('Original process average from monte carlo');
MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
disp('Original process average from our simulation');
ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates
disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
MaxCutOff=30;
NoOfBins=round(500*gamma^2*1*sigma0);
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');
title('Graphical Comparison of Ito-Hermite Method and Monte Carlo Method Density')
legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
end
Here is the new function for calculation of drift, its derivatives and volatility in terms of hermites. It is named order4 but it is actually order dt^3 upto four hermite polynomials.
function [wMu0dt,dwMu0dtdw,d2wMu0dtdw2,c1,c2,c3,c4] = CalculateDriftAndVolOrder4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt)
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
Fp2=Fp1/(1-gamma);
wMu0dt(wnStart:Nn)=(YqCoeff0(1,1,2,1).*yy(wnStart:Nn).^Fp1(1,1,2,1)+ ...
YqCoeff0(1,2,1,1).*yy(wnStart:Nn).^Fp1(1,2,1,1)+ ...
YqCoeff0(2,1,1,1).*yy(wnStart:Nn).^Fp1(2,1,1,1))*dt + ...
(YqCoeff0(1,1,3,1).*yy(wnStart:Nn).^Fp1(1,1,3,1)+ ...
YqCoeff0(1,2,2,1).*yy(wnStart:Nn).^Fp1(1,2,2,1)+ ...
YqCoeff0(2,1,2,1).*yy(wnStart:Nn).^Fp1(2,1,2,1)+ ...
YqCoeff0(1,3,1,1).*yy(wnStart:Nn).^Fp1(1,3,1,1)+ ...
YqCoeff0(2,2,1,1).*yy(wnStart:Nn).^Fp1(2,2,1,1)+ ...
YqCoeff0(3,1,1,1).*yy(wnStart:Nn).^Fp1(3,1,1,1))*dt^2 + ...
(YqCoeff0(1,1,4,1).*yy(wnStart:Nn).^Fp1(1,1,4,1)+ ...
YqCoeff0(1,2,3,1).*yy(wnStart:Nn).^Fp1(1,2,3,1)+ ...
YqCoeff0(2,1,3,1).*yy(wnStart:Nn).^Fp1(2,1,3,1)+ ...
YqCoeff0(1,3,2,1).*yy(wnStart:Nn).^Fp1(1,3,2,1)+ ...
YqCoeff0(2,2,2,1).*yy(wnStart:Nn).^Fp1(2,2,2,1)+ ...
YqCoeff0(3,1,2,1).*yy(wnStart:Nn).^Fp1(3,1,2,1)+ ...
YqCoeff0(1,4,1,1).*yy(wnStart:Nn).^Fp1(1,4,1,1)+ ...
YqCoeff0(2,3,1,1).*yy(wnStart:Nn).^Fp1(2,3,1,1)+ ...
YqCoeff0(3,2,1,1).*yy(wnStart:Nn).^Fp1(3,2,1,1)+ ...
YqCoeff0(4,1,1,1).*yy(wnStart:Nn).^Fp1(4,1,1,1))*dt^3;
dwMu0dtdw(wnStart:Nn)=(YqCoeff0(1,1,2,1).*Fp1(1,1,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,1,2,1))+ ...
YqCoeff0(1,2,1,1).*Fp1(1,2,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,2,1,1))+ ...
YqCoeff0(2,1,1,1).*Fp1(2,1,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,1,1,1)))*dt + ...
(YqCoeff0(1,1,3,1).*Fp1(1,1,3,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,1,3,1))+ ...
YqCoeff0(1,2,2,1).*Fp1(1,2,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,2,2,1))+ ...
YqCoeff0(2,1,2,1).*Fp1(2,1,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,1,2,1))+ ...
YqCoeff0(1,3,1,1).*Fp1(1,3,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,3,1,1))+ ...
YqCoeff0(2,2,1,1).*Fp1(2,2,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,2,1,1))+ ...
YqCoeff0(3,1,1,1).*Fp1(3,1,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(3,1,1,1)))*dt^2 + ...
(YqCoeff0(1,1,4,1).*Fp1(1,1,4,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,1,4,1))+ ...
YqCoeff0(1,2,3,1).*Fp1(1,2,3,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,2,3,1))+ ...
YqCoeff0(2,1,3,1).*Fp1(2,1,3,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,1,3,1))+ ...
YqCoeff0(1,3,2,1).*Fp1(1,3,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,3,2,1))+ ...
YqCoeff0(2,2,2,1).*Fp1(2,2,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,2,2,1))+ ...
YqCoeff0(3,1,2,1).*Fp1(3,1,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(3,1,2,1))+ ...
YqCoeff0(1,4,1,1).*Fp1(1,4,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,4,1,1))+ ...
YqCoeff0(2,3,1,1).*Fp1(2,3,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,3,1,1))+ ...
YqCoeff0(3,2,1,1).*Fp1(3,2,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(3,2,1,1))+ ...
YqCoeff0(4,1,1,1).*Fp1(4,1,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(4,1,1,1)))*dt^3;
d2wMu0dtdw2(wnStart:Nn)=(1-gamma).*((YqCoeff0(1,1,2,1).*Fp1(1,1,2,1).*(-1+Fp2(1,1,2,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(1,1,2,1))+ ...
YqCoeff0(1,2,1,1).*Fp1(1,2,1,1).*(-1+Fp2(1,2,1,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(1,2,1,1))+ ...
YqCoeff0(2,1,1,1).*Fp1(2,1,1,1).*(-1+Fp2(2,1,1,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(2,1,1,1)))*dt + ...
(YqCoeff0(1,1,3,1).*Fp1(1,1,3,1).*(-1+Fp2(1,1,3,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(1,1,3,1))+ ...
YqCoeff0(1,2,2,1).*Fp1(1,2,2,1).*(-1+Fp2(1,2,2,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(1,2,2,1))+ ...
YqCoeff0(2,1,2,1).*Fp1(2,1,2,1).*(-1+Fp2(2,1,2,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(2,1,2,1))+ ...
YqCoeff0(1,3,1,1).*Fp1(1,3,1,1).*(-1+Fp2(1,3,1,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(1,3,1,1))+ ...
YqCoeff0(2,2,1,1).*Fp1(2,2,1,1).*(-1+Fp2(2,2,1,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(2,2,1,1))+ ...
YqCoeff0(3,1,1,1).*Fp1(3,1,1,1).*(-1+Fp2(3,1,1,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(3,1,1,1)))*dt^2 + ...
(YqCoeff0(1,1,4,1).*Fp1(1,1,4,1).*(-1+Fp2(1,1,4,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(1,1,4,1))+ ...
YqCoeff0(1,2,3,1).*Fp1(1,2,3,1).*(-1+Fp2(1,2,3,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(1,2,3,1))+ ...
YqCoeff0(2,1,3,1).*Fp1(2,1,3,1).*(-1+Fp2(2,1,3,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(2,1,3,1))+ ...
YqCoeff0(1,3,2,1).*Fp1(1,3,2,1).*(-1+Fp2(1,3,2,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(1,3,2,1))+ ...
YqCoeff0(2,2,2,1).*Fp1(2,2,2,1).*(-1+Fp2(2,2,2,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(2,2,2,1))+ ...
YqCoeff0(3,1,2,1).*Fp1(3,1,2,1).*(-1+Fp2(3,1,2,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(3,1,2,1))+ ...
YqCoeff0(1,4,1,1).*Fp1(1,4,1,1).*(-1+Fp2(1,4,1,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(1,4,1,1))+ ...
YqCoeff0(2,3,1,1).*Fp1(2,3,1,1).*(-1+Fp2(2,3,1,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(2,3,1,1))+ ...
YqCoeff0(3,2,1,1).*Fp1(3,2,1,1).*(-1+Fp2(3,2,1,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(3,2,1,1))+ ...
YqCoeff0(4,1,1,1).*Fp1(4,1,1,1).*(-1+Fp2(4,1,1,1)).*((1-gamma).*w(wnStart:Nn)).^(-2+Fp2(4,1,1,1)))*dt^3);
c1(wnStart:Nn)=YqCoeff0(1,1,1,2).*yy(wnStart:Nn).^Fp1(1,1,1,2).*sqrt(dt)+ ...
(YqCoeff0(1,1,2,2).*yy(wnStart:Nn).^Fp1(1,1,2,2)+YqCoeff0(1,2,1,2).*yy(wnStart:Nn).^Fp1(1,2,1,2)+ ...
YqCoeff0(2,1,1,2).*yy(wnStart:Nn).^Fp1(2,1,1,2)).*dt^1.5+ ...
(YqCoeff0(1,1,3,2).*yy(wnStart:Nn).^Fp1(1,1,3,2)+YqCoeff0(1,2,2,2).*yy(wnStart:Nn).^Fp1(1,2,2,2)+ ...
YqCoeff0(2,1,2,2).*yy(wnStart:Nn).^Fp1(2,1,2,2)+YqCoeff0(1,3,1,2).*yy(wnStart:Nn).^Fp1(1,3,1,2)+ ...
YqCoeff0(2,2,1,2).*yy(wnStart:Nn).^Fp1(2,2,1,2)+YqCoeff0(3,1,1,2).*yy(wnStart:Nn).^Fp1(3,1,1,2)).*dt^2.5;
c2(wnStart:Nn)=YqCoeff0(1,1,1,3).*yy(wnStart:Nn).^Fp1(1,1,1,3) *dt + ...
(YqCoeff0(1,1,2,3).*yy(wnStart:Nn).^Fp1(1,1,2,3)+YqCoeff0(1,2,1,3).*yy(wnStart:Nn).^Fp1(1,2,1,3)+ ...
YqCoeff0(2,1,1,3).*yy(wnStart:Nn).^Fp1(2,1,1,3)).*dt^2+ ...
(YqCoeff0(1,1,3,3).*yy(wnStart:Nn).^Fp1(1,1,3,3)+YqCoeff0(1,2,2,3).*yy(wnStart:Nn).^Fp1(1,2,2,3)+ ...
YqCoeff0(2,1,2,3).*yy(wnStart:Nn).^Fp1(2,1,2,3) + YqCoeff0(1,3,1,3).*yy(wnStart:Nn).^Fp1(1,3,1,3)+ ...
YqCoeff0(2,2,1,3).*yy(wnStart:Nn).^Fp1(2,2,1,3)+YqCoeff0(3,1,1,3).*yy(wnStart:Nn).^Fp1(3,1,1,3)).*dt^3;
c3(wnStart:Nn)=YqCoeff0(1,1,1,4).*yy(wnStart:Nn).^Fp1(1,1,1,4)*dt^1.5 + ...
(YqCoeff0(1,1,2,4).*yy(wnStart:Nn).^Fp1(1,1,2,4)+YqCoeff0(1,2,1,4).*yy(wnStart:Nn).^Fp1(1,2,1,4)+ ...
YqCoeff0(2,1,1,4).*yy(wnStart:Nn).^Fp1(2,1,1,4))*dt^2.5;
c4(wnStart:Nn)=YqCoeff0(1,1,1,5).*yy(wnStart:Nn).^Fp1(1,1,1,5)*dt^2.0;
end
Here is another re-worked function that you will need.
function [w0,Sigma00,Sigma11,Sigma22,Sigma33,Sigma1,Sigma2,Sigma3,Sigma0] = CalculateHermiteRepresentationOfDensityOrderThreeNew(w,Z,wnStart,Nn,dNn,NnMidh,NnMidl)
%This program find hermite representation of the SDE variable w to
%tjird Order as a function of three hermite polynomials.
%This program is not perfectly precise and has very small errors so you
%will have to be careful using it everywhere.
%Here wnStart is the starting grid point of the SDE vector
%mostly wnStart=1.
%Please note we are representing SDE Variable w as a function of hermite
%polynomials. It is not the density which can however be found by change of variable
%derivative of w with respect to Z using standard gaussian variable.
%Nn is the index of last density grid point.
%So w(wnStart:Nn) represents the SDE density variable grid.
% Z(wnStart:Nn) represents the associated value of standard Gaussian on
% the grid. w and Z can be associated by matching appropriate CDF points
% on the grid.
%I suppose the grid is uniform in Z and dNn is the spacing between and
% Z(n) and Z(n+1) two points on the grid.
Z0=0;
w0= InterpolateOrderN6(6,Z0,Z(NnMidh+2),Z(NnMidh+1),Z(NnMidh),Z(NnMidl),Z(NnMidl-1),Z(NnMidl-2),w(NnMidh+2),w(NnMidh+1),w(NnMidh),w(NnMidl),w(NnMidl-1),w(NnMidl-2));
%Below Calculate the vol associated with each point on the grid
for nn=wnStart:Nn
Sigma0(nn)=(w(nn)-w0)./(Z(nn)-Z0);
end
%Below calculate the associated Sigma00 which is a single volatility associated with
%the first orthogonal hermite polynomial.
Sigma00= InterpolateOrderN(4,Z0,Z(NnMidh+1),Z(NnMidh),Z(NnMidl),Z(NnMidl-1),0,Sigma0(NnMidh+1),Sigma0(NnMidh),Sigma0(NnMidl),Sigma0(NnMidl-1),0);
%Below calculate the associated first derivative of Sigma0 with respect to
%Z on the grid.
dSigma0dZ(wnStart)=(-3*Sigma0(wnStart)+4*Sigma0(wnStart+1)-1*Sigma0(wnStart+2))/(2*dNn);
dSigma0dZ(wnStart+1:Nn-1)=(-1*Sigma0(wnStart:Nn-2)+1*Sigma0(wnStart+2:Nn))/(2*dNn);
dSigma0dZ(Nn)=(3*Sigma0(Nn)-4*Sigma0(Nn-1)+1*Sigma0(Nn-2))/(2*dNn);
dSigma0dZ(wnStart+2:Nn-2)= (1*Sigma0(wnStart:Nn-4)-8*Sigma0(wnStart+1:Nn-3)+0*Sigma0(wnStart+2:Nn-2)+8*Sigma0(wnStart+3:Nn-1)-1*Sigma0(wnStart+4:Nn))/(12*dNn);
dSigma0dZ0 = (-9*Sigma0(NnMidl-2)+125*Sigma0(NnMidl-1)-2250*Sigma0(NnMidl)+2250*Sigma0(NnMidh)-125*Sigma0(NnMidh+1)+9*Sigma0(NnMidh+2))/(3840*(.5*dNn));
%Following is the operation related to equation(4) in post 860.
SSigma1(NnMidh)=(.5*dSigma0dZ0+.5*dSigma0dZ(NnMidh))*dNn/2;
for nn=NnMidh+1:Nn
SSigma1(nn)=SSigma1(nn-1)+(.5*dSigma0dZ(nn-1)+.5*dSigma0dZ(nn))*dNn;
end
SSigma1(NnMidl)=-(.5*dSigma0dZ0+.5*dSigma0dZ(NnMidl))*dNn/2;
for nn=NnMidl-1:-1:wnStart
SSigma1(nn)=SSigma1(nn+1)-(.5*dSigma0dZ(nn+1)+.5*dSigma0dZ(nn))*dNn;
end
Sigma1(wnStart:Nn)=SSigma1(wnStart:Nn)./Z(wnStart:Nn);
Sigma11= InterpolateOrderN6(6,Z0,Z(NnMidh+2),Z(NnMidh+1),Z(NnMidh),Z(NnMidl),Z(NnMidl-1),Z(NnMidl-2),Sigma1(NnMidh+2),Sigma1(NnMidh+1),Sigma1(NnMidh),Sigma1(NnMidl),Sigma1(NnMidl-1),Sigma1(NnMidl-2));
%Now below we represent the SDE variable w in terms of first two hermite[/font][/size]
%polynomials and this associated value of w produced by hermite polynomials
%coefficients is called w1.
% w1(wnStart:Nn)=wMean-Sigma00*ZMean-Sigma1(wnStart:Nn).*(ZMean.^2-1)+ ...
% Sigma00*Z(wnStart:Nn)+ ...
% Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
w0(wnStart:Nn)=w0+Sigma0(wnStart:Nn).*Z(wnStart:Nn);
w1=w0;
w2=w0;
%w1(wnStart:Nn)=w0+ Sigma00*Z(wnStart:Nn)+ ...
% Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
%w2(wnStart:Nn)=w0+Sigma00*Z(wnStart:Nn)+ ...
% Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
dSigma1dZ(wnStart)=(-3*Sigma1(wnStart)+4*Sigma1(wnStart+1)-1*Sigma1(wnStart+2))/(2*dNn);
dSigma1dZ(wnStart+1:Nn-1)=(-1*Sigma1(wnStart:Nn-2)+1*Sigma1(wnStart+2:Nn))/(2*dNn);
dSigma1dZ(Nn)=(3*Sigma1(Nn)-4*Sigma1(Nn-1)+1*Sigma1(Nn-2))/(2*dNn);
dSigma1dZ0 = (-9*Sigma1(NnMidl-2)+125*Sigma1(NnMidl-1)-2250*Sigma1(NnMidl)+2250*Sigma1(NnMidh)-125*Sigma1(NnMidh+1)+9*Sigma1(NnMidh+2))/(3840*(.5*dNn));
SSigma2(NnMidh)=(.5*dSigma1dZ0+.5*dSigma1dZ(NnMidh))*dNn/2;
for nn=NnMidh+1:Nn
SSigma2(nn)=SSigma2(nn-1)+(.5*dSigma1dZ(nn-1)+.5*dSigma1dZ(nn))*dNn;
end
SSigma2(NnMidl)=-(.5*dSigma1dZ0+.5*dSigma1dZ(NnMidl))*dNn/2;
for nn=NnMidl-1:-1:wnStart
SSigma2(nn)=SSigma2(nn+1)-(.5*dSigma1dZ(nn+1)+.5*dSigma1dZ(nn))*dNn;
end
Sigma2(wnStart:Nn)=SSigma2(wnStart:Nn)./Z(wnStart:Nn);
Sigma22= InterpolateOrderN6(6,Z0,Z(NnMidh+2),Z(NnMidh+1),Z(NnMidh),Z(NnMidl),Z(NnMidl-1),Z(NnMidl-2),Sigma2(NnMidh+2),Sigma2(NnMidh+1),Sigma2(NnMidh),Sigma2(NnMidl),Sigma2(NnMidl-1),Sigma2(NnMidl-2));
dSigma2dZ(wnStart)=(-3*Sigma2(wnStart)+4*Sigma2(wnStart+1)-1*Sigma2(wnStart+2))/(2*dNn);
dSigma2dZ(wnStart+1:Nn-1)=(-1*Sigma2(wnStart:Nn-2)+1*Sigma2(wnStart+2:Nn))/(2*dNn);
dSigma2dZ(Nn)=(3*Sigma2(Nn)-4*Sigma2(Nn-1)+1*Sigma2(Nn-2))/(2*dNn);
dSigma2dZ0 = (-9*Sigma2(NnMidl-2)+125*Sigma2(NnMidl-1)-2250*Sigma2(NnMidl)+2250*Sigma2(NnMidh)-125*Sigma2(NnMidh+1)+9*Sigma2(NnMidh+2))/(3840*(.5*dNn));
SSigma3(NnMidh)=(.5*dSigma2dZ0+.5*dSigma2dZ(NnMidh))*dNn/2;
for nn=NnMidh+1:Nn
SSigma3(nn)=SSigma3(nn-1)+(.5*dSigma2dZ(nn-1)+.5*dSigma2dZ(nn))*dNn;
end
SSigma3(NnMidl)=-(.5*dSigma2dZ0+.5*dSigma2dZ(NnMidl))*dNn/2;
for nn=NnMidl-1:-1:wnStart
SSigma3(nn)=SSigma3(nn+1)-(.5*dSigma2dZ(nn+1)+.5*dSigma2dZ(nn))*dNn;
end
Sigma3(wnStart:Nn)=SSigma3(wnStart:Nn)./Z(wnStart:Nn);
Sigma33= InterpolateOrderN6(6,Z0,Z(NnMidh+2),Z(NnMidh+1),Z(NnMidh),Z(NnMidl),Z(NnMidl-1),Z(NnMidl-2),Sigma3(NnMidh+2),Sigma3(NnMidh+1),Sigma3(NnMidh),Sigma3(NnMidl),Sigma3(NnMidl-1),Sigma3(NnMidl-2));
end
Here is another function
function [C0,C1,C2,C3,w0,C20,wMean,ZMean,nnMean,Cm1,Cm2,Cm3] = CalculateHermitesOrderThree(w,Z,Prob,wnStart,Nn,dNn,NnMidh,NnMidl)
%-------------------------------------------------------------------------
%Below hermite coefficients calculation starts. InterpolateOrderN is a
%program that interpolates the values of w0 at Z=0 since in our set up grid
%points are located at Z(NnMidl)=-.1 and Z(NnMidh)=.1 and there are no grid
%points at Z=0.
Z0=0;
%w0= InterpolateOrderN(4,Z0,Z(NnMidh),Z(NnMidh+1),Z(NnMidl),Z(NnMidl-1),0,w(NnMidh),w(NnMidh+1),w(NnMidl),w(NnMidl-1),0);
w0= InterpolateOrderN6(6,Z0,Z(NnMidh+2),Z(NnMidh+1),Z(NnMidh),Z(NnMidl),Z(NnMidl-1),Z(NnMidl-2),w(NnMidh+2),w(NnMidh+1),w(NnMidh),w(NnMidl),w(NnMidl-1),w(NnMidl-2));
%dwdZ0=(w(NnMidh)-w(NnMidl))./dNn;
%dwdZ0= (1*w(NnMidl-1)-27*w(NnMidl)+27*w(NnMidh)-1*w(NnMidh+1))/(48*(.5*dNn));
dwdZ0 = (-9*w(NnMidl-2)+125*w(NnMidl-1)-2250*w(NnMidl)+2250*w(NnMidh)-125*w(NnMidh+1)+9*w(NnMidh+2))/(3840*(.5*dNn));
%d2wdZ20 = (-1*w(NnMidl-1)+81*w(NnMidl)-160*w0+81*w(NnMidh)-1*w(NnMidh+1))/(72*(.5*dNn)^2);
%d2wdZ20 = (-1*w(NnMidl-1)+81*w(NnMidl)-160*w0+81*w(NnMidh)-1*w(NnMidh+1))/(72*(.5*dNn)^2)
d2wdZ20 = (-5*w(NnMidl-2)+39*w(NnMidl-1)-34*w(NnMidl)-34*w(NnMidh)+39*w(NnMidh+1)-5*w(NnMidh+2))/(192*(.5*dNn)^2);
%d4wdZ40 = (-3*w(NnMidl-2)+65*w(NnMidl-1)-510*w(NnMidl)+896*w0-510*w(NnMidh)+65*w(NnMidh+1)-3*w(NnMidh+2))/(240*(.5*dNn)^2)
d4wdZ40 = (1*w(NnMidl-2)-3*w(NnMidl-1)+2*w(NnMidl)+2*w(NnMidh)-3*w(NnMidh+1)+1*w(NnMidh+2))/(32*(.5*dNn)^2);
%str=input('Look at derivatives at zero');
%Calculate first hermite polynomial coefficients C1 below.
dwdZ(wnStart)=(-3*w(wnStart)+4*w(wnStart+1)-1*w(wnStart+2))/(2*dNn);
dwdZ(wnStart+1:Nn-1)=(-1*w(wnStart:Nn-2)+1*w(wnStart+2:Nn))/(2*dNn);
dwdZ(Nn)=(3*w(Nn)-4*w(Nn-1)+1*w(Nn-2))/(2*dNn);
dwdZ(wnStart+2:Nn-2)= (1*w(wnStart:Nn-4)-8*w(wnStart+1:Nn-3)+0*w(wnStart+2:Nn-2)+8*w(wnStart+3:Nn-1)-1*w(wnStart+4:Nn))/(12*dNn);
CSum(wnStart:Nn)=0.0;
C1(wnStart:Nn)=0;
C1(NnMidh)=.5*dwdZ0+.5*dwdZ(NnMidh);
CSum(NnMidh)=C1(NnMidh).*Z(NnMidh);
for nn=NnMidh+1:Nn
CSum(nn)=(C1(nn-1).*Z(nn-1)+(.5*dwdZ(nn-1)+.5*dwdZ(nn))*dNn);
C1(nn)=CSum(nn)./Z(nn);
end
C1(NnMidl)=.5*dwdZ0+.5*dwdZ(NnMidl);
CSum(NnMidl)=C1(NnMidl).*Z(NnMidl);
for nn=NnMidl-1:-1:wnStart
CSum(nn)=(C1(nn+1).*Z(nn+1)-(.5*dwdZ(nn+1)+.5*dwdZ(nn))*dNn);
C1(nn)=CSum(nn)./Z(nn);
end
%Calculate second hermite coefficients below. Here C20 are coefficients of
%second hermite polynomial at Z=0
C20=.5*d2wdZ20;
C40=1*1/24*d4wdZ40;
C0=w0+C20-3*C40;
[dwdZ,d2wdZ2,d3wdZ3,d4wdZ4] = First4Derivatives2ndOrderEqSpaced(wnStart,Nn,dNn,w,Z);
d2wdZ2(wnStart+2:Nn-2)= (-1*w(wnStart:Nn-4)+16*w(wnStart+1:Nn-3)-30*w(wnStart+2:Nn-2)+16*w(wnStart+3:Nn-1)-1*w(wnStart+4:Nn))/(12*dNn^2);
C2Sum1(wnStart:Nn)=0.0;
C2(wnStart:Nn)=0;
C2Sum1(NnMidh)=d2wdZ20*dNn/4+d2wdZ2(NnMidh)*dNn/4;
C2(NnMidh)=(0*dNn/4 + C2Sum1(NnMidh)*dNn/4)/Z(NnMidh).^2;
for nn=NnMidh+1:Nn
C2Sum1(nn)=C2Sum1(nn-1)+(d2wdZ2(nn-1)+d2wdZ2(nn))*dNn/2;
C2(nn)=(C2Sum1(nn-1)*dNn/2+C2Sum1(nn)*dNn/2)./(Z(nn).^2);
end
C2Sum1(NnMidl)=d2wdZ20*dNn/4+d2wdZ2(NnMidl)*dNn/4;
C2(NnMidl)=(0*dNn/4 + C2Sum1(NnMidl)*dNn/4)/Z(NnMidl)^2;
for nn=NnMidl-1:-1:wnStart
C2Sum1(nn)=C2Sum1(nn+1)+(d2wdZ2(nn+1)+d2wdZ2(nn))*dNn/2;
C2(nn)=(C2Sum1(nn+1)*dNn/2+C2Sum1(nn)*dNn/2)./(Z(nn).^2);
end
for nn=NnMidh:Nn
C3(nn)=((w(nn)-w0)-C1(nn).*Z(nn)-C2(nn).*(Z(nn)^2-1)-C20+3*C40)./(Z(nn).^3-3*Z(nn));
end
for nn=NnMidl:-1:wnStart
C3(nn)=((w(nn)-w0)-C1(nn).*Z(nn)-C2(nn).*(Z(nn)^2-1)-C20+3*C40)./(Z(nn).^3-3*Z(nn));
end
w2(wnStart:Nn)=C0+C1(wnStart:Nn).*Z(wnStart:Nn)+C2(wnStart:Nn).*(Z(wnStart:Nn).^2-1)+C3(wnStart:Nn).*(Z(wnStart:Nn).^3-3*Z(wnStart:Nn));
%Mean is calculated and returned but it has nothing to do with calculation
%of hermite coefficients. It is for future use.
%--------------------------------------------------------------------
wMean=sum(w(wnStart:Nn).*Prob(wnStart:Nn)); %Calculate the mean. This is in general
%nnMean=NnMidl;
%for nn=wnStart:Nn-1
% if( (w(nn+1)>wMean) && (w(nn)<=wMean))
% %ZMean=Z(nn)+(wMean-w(nn))/(w(nn+1)-w(nn))*(Z(nn+1)-Z(nn));
% nnMean=nn;
% end
%end
MeanFoundFlag=0;
nnMean=NnMidl;
for nn=NnMidl:Nn-1
if( (w(nn+1)>wMean) && (w(nn)<=wMean)&& (MeanFoundFlag==0))
%ZMean=Z(nn)+(wMean-w(nn))/(w(nn+1)-w(nn))*(Z(nn+1)-Z(nn));
nnMean=nn;
MeanFoundFlag=1;
end
end
for nn=NnMidl-1:-1:wnStart
if( (w(nn+1)>wMean) && (w(nn)<=wMean)&& (MeanFoundFlag==0))
%ZMean=Z(nn)+(wMean-w(nn))/(w(nn+1)-w(nn))*(Z(nn+1)-Z(nn));
nnMean=nn;
MeanFoundFlag=1;
end
end
if(abs(wMean-w(nnMean))<abs(wMean-w(nnMean+1)))
ZMean= InterpolateOrderN(5,wMean,w(nnMean-2),w(nnMean-1),w(nnMean),w(nnMean+1),w(nnMean+2),Z(nnMean-2),Z(nnMean-1),Z(nnMean),Z(nnMean+1),Z(nnMean+2));
CSum=InterpolateOrderN(5,wMean,w(nnMean-2),w(nnMean-1),w(nnMean),w(nnMean+1),w(nnMean+2),CSum(nnMean-2),CSum(nnMean-1),CSum(nnMean),CSum(nnMean+1),CSum(nnMean+2));
%C1=InterpolateOrderN(5,wMean,w(nnMean-2),w(nnMean-1),w(nnMean),w(nnMean+1),w(nnMean+2),C1(nnMean-2),C1(nnMean-1),C1(nnMean),C1(nnMean+1),C1(nnMean+2));
Cm2=InterpolateOrderN(5,wMean,w(nnMean-2),w(nnMean-1),w(nnMean),w(nnMean+1),w(nnMean+2),C2(nnMean-2),C2(nnMean-1),C2(nnMean),C2(nnMean+1),C2(nnMean+2));
end
if(abs(wMean-w(nnMean))>=abs(wMean-w(nnMean+1)))
ZMean= InterpolateOrderN(5,wMean,w(nnMean-1),w(nnMean),w(nnMean+1),w(nnMean+2),w(nnMean+3),Z(nnMean-1),Z(nnMean),Z(nnMean+1),Z(nnMean+2),Z(nnMean+3));
CSum=InterpolateOrderN(5,wMean,w(nnMean-1),w(nnMean),w(nnMean+1),w(nnMean+2),w(nnMean+3),CSum(nnMean-1),CSum(nnMean),CSum(nnMean+1),CSum(nnMean+2),CSum(nnMean+3));
%Cm1=InterpolateOrderN(5,wMean,w(nnMean-1),w(nnMean),w(nnMean+1),w(nnMean+2),w(nnMean+3),C1(nnMean-1),C1(nnMean),C1(nnMean+1),C1(nnMean+2),C1(nnMean+3));
Cm2=InterpolateOrderN(5,wMean,w(nnMean-1),w(nnMean),w(nnMean+1),w(nnMean+2),w(nnMean+3),C2(nnMean-1),C2(nnMean),C2(nnMean+1),C2(nnMean+2),C2(nnMean+3));
end
%ZMean=InterpolateOrderN(4,wMean,w(nnMean-1),w(nnMean),w(nnMean+1),w(nnMean+2),0,Z(nnMean-1),Z(nnMean),Z(nnMean+1),Z(nnMean+2),0);
%CSum=InterpolateOrderN(4,wMean,w(nnMean-1),w(nnMean),w(nnMean+1),w(nnMean+2),0,CSum(nnMean-1),CSum(nnMean),CSum(nnMean+1),CSum(nnMean+2),0);
Cm1=CSum./ZMean;
% if ((ZMean>=0)&&(ZMean<dNn/4))
% Cm1=dwdZ0;
% end
% if ((ZMean>dNn/4)&&(ZMean<dNn/2))
% Cm1=(dwdZ0*dNn/4+(ZMean-dNn/4)*dwdZ(NnMidh))/ZMean;
% end
%
% if(ZMean>Z(NnMidh))
% if(ZMean<Z(nnMean)+dNn/2)
% Cm1=(C1(nnMean).*Z(nnMean)+dwdZ(nnMean).*(ZMean-Z(nnMean)))./ZMean;
% end
% if(ZMean>Z(nnMean)+dNn/2)
% Cm1=(C1(nnMean).*Z(nnMean)+dwdZ(nnMean).*dNn/2+dwdZ(nnMean+1).*(ZMean-Z(nnMean)-dNn/2))./ZMean;
% end
% end
%
% if ((ZMean<0)&&(ZMean>-dNn/4))
% Cm1=dwdZ0;
% end
% if ((ZMean<-dNn/4)&&(ZMean>-dNn/2))
% Cm1=(dwdZ0*-dNn/4+(ZMean+dNn/4)*dwdZ(NnMidl))/ZMean;
% end
%
% if(ZMean<Z(NnMidl))
% if(ZMean>Z(nnMean+1)-dNn/2)
% Cm1=(C1(nnMean+1).*Z(nnMean+1)+dwdZ(nnMean).*(ZMean-Z(nnMean+1)))./ZMean;
% end
% if(ZMean<Z(nnMean)-dNn/2)
% Cm1=(C1(nnMean+1).*Z(nnMean+1)+dwdZ(nnMean+1).*-dNn/2+dwdZ(nnMean).*(ZMean-Z(nnMean+1)+dNn/2))./ZMean;
% end
% end
Cm3=((wMean-w0)-Cm1.*ZMean-Cm2.*(ZMean.^2-1)-C20+3*C40)./(ZMean.^3-3*ZMean);
%Below we plot a graphical comparison of the original SDE variable w and
%the variable w1 reproduced from the coefficients of the hermite
%polynomials
% plot(Z(wnStart:Nn),w(wnStart:Nn),'g',Z(wnStart:Nn),w2(wnStart:Nn),'r')
%Sigma00
C1
C2
C0
% str=input('Look at the graphical comparison of w and w1 as a function of standard gaussian Z');
end
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal