Serving the Quantitative Finance Community

Amin
Topic Author
Posts: 2894
Joined: July 14th, 2002, 3:00 am

### Re: Breakthrough in the theory of stochastic differential equations and their simulation

Friends, I continued coding the solution of 2D SV fokker-planck equation last night. I am sure I will complete the coding part today and then post the matlab program in another two days. I have taken the same pair of SDEs for 2D Fokker-Planck that I used for 2D SV monte carlo program that I posted a few days ago. I have first started with FPE in bessel coordinates but later I will post a new program in original coordinates as well.

Amin
Topic Author
Posts: 2894
Joined: July 14th, 2002, 3:00 am

### Re: Breakthrough in the theory of stochastic differential equations and their simulation

I will continue to write about torture and victimization experiences I face and write everyday about it in off-topic but I would try to stay away from saying anything about mind control apparatus and I really hope that good people will intervene to end my mind control and that of a large number of other intelligent victims and let them start their lives again with their human dignity.

Amin
Topic Author
Posts: 2894
Joined: July 14th, 2002, 3:00 am

### Re: Breakthrough in the theory of stochastic differential equations and their simulation

Friends, I have been able to run my first prototype program with zero correlations but it currently seems to have a small mismatch with monte carlo but I am quite close. I hope that I would be able to fix any errors and post the program in a few days.

Amin
Topic Author
Posts: 2894
Joined: July 14th, 2002, 3:00 am

### Re: Breakthrough in the theory of stochastic differential equations and their simulation

Friends my program with SV SDEs without correlations is working well for many SDEs I initially tried. The program was already fine but the way I was averaging density of assets to integrate over SV density to find marginal density of asset had some problems since I was integrating along Z-lines and not along the asset value and that was a wrong approach and changed the density shape. In another 2-3 days I would post this algorithm since I still have to improve the density integration part somewhat, try to check the algorithm on wide range of scenarios and then write comments but the good news is that the program without correlations is working perfectly fine. I have taken both SDEs in bessel coordinates initially and then will switch to original coordinates later. I would post the new program in another 2-3 days.

Amin
Topic Author
Posts: 2894
Joined: July 14th, 2002, 3:00 am

### Re: Breakthrough in the theory of stochastic differential equations and their simulation

Friends, here are a few sample graphs for the asset price density in Stochastic volatility set up from the new program. You can see some very minor discrepancy in the graphs. I will be changing my interpolation routine tomorrow and hope that this discrepancy would end or at least decrease even further. Here are the graphs. You can notice fat tails of asset price distributions due to SV. I will further work on improving the programs and post them here in 2-3 days.

Amin
Topic Author
Posts: 2894
Joined: July 14th, 2002, 3:00 am

### Re: Breakthrough in the theory of stochastic differential equations and their simulation

Friends, I spent a lot of time figuring out how to solve the 2D SV PDE with correlations but could not succeed. I was thinking if I could easily solve the correlated problem, I will post both the programs together but it seems to take some more time therefore I have decided to post the work on non-correlated SV diffusions that I have already done. But it will still take at least a day since I have to write one or two functions to convert asset prices from two dim coordinates to 1 dim coordinates. I have written a function already that I had used to make the graphs I posted on internet but that function is not robust and I have better ideas now to write a simple function. It will not be late now since I will first complete and post the non-correlated SV program and then do any further work.

Amin
Topic Author
Posts: 2894
Joined: July 14th, 2002, 3:00 am

### Re: Breakthrough in the theory of stochastic differential equations and their simulation

After posting the non-correlated SV solution of PDE, I plan to work on a monte carlo simulation program that combines stochastic volatility, correlation between asset and SV diffusions, and enhanced quadratic volatility type structure for asset volatility so that we can have interesting shapes for implied volatility like bimodal implied volatility densities and many other. I will add enhanced quadratic volatility structure in asset SDE on top of what I have already done in previous SV monte carlo simulation program that I posted about ten days ago. I hope the program will be useful for friends who want more interesting shapes in implied volatility.

Amin
Topic Author
Posts: 2894
Joined: July 14th, 2002, 3:00 am

### Re: Breakthrough in the theory of stochastic differential equations and their simulation

Friends here is the PDE solution of non-correlated set of SV SDEs. I will post more detailed notes later.

Here is the main program.
.
function [] = SVSDEAnalyticAndMonteCarlo2DWmtNew()

% dX(t)=(a X(t)^alpha1 + b X(t)^alpha2) dt + sqrt(1-rho^2) sigma1 V(t)^gammaV X(t)^gammaX dZ1(t) + rho sigma1 V(t)^gammaV X(t)^gammaX dZ2(t)
% dV(t)=(mu1 V(t)^beta1 + mu2 V(t)^beta2)*dt + sigma0 V(t)^gamma dZ2(t)

%%%%%%%%%%%%%
X0=1.0;
alpha1=0;
alpha2=1;
%kappaX=2.0;
% =.5;
kappaX=0.0;%If you give drift to density, you may have to decrease the step size.
thetaX=0.750;
a=kappaX*thetaX;
b=-kappaX;
%a=0;
%b=0;
rho=-.7*0;
sigma1=.5;
gammaV=.5;
gammaX=.75;
%%%%%%%%%%%%%
V0=.850;
thetaV=.25;
kappaV=2.0;
mu1=kappaV*thetaV;
mu2=-kappaV;
beta1=0;
beta2=1;
gamma=.95;
sigma0=.75;
%%%%%%%%%%%%

dtM=.03125;
TtM=32;%%%32 steps mean one year.

dt=dtM/16;
Tt=TtM*16;
T=TtM*dtM;

dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=50;  % 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;

dMm=.25/1;   % Normal density subdivisions width. would change with number of subdivisions
Mm=40;  % No of normal density subdivisions
MmMidl=20;%One half density Subdivision left from mid of normal density(low)
MmMidh=21;%One half density subdivision right from the mid of normal density(high)
MmMid=4.0;

w(1:Nn)=V0^(1-gamma)/(1-gamma);
y(1:Nn,1:Mm)=X0^(1-gammaX)/(1-gammaX);

Z(1:Nn)=(((1:Nn)-5.5)*dNn-NnMid);

Z0(1:Mm)=(((1:Mm)-4.5)*dMm-MmMid);
Z
Z0
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);
ZProb(2:Nn-1)=normcdf(.5*Z(2:Nn-1)+.5*Z(3:Nn),0,1)-normcdf(.5*Z(2:Nn-1)+.5*Z(1:Nn-2),0,1);

ZCProb(1:Nn)=normcdf(Z(1:Nn));

Z0Prob(1)=normcdf(.5*Z0(1)+.5*Z0(2),0,1)-normcdf(.5*Z0(1)+.5*Z0(2)-dMm,0,1);
Z0Prob(Mm)=normcdf(.5*Z0(Mm)+.5*Z0(Mm-1)+dMm,0,1)-normcdf(.5*Z0(Mm)+.5*Z0(Mm-1),0,1);
Z0Prob(2:Mm-1)=normcdf(.5*Z0(2:Mm-1)+.5*Z0(3:Mm),0,1)-normcdf(.5*Z0(2:Mm-1)+.5*Z0(1:Mm-2),0,1);

Z0CProb(1:Mm)=normcdf(Z0(1:Mm));

%%%%%%%%%%%%%%%%%%%%%%%%5
wnStart=1;
ynStart=1;
dwdZ(wnStart:Nn)=0.0;
tic
for tt=1:Tt

V(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
X(wnStart:Nn,ynStart:Mm)=((1-gammaX)*y(wnStart:Nn,ynStart:Mm)).^(1/(1-gammaX));

wMu0dt(wnStart:Nn)= (mu1*V(wnStart:Nn).^(beta1-gamma)+mu2*V(wnStart:Nn).^(beta2-gamma) - ...
.5*gamma*sigma0^2*V(wnStart:Nn).^(gamma-1))*dt + ...
(mu1*(beta1-gamma).*V(wnStart:Nn).^(beta1-gamma-1)+mu2*(beta2-gamma)*V(wnStart:Nn).^(beta2-gamma-1) - ...
.5*gamma*sigma0^2*(gamma-1)*V(wnStart:Nn).^(gamma-2)).* ...
(mu1*V(wnStart:Nn).^(beta1-gamma)+mu2*V(wnStart:Nn).^(beta2-gamma))*dt^2/2+ ...
(mu1*(beta1-gamma).*(beta1-gamma-1).*V(wnStart:Nn).^(beta1-gamma-2)+mu2*(beta2-gamma-1).*(beta2-gamma)*V(wnStart:Nn).^(beta2-gamma-2) - ...
.5*gamma*sigma0^2*(gamma-1)*(gamma-2).*V(wnStart:Nn).^(gamma-3)).* ...
.5.*sigma0^2.*V(wnStart:Nn).^(2*gamma)*dt^2/2;
for nn=wnStart:Nn
yMu0dt(nn,ynStart:Mm)= (a*X(nn,ynStart:Mm).^(alpha1-gammaX)+b*X(nn,ynStart:Mm).^(alpha2-gammaX) - ...
.5*gammaX*sigma1^2*V(nn).^(2*gammaV).*X(nn,ynStart:Mm).^(gammaX-1))*dt + ...
(a*(alpha1-gammaX).*X(nn,ynStart:Mm).^(alpha1-gammaX-1)+b*(alpha2-gammaX)*X(nn,ynStart:Mm).^(alpha2-gammaX-1) - ...
.5*gammaX*sigma1^2*V(nn).^(2*gammaV).*(gammaX-1)*X(nn,ynStart:Mm).^(gammaX-2)).* ...
(a*X(nn,ynStart:Mm).^(alpha1)+b*X(nn,ynStart:Mm).^(alpha2))*dt^2/2+ ...
-(.5*gammaX*sigma1^2*(2*gammaV).*V(nn).^(2*gammaV-1).*X(nn,ynStart:Mm).^(gammaX-1)).* ...
(mu1*V(nn).^(beta1)+mu2*V(nn).^(beta2))*dt^2/2 + ...
(a*(alpha1-gammaX).*(alpha1-gammaX-1).*X(nn,ynStart:Mm).^(alpha1-gammaX-2)+b*(alpha2-gammaX-1).*(alpha2-gammaX)*X(nn,ynStart:Mm).^(alpha2-gammaX-2) - ...
.5*gammaX*sigma1^2*(gammaX-1)*(gammaX-2)*V(nn).^(2*gammaV).*X(nn,ynStart:Mm).^(gammaX-3)).* ...
.5*sigma1^2*V(nn).^(2*gammaV).*X(nn,ynStart:Mm).^(2*gammaX)*dt^2/2+ ...
-(.5*gammaX*sigma1^2*(2*gammaV).*(2*gammaV-1).*V(nn).^(2*gammaV-2).*X(nn,ynStart:Mm).^(gammaX-1)).* ...
(sigma0^2*V(nn).^(2*gamma))*dt^2/2;
end

[dwdZ,d2wdZ2] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z);

for nn=wnStart:Nn
yA(ynStart:Mm)=y(nn,ynStart:Mm);
[dydZ0A,d2ydZ02A] = First2Derivatives2ndOrderEqSpacedA(ynStart,Mm,dMm,yA,Z0);
dydZ0(nn,ynStart:Mm)=dydZ0A(ynStart:Mm);
d2ydZ02(nn,ynStart:Mm)=d2ydZ02A(ynStart:Mm);
d2ydZ02(nn,Mm-1:Mm)=0.0;
end

%  for nn=ynStart:Mm
%      yA(wnStart:Nn)=y(wnStart:Nn,nn);
%      [dydZA,d2ydZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,yA,Z);
%      dydZ(wnStart:Nn,nn)=dydZA(wnStart:Nn);
%      d2ydZ2(wnStart:Nn,nn)=d2ydZ2A(wnStart:Nn);
%      %d2ydZ02(nn,Mm-1:Mm)=0.0;
%  end
if(tt>7)
dwb(wnStart:Nn)= wMu0dt(wnStart:Nn)+.5*sigma0^2*dt*Z(wnStart:Nn).*dwdZ(wnStart:Nn).^(-1)+ ...
.5*sigma0^2*dt*(dwdZ(wnStart:Nn).^(-2)).*d2wdZ2(wnStart:Nn);
w(wnStart:Nn)=w(wnStart:Nn)+ dwb(wnStart:Nn);
elseif((tt<=7)&&(tt>1))
dwb(wnStart:Nn)= wMu0dt(wnStart:Nn)+.5*sigma0^2*dt*Z(wnStart:Nn).*dwdZ(wnStart:Nn).^(-1);

w(wnStart:Nn)=w(wnStart:Nn)+ dwb(wnStart:Nn);

elseif((tt==1))
dwb(wnStart:Nn)= wMu0dt(wnStart:Nn)+sigma0*sqrt(dt).*Z(wnStart:Nn);
w(wnStart:Nn)=w(wnStart:Nn)+ dwb(wnStart:Nn);

end

Vnew(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));

dwdZPrev(wnStart:Nn)=dwdZ(wnStart:Nn);
d2wdZ2Prev(wnStart:Nn)=d2wdZ2(wnStart:Nn);
%dwdZPrev(Nn)=dwdZPrev(Nn)*.999;
[dwdZNew,d2wdZ2New] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z);

Ratio(wnStart:Nn,ynStart:Mm)=(dwdZPrev(wnStart:Nn)*dNn+.25*d2wdZ2Prev(wnStart:Nn)*dNn^2)/(dwdZNew(wnStart:Nn)*dNn+.25*d2wdZ2New(wnStart:Nn)*dNn^2);

tt

if(tt>=3)

for nn=wnStart:Nn

yv(nn,ynStart:Mm)=.5*sigma1^2*((V(nn).^(2*gammaV)).*dt).* ...
(Z0(ynStart:Mm).*dydZ0(nn,ynStart:Mm).^(-1).* Ratio(nn,ynStart:Mm).^1+ ...
(dydZ0(nn,ynStart:Mm).^(-2)).*d2ydZ02(nn,ynStart:Mm).* Ratio(nn,ynStart:Mm).^1);% + ...
end
y(wnStart:Nn,ynStart:Mm)=y(wnStart:Nn,ynStart:Mm)+yMu0dt(wnStart:Nn,ynStart:Mm)+ ...
yv(wnStart:Nn,ynStart:Mm);

elseif ((tt<=2)&&(tt>1))
for nn=wnStart:Nn
Ratio(nn,ynStart:Mm)=(dwdZPrev(nn)./dwdZNew(nn));
yv(nn,ynStart:Mm)=.5*V(nn).^(2*gammaV).*dt.* ...
(Z0(ynStart:Mm).*dydZ0(nn,ynStart:Mm).^(-1)).*Ratio(nn,ynStart:Mm);
end

y(wnStart:Nn,ynStart:Mm)=y(wnStart:Nn,ynStart:Mm)+yMu0dt(wnStart:Nn,ynStart:Mm)+yv(wnStart:Nn,ynStart:Mm);
elseif(tt==1)
for mm=1:Mm
yv(wnStart:Nn,mm)=sigma1*V(:).^(gammaV).*sqrt(dt)*Z0(mm);%- ...
end
y(wnStart:Nn,ynStart:Mm)=y(wnStart:Nn,ynStart:Mm)+yMu0dt(wnStart:Nn,ynStart:Mm)+yv(wnStart:Nn,ynStart:Mm);
end

end
Vw(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
DfVw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1

DfVw(nn) = (Vw(nn + 1) - Vw(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end

pVw(1:Nn)=0;
for nn = wnStart:Nn-1

pVw(nn) = (normpdf(Z(nn),0, 1))/abs(DfVw(nn));
end

%Below algorithm caclulates minX and maxX which are the end points
%of new grid constructed to calculate the marginal probability of asset X.
%Any values larger than maxX are not considered towards calculation of
%probability density of X. The current grid has 400 cells equidistributed
%between minX and maxX.
minX=10;
maxX=.001;
for mm=ynStart:Mm
%Xy(mm)=sum(((1-gammaX).*y(wnStart:Nn,mm)).^(1/(1-gammaX)).*ZProb(wnStart:Nn));
Xy_limits(mm)=sum(((1-gammaX).*y(45,mm)).^(1/(1-gammaX)));

end

Xy_limits
if(Xy_limits(4)<minX)
minX=real(Xy_limits(4));
end
if(Xy_limits(33)>maxX)
maxX=real(Xy_limits(31));
end

minX=minX-.025;
maxX=maxX+2.50;
minX=max(.001,minX);
minX
maxX
%minX=.001;
%maxX=6.0;
Mm0=400;
plot(X(1:Mm0),XProb(1:Mm0),'r');
XCDF(1:Mm0)=0.0;
XCDF(1)=XProb(1);
for mm0=2:Mm0
XCDF(mm0)=XCDF(mm0-1)+XProb(mm0);
end
%XCDF=XCDF(:)./sum(XProb(:));
sum(XProb(:))
Zb(1:Mm0)=norminv(XCDF(1:Mm0));%is calculated at boundaries of cells
ZX(2:Mm0)=(Zb(2:Mm0)+Zb(1:Mm0-1))/2.0;%is calculated at center of grid cells.
dXdZX(3:Mm0-1)=(X(4:Mm0)-X(2:Mm0-2))./(ZX(4:Mm0)-ZX(2:Mm0-2))
[dXdZX(3)] = InterpolateOrderN6(4,ZX(3),ZX(4),ZX(5),ZX(6),ZX(7),0,0,dXdZX(4),dXdZX(5),dXdZX(6),dXdZX(7),0,0);
[dXdZX(2)] = InterpolateOrderN6(4,ZX(2),ZX(4),ZX(5),ZX(6),ZX(7),0,0,dXdZX(4),dXdZX(5),dXdZX(6),dXdZX(7),0,0);

%for nn = wnStart:Nn-1
for mm=2:398
pX(mm) = (normpdf(ZX(mm),0, 1))/abs(dXdZX(mm));
end
ZX
pX
pX(isnan(pX)==1)=0.0;
pX(isinf(pX)==1)=0.0;
%pX
XProb
XX=X;
plot(XX(2:380),pX(2:380),'r');
str=input('Look at graph');

for mm=ynStart:Mm
%Xy(mm)=sum(((1-gammaX).*y(wnStart:Nn,mm)).^(1/(1-gammaX)).*ZProb(wnStart:Nn));
Xy(mm)=sum(((1-gammaX).*y(:,mm)).^(1/(1-gammaX)).*ZProb(:));
end
DfXy(ynStart:Mm)=0;
for mm=ynStart+2:Mm-2

DfXy(mm) = (Xy(mm + 1) - Xy(mm - 1))/(Z0(mm + 1) - Z0(mm - 1));
%Change of variable derivative for densities
end

pXy(ynStart:Mm)=0;
for mm = ynStart+2:Mm-2

pXy(mm) = (normpdf(Z0(mm),0, 1))/abs(DfXy(mm));
end

toc

%Xy
ItoHermiteMean=sum(Vw(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1))
TrueMean=thetaV+(V0-thetaV)*exp(-kappaV*dt*Tt)

ItoHermiteMeanX=sum(Xy(ynStart+1:Mm-1).*Z0Prob(ynStart+1:Mm-1))
ItoHermiteMeanX2=sum(X(:).*XProb(:))
TrueMean=thetaX+(X0-thetaX)*exp(-kappaX*dt*Tt)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%rng(29079137, 'twister')
rng(15898837, 'twister')
paths=100000;
V(1:paths)=V0;  %Original process monte carlo.
X=0.0;
X(1:paths)=X0;
Random1(1:paths)=0;
Random2(1:paths)=0;

for ttM=1:TtM
Random1=randn(size(Random1));
Random2=randn(size(Random2));

X(1:paths)=X(1:paths)+ ...
(a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dtM + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) * sqrt(dtM) + ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*sqrt(dtM) + ...
(a*alpha1* X(1:paths).^(alpha1-1)+b*alpha2* X(1:paths).^(alpha2-1)).* ...
(((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dtM^2/2)+ ...
(sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) *(1-1/sqrt(3)).*dtM^1.5+ ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*(1-1/sqrt(3)).*dtM^1.5))+ ...
.5*(a*alpha1*(alpha1-1)* X(1:paths).^(alpha1-2)+b*alpha2*(alpha2-1).* X(1:paths).^(alpha2-2)).* ...
( sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX)) *dtM^2/2 + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random1(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*(Random1(1:paths).^2-1) * dtM/2 + ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths)*dtM/2)+ ...
.5*sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
(sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random1(1:paths).*1/sqrt(3).*dtM^1.5)+ ...
sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random1(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
sigma0*V(1:paths).^gamma.*Random1(1:paths).*Random2(1:paths)*dtM/2)+ ...
.5*sqrt(1-rho^2)* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
(sigma0^2*V(1:paths).^(2*gamma).*Random1(1:paths)*1/sqrt(3)*dtM^1.5)+ ...
sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random1(1:paths)*1/sqrt(3)*dtM^1.5+ ...
rho* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random2(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths) * dtM/2 + ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*(Random2(1:paths).^2-1)*dtM/2)+ ...
.5*rho* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
(sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random2(1:paths).*1/sqrt(3).*dtM^1.5)+ ...
rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
sigma0*V(1:paths).^gamma.*(Random2(1:paths).^2-1)*dtM/2)+ ...
.5*rho* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
sigma0^2.*V(1:paths).^(2*gamma).*Random2(1:paths) * 1/sqrt(3).* dtM^1.5+ ...
rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random2(1:paths)*1/sqrt(3)*dtM^1.5;

V(1:paths)=V(1:paths)+ ...
(mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dtM + ...
sigma0*V(1:paths).^gamma .*Random2(1:paths)*sqrt(dtM) + ...
(mu1.*beta1*V(1:paths).^(beta1-1) + mu2.*beta2.*V(1:paths).^(beta2-1)).* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dtM^2/2 + ...
sigma0*V(1:paths).^gamma .*Random2(1:paths)*(1-1/sqrt(3))*dtM^1.5) + ...
.5*(mu1.*beta1.*(beta1-1).*V(1:paths).^(beta1-2) + mu2.*beta2.*(beta2-1).*V(1:paths).^(beta2-2)).* ...
sigma0^2.*V(1:paths).^(2*gamma).*dtM^2/2 + ...
sigma0*gamma*V(1:paths).^(gamma-1) .* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths).*1/sqrt(3)*dtM^1.5 + ...
sigma0.*V(1:paths).^gamma .*(Random2(1:paths).^2-1)*dtM/2) + ...
.5*sigma0*gamma*(gamma-1).*V(1:paths).^(gamma-2) .* ...
sigma0^2.*V(1:paths).^(2*gamma) .*Random2(1:paths).*1/sqrt(3)*dtM^1.5;

end

SVolMeanAnalytic=thetaV+(V0-thetaV)*exp(-kappaV*dt*Tt)
SVolMean=sum(V(1:paths))/paths
AssetMeanAnalytic=thetaX+exp(-kappaX*T)*(X0-thetaX)
AssetMean=sum(X(1:paths))/paths

MaxCutOff=30;
%NoOfBins=round(600*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too
NoOfBins=1000;
[XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti_NEW(X,paths,NoOfBins,MaxCutOff );
%plot(IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g',Xy(ynStart+2:Mm-2),pXy(ynStart+2:Mm-2),'r');
plot(IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g',XX(2:398),pX(2:398),'r');
%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(sprintf('X0 = %.4f,thetaX=%.3f,kappaX=%.2f,gammaX=%.3f,gammaV=%.3f,rho=%.2f,sigma1=%.2f,T=%.2f,dt=%.5f,AM=%.4f,M=%.4f,TM=%.4f', X0,thetaX,kappaX,gammaX,gammaV,rho,sigma1,T,dt,ItoHermiteMeanX2,AssetMean,AssetMeanAnalytic));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

%legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')

str=input('Green line is density of Asset SDE.');

MaxCutOff=30;
%NoOfBins=round(600*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too
NoOfBins=2000;
[VDensity,IndexOutV,IndexMaxV] = MakeDensityFromSimulation_Infiniti_NEW(V,paths,NoOfBins,MaxCutOff );
plot(IndexOutV(1:IndexMaxV),VDensity(1:IndexMaxV),'g',Vw(wnStart+1:Nn-1),pVw(wnStart+1:Nn-1),'r');
%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');
%plot(v(wnStart+1:Nn-1),pv(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');

title(sprintf('V0 = %.4f,thetaV=%.3f,kappaV=%.2f,gamma=%.3f,sigma0=%.2f,T=%.2f,dt=%.5f,SVolMean=%.4f,SVolMeanAnalytic=%.4f', V0,thetaV,kappaV,gamma,sigma0,T,dt,SVolMean,SVolMeanAnalytic));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

%legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')

str=input('red line is density of Asset SDE.');

end

.
You will need three sub-functions. The first one.
.
function [X,XProb] = AddXAcrossSVolGrids(y,wnStart,Nn,ynStart,Mm,dMm,Mm0,ZProb,gammaX,minX,maxX)

%Mm0=75;
%maxX=0;
%minX=10;

Xy(wnStart:Nn,ynStart:Mm)=((1-gammaX).*y(wnStart:Nn,ynStart:Mm)).^(1/(1-gammaX));
Xy(real(Xy)<0)=0;
%for nn=1:Nn
%    for mm=1:Mm
%        maxX=max(maxX,Xy(nn,mm));
%        minX=min(minX,Xy(nn,mm));
%    end
%end

maxX
minX
str=input('Look at numbers');

X(1)=minX;
X(Mm0)=maxX;
dX=(X(Mm0)-X(1))/(Mm0-1);
X(1:Mm0)=X(1)+(0:Mm0-1)*dX;
Xa(1:Mm0)=X(1:Mm0)-dX/2;
Xb(1:Mm0)=X(1:Mm0)+dX/2;
Xa(1)=X(1)-dX/20;
if(Xa(1)<0)
Xa(1)=.0000001;
end
XProb(1:Mm0)=0.0;

MmA=200;
dMmA=.05;
NnMid=4.0;
Z0A(1:MmA)=(((1:MmA)-20.5)*dMmA-NnMid);
Z0Prob(1:MmA)=0.0;
Z0Prob(1)=normcdf(.5*Z0A(1)+.5*Z0A(2),0,1)-normcdf(.5*Z0A(1)+.5*Z0A(2)-dMmA,0,1);
Z0Prob(MmA)=normcdf(.5*Z0A(MmA)+.5*Z0A(MmA-1)+dMmA,0,1)-normcdf(.5*Z0A(MmA)+.5*Z0A(MmA-1),0,1);
Z0Prob(2:MmA-1)=normcdf(.5*Z0A(2:MmA-1)+.5*Z0A(3:MmA),0,1)-normcdf(.5*Z0A(2:MmA-1)+.5*Z0A(1:MmA-2),0,1);

for nn=wnStart+2:Nn-2
Xy0(ynStart:Mm)=Xy(nn,ynStart:Mm);
XyI=InterpolateFPEGrid40To200(Xy0,ynStart,Mm);
[Xya0,Xyb0] = CalculateGridStartsAndEndsBoundedEnds(XyI,Z0A,ynStart,MmA,dMmA);
Xya(nn,ynStart:MmA)=real(Xya0(ynStart:MmA));
Xyb(nn,ynStart:MmA)=real(Xyb0(ynStart:MmA));
Xya(nn,Xya<0)=0.0;
Xyb(nn,Xyb<0)=0.0;

end

for nn=wnStart+2:Nn-2

%DfXy(nn,ynStart:Mm)=0;
%for mm=ynStart+1:Mm-1

%    dXydZ0(nn, mm) = (Xy(nn,mm + 1) - Xy(nn,mm - 1))/(Z0(mm + 1) - Z0(mm - 1));
%Change of variable derivative for densities
%end

mm=1;
mm0=1;
while((mm<MmA)&&(mm0<Mm0))
%while((mm<Mm)||(mm0<Mm0)||((mm==Mm)&&(Xyb(mm)<Xb(Mm0))))
%if((mm==Mm)&&(Xyb(mm)<Xb(Mm0))

if((Xya(nn,mm)>=Xa(mm0))&&(Xyb(nn,mm)<Xb(mm0)))
XProb(mm0)=XProb(mm0)+Z0Prob(mm).*ZProb(nn);
if(mm<MmA)
mm=mm+1;
end
elseif((Xya(nn,mm)>=Xa(mm0))&&(Xyb(nn,mm)==Xb(mm0)))
XProb(mm0)=XProb(mm0)+Z0Prob(mm).*ZProb(nn);
if(mm<MmA)
mm=mm+1;
end
if(mm0<Mm0)
mm0=mm0+1;
end
elseif((Xya(nn,mm)>=Xa(mm0))&&(Xyb(nn,mm)>Xb(mm0)) &&(Xya(nn,mm)<=Xb(mm0)))
XProb(mm0)=XProb(mm0)+Z0Prob(mm).*ZProb(nn).*(Xb(mm0)-Xya(nn,mm))./(Xyb(nn,mm)-Xya(nn,mm));
if(mm0<Mm0)
mm0=mm0+1;
end
elseif((Xya(nn,mm)<Xa(mm0))&&(Xyb(nn,mm)<Xb(mm0))&&(Xyb(nn,mm)>Xa(mm0)))
XProb(mm0)=XProb(mm0)+Z0Prob(mm).*ZProb(nn).*(Xyb(nn,mm)-Xa(mm0))./(Xyb(nn,mm)-Xya(nn,mm));
if(mm<MmA)
mm=mm+1;
end
elseif((Xya(nn,mm)<Xa(mm0))&&(Xyb(nn,mm)==Xb(mm0))&&(Xyb(nn,mm)>Xa(mm0)))
XProb(mm0)=XProb(mm0)+Z0Prob(mm).*ZProb(nn).*(Xyb(nn,mm)-Xa(mm0))./(Xyb(nn,mm)-Xya(nn,mm));
if(mm<MmA)
mm=mm+1;
end
if(mm0<Mm0)
mm0=mm0+1
end
elseif((Xya(nn,mm)<Xa(mm0))&&(Xyb(nn,mm)>Xb(mm0)))
XProb(mm0)=XProb(mm0)+Z0Prob(mm).*ZProb(nn).*(Xb(mm0)-Xa(mm0))./(Xyb(nn,mm)-Xya(nn,mm));
if(mm0<Mm0)
mm0=mm0+1;
end
elseif((Xya(nn,mm)<Xa(mm0))&&(Xyb(nn,mm)<=Xa(mm0)))
if(mm<MmA)
mm=mm+1;
end
elseif((Xya(nn,mm)>=Xb(mm0))&&(Xyb(nn,mm)>Xb(mm0)))
if(mm0<Mm0)
mm0=mm0+1;
end
end
end

end
end


.
.
2nd sub-function.
.
function [wI] = InterpolateFPEGrid40To200(w,wnStart,Nn)
%UNTITLED2 Summary of this function goes here
%   Detailed explanation goes here

dNn=.25;
dMm=.05;
Nn=40;
Mm=200;
NnMid=4.0;
Z(1:Nn)=(((1:Nn)-4.5)*dNn-NnMid);
Z2(1:Mm)=(((1:Mm)-20.5)*dMm-NnMid);
wI(1:Mm)=0.0;

for nn=wnStart+3:Nn-4
nn;
Z0=(((nn)-4.5)*dNn-NnMid);
%Z01=Z0+.025;
%Z02=Z0+.075;
%Z03=Z0+.125;
%Z04=Z0+.175;

Z01=Z0+.05;
Z02=Z0+.1;
Z03=Z0+.15;
Z04=Z0+.2;

mm0=round((Z0+NnMid)/dMm+20.5);
mm1=round((Z01+NnMid)/dMm+20.5);
mm2=round((Z02+NnMid)/dMm+20.5);
mm3=round((Z03+NnMid)/dMm+20.5);
mm4=round((Z04+NnMid)/dMm+20.5);

nn1=nn-3;
nn2=nn-2;
nn3=nn-1;
nn4=nn;
nn5=nn+1;
nn6=nn+2;
nn7=nn+3;
nn8=nn+4;

wI(mm0)=w(nn);
wI(mm1) =InterpolateOrderN8(8,Z2(mm1),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),Z(nn7),Z(nn8),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6),w(nn7),w(nn8));
wI(mm2) =InterpolateOrderN8(8,Z2(mm2),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),Z(nn7),Z(nn8),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6),w(nn7),w(nn8));
wI(mm3) =InterpolateOrderN8(8,Z2(mm3),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),Z(nn7),Z(nn8),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6),w(nn7),w(nn8));
wI(mm4) =InterpolateOrderN8(8,Z2(mm4),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),Z(nn7),Z(nn8),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6),w(nn7),w(nn8));

end

nn=wnStart;
nn;
nn1=nn;
nn2=nn+1;
nn3=nn+2;
nn4=nn+3;
nn5=nn+4;
nn6=nn+5;

Z0=(((nn)-4.5)*dNn-NnMid);

Z01=Z0-.1;
Z02=Z0-.05;
Z00=Z0;
Z03=Z0+.05;
Z04=Z0+.1;
Z05=Z0+.15;
Z06=Z0+.2;

mm0=round((Z0+NnMid)/dMm+20.5);
mm1=round((Z01+NnMid)/dMm+20.5);
mm2=round((Z02+NnMid)/dMm+20.5);
mm3=round((Z03+NnMid)/dMm+20.5);
mm4=round((Z04+NnMid)/dMm+20.5);
mm5=round((Z05+NnMid)/dMm+20.5);
mm6=round((Z06+NnMid)/dMm+20.5);

wI(mm0) =w(nn);
wI(mm1) =InterpolateOrderN6(6,Z2(mm1),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm2) =InterpolateOrderN6(6,Z2(mm2),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm3) =InterpolateOrderN6(6,Z2(mm3),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm4) =InterpolateOrderN6(6,Z2(mm4),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm5) =InterpolateOrderN6(6,Z2(mm5),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm6) =InterpolateOrderN6(6,Z2(mm6),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));

nn=wnStart+1;
nn;
nn1=nn-1;
nn2=nn;
nn3=nn+1;
nn4=nn+2;
nn5=nn+3;
nn6=nn+4;

Z0=(((nn)-4.5)*dNn-NnMid);
Z01=Z0+.05;
Z02=Z0+.1;
Z03=Z0+.15;
Z04=Z0+.2;

mm0=round((Z0+NnMid)/dMm+20.5);
mm1=round((Z01+NnMid)/dMm+20.5);
mm2=round((Z02+NnMid)/dMm+20.5);
mm3=round((Z03+NnMid)/dMm+20.5);
mm4=round((Z04+NnMid)/dMm+20.5);

wI(mm0)=w(nn);
wI(mm1) =InterpolateOrderN6(6,Z2(mm1),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm2) =InterpolateOrderN6(6,Z2(mm2),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm3) =InterpolateOrderN6(6,Z2(mm3),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm4) =InterpolateOrderN6(6,Z2(mm4),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));

nn=wnStart+2;
nn;
nn1=nn-2;
nn2=nn-1;
nn3=nn;
nn4=nn+1;
nn5=nn+2;
nn6=nn+3;
nn7=nn+4;

Z0=(((nn)-4.5)*dNn-NnMid);
Z01=Z0+.05;
Z02=Z0+.1;
Z03=Z0+.15;
Z04=Z0+.2;

mm0=round((Z0+NnMid)/dMm+20.5);
mm1=round((Z01+NnMid)/dMm+20.5);
mm2=round((Z02+NnMid)/dMm+20.5);
mm3=round((Z03+NnMid)/dMm+20.5);
mm4=round((Z04+NnMid)/dMm+20.5);

wI(mm0)=w(nn);
wI(mm1) =InterpolateOrderN8(7,Z2(mm1),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),Z(nn7),0,w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6),w(nn7),0);
wI(mm2) =InterpolateOrderN8(7,Z2(mm2),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),Z(nn7),0,w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6),w(nn7),0);
wI(mm3) =InterpolateOrderN8(7,Z2(mm3),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),Z(nn7),0,w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6),w(nn7),0);
wI(mm4) =InterpolateOrderN8(7,Z2(mm4),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),Z(nn7),0,w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6),w(nn7),0);

nn=Nn-3;
nn;
nn1=nn-3;
nn2=nn-2;
nn3=nn-1;
nn4=nn;
nn5=nn+1;
nn6=nn+2;
nn7=nn+3;

Z0=(((nn)-4.5)*dNn-NnMid);
Z01=Z0+.05;
Z02=Z0+.1;
Z03=Z0+.15;
Z04=Z0+.2;

mm0=round((Z0+NnMid)/dMm+20.5);
mm1=round((Z01+NnMid)/dMm+20.5);
mm2=round((Z02+NnMid)/dMm+20.5);
mm3=round((Z03+NnMid)/dMm+20.5);
mm4=round((Z04+NnMid)/dMm+20.5);

wI(mm0)=w(nn);
wI(mm1) =InterpolateOrderN8(7,Z2(mm1),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),Z(nn7),0,w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6),w(nn7),0);
wI(mm2) =InterpolateOrderN8(7,Z2(mm2),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),Z(nn7),0,w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6),w(nn7),0);
wI(mm3) =InterpolateOrderN8(7,Z2(mm3),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),Z(nn7),0,w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6),w(nn7),0);
wI(mm4) =InterpolateOrderN8(7,Z2(mm4),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),Z(nn7),0,w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6),w(nn7),0);

nn=Nn-2;
nn;
nn1=nn-3;
nn2=nn-2;
nn3=nn-1;
nn4=nn;
nn5=nn+1;
nn6=nn+2;

Z0=(((nn)-4.5)*dNn-NnMid);
Z01=Z0+.05;
Z02=Z0+.1;
Z03=Z0+.15;
Z04=Z0+.2;

mm0=round((Z0+NnMid)/dMm+20.5);
mm1=round((Z01+NnMid)/dMm+20.5);
mm2=round((Z02+NnMid)/dMm+20.5);
mm3=round((Z03+NnMid)/dMm+20.5);
mm4=round((Z04+NnMid)/dMm+20.5);

wI(mm0)=w(nn);
wI(mm1) =InterpolateOrderN6(6,Z2(mm1),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm2) =InterpolateOrderN6(6,Z2(mm2),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm3) =InterpolateOrderN6(6,Z2(mm3),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm4) =InterpolateOrderN6(6,Z2(mm4),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));

nn=Nn-1;
nn;
nn1=nn-4;
nn2=nn-3;
nn3=nn-2;
nn4=nn-1;
nn5=nn;
nn6=nn+1;

Z0=(((nn)-5.5)*dNn-NnMid);
Z01=Z0+.05;
Z02=Z0+.1;
Z03=Z0+.15;
Z04=Z0+.2;
Z05=Z0+.25;
Z06=Z0+.3;
Z07=Z0+.35;

mm0=round((Z0+NnMid)/dMm+20.5);
mm1=round((Z01+NnMid)/dMm+20.5);
mm2=round((Z02+NnMid)/dMm+20.5);
mm3=round((Z03+NnMid)/dMm+20.5);
mm4=round((Z04+NnMid)/dMm+20.5);
mm5=round((Z05+NnMid)/dMm+20.5);
mm6=round((Z06+NnMid)/dMm+20.5);
mm7=round((Z07+NnMid)/dMm+20.5);

wI(mm0)=w(nn);
wI(mm5)=w(nn+1);
wI(mm1) =InterpolateOrderN6(6,Z2(mm1),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm2) =InterpolateOrderN6(6,Z2(mm2),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm3) =InterpolateOrderN6(6,Z2(mm3),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm4) =InterpolateOrderN6(6,Z2(mm4),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
%wI(mm5) =InterpolateOrderN6(6,Z2(mm5),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm6) =InterpolateOrderN6(6,Z2(mm6),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
wI(mm7) =InterpolateOrderN6(6,Z2(mm7),Z(nn1),Z(nn2),Z(nn3),Z(nn4),Z(nn5),Z(nn6),w(nn1),w(nn2),w(nn3),w(nn4),w(nn5),w(nn6));
end

.
Thirs sub-function.
.
function [GridStarts,GridEnds] = CalculateGridStartsAndEndsBoundedEnds(w,Z,wnStart,Nn,dNn)
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here

for nn=wnStart+2:Nn-3
x0(nn-wnStart-2+1)=Z(nn)+.5*dNn;
x1(nn-wnStart-2+1)=Z(nn-2);
x2(nn-wnStart-2+1)=Z(nn-1);
x3(nn-wnStart-2+1)=Z(nn);
x4(nn-wnStart-2+1)=Z(nn+1);
x5(nn-wnStart-2+1)=Z(nn+2);
x6(nn-wnStart-2+1)=Z(nn+3);

y1(nn-wnStart-2+1)=w(nn-2);
y2(nn-wnStart-2+1)=w(nn-1);
y3(nn-wnStart-2+1)=w(nn);
y4(nn-wnStart-2+1)=w(nn+1);
y5(nn-wnStart-2+1)=w(nn+2);
y6(nn-wnStart-2+1)=w(nn+3);
end

N=6;
[y0] = InterpolateOrderN6(N,x0,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
GridEnds(wnStart+2:Nn-3)=y0(1:Nn-3-wnStart-2+1);

x0=Z(wnStart+1)+.5*dNn;
x1=Z(wnStart);
x2=Z(wnStart+1);
x3=Z(wnStart+2);
x4=Z(wnStart+3);
x5=Z(wnStart+4);
x6=0;

y1=w(wnStart);
y2=w(wnStart+1);
y3=w(wnStart+2);
y4=w(wnStart+3);
y5=w(wnStart+4);
y6=0;
N=5;
[y0] = InterpolateOrderN6(N,x0,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);

GridEnds(wnStart+1)=y0;

x01=Z(wnStart)+.5*dNn;
x02=Z(wnStart)-.5*dNn;
x1=Z(wnStart);
x2=Z(wnStart+1);
x3=Z(wnStart+2);
x4=Z(wnStart+3);

y1=w(wnStart);
y2=w(wnStart+1);
y3=w(wnStart+2);
y4=w(wnStart+3);

N=4;
[y01] = InterpolateOrderN6(N,x01,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
[y02] = InterpolateOrderN6(N,x02,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
GridEnds(wnStart)=y01;
GridStarts(wnStart)=y02;

x0=Z(Nn-2)+.5*dNn;
x1=Z(Nn-4);
x2=Z(Nn-3);
x3=Z(Nn-2);
x4=Z(Nn-1);
x5=Z(Nn);
x6=0;

y1=w(Nn-4);
y2=w(Nn-3);
y3=w(Nn-2);
y4=w(Nn-1);
y5=w(Nn);
y6=0;
N=5;
[y0] = InterpolateOrderN6(N,x0,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);

GridEnds(Nn-2)=y0;

x01=Z(Nn-1)+.5*dNn;
x02=Z(Nn)+.5*dNn;
x1=Z(Nn-3);
x2=Z(Nn-2);
x3=Z(Nn-1);
x4=Z(Nn);
x5=0;
x6=0;

y1=w(Nn-3);
y2=w(Nn-2);
y3=w(Nn-1);
y4=w(Nn);
y5=0;
y6=0;
N=4;
[y01] = InterpolateOrderN6(N,x01,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
[y02] = InterpolateOrderN6(N,x02,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);

GridEnds(Nn-1)=y01;
GridEnds(Nn)=y02;

%GridEnds(Nn)=Inf;
%GridStarts(wnStart)=-Inf;
GridStarts(2:Nn)=GridEnds(1:Nn-1);

end

.
.
Following is the output of the function.
ItoHermiteMean =
0.331502245770817
TrueMean =
0.331201169941968
ItoHermiteMeanX =
1.034410761585111
ItoHermiteMeanX2 =
1.000326543369098
TrueMean =
1
SVolMeanAnalytic =
0.331201169941968
SVolMean =
0.331317893621948
AssetMeanAnalytic =
1
AssetMean =
0.999501853160871
IndexMax =
1001
Green line is density of Asset SDE.
IndexMax =
2001
red line is density of Asset SDE.

Here are the two graphs for SV SDEs you will see.

.

Amin
Topic Author
Posts: 2894
Joined: July 14th, 2002, 3:00 am

### Re: Breakthrough in the theory of stochastic differential equations and their simulation

I want to write some notes and comments for friends who might like to play with the program I posted yesterday.
1. The grid for asset within each stochastic volatility subdivision has 40 grid points. Mm=40. There are 50 SV grid points so there are a total of 50 X 40 grid points. I was earlier using a 50 X 50 grid but my solution would break down and become unstable after 10-20 time steps. For several days I could not know the reason. While attempting to get the program to run, I increased the grid points in asset from 50 to 100 and the grid became very unstable and could not do more than 2-3 points. This prompted me to think if I could decrease grid points in reverse direction to increasing points, the grid might become stable and I tried a grid with size 40 and only then grid became reasonably stable in time with most SDEs. With precisely the same algorithm on grid size 50, I could not go beyond 20-25 small steps in time. There is a lot to be known about stability of algorithm and I am sure numerical PDE experts could discover far better algorithms that can take an arbitrary grid and quite large step sizes in time.

2. In the main simulation of asset price part of the algorithm, I will copy the equation from the program here as

yv(nn,ynStart:Mm)=.5*sigma1^2*((V(nn).^(2*gammaV)).*dt).* ...
(Z0(ynStart:Mm).*dydZ0(nn,ynStart:Mm).^(-1).* Ratio(nn,ynStart:Mm).^1+ ...
(dydZ0(nn,ynStart:Mm).^(-2)).*d2ydZ02(nn,ynStart:Mm).* Ratio(nn,ynStart:Mm).^1);

The whole line is derived from the solution of  similar 1st order PDEs but it is multiplied with  Ratio(nn,ynStart:Mm) where
Ratio(nn,ynStart:Mm)=(dwdZPrev(nn)./dwdZNew(nn))
When I wrote the algorithm without ratio, there was large slack and then I had the intuitive idea that both sides of the square are changing so while calculating the change in length of conditional side, we may have to account for how much it has changed in the unconditional side already which is given by Ratio(nn,ynStart:Mm). I did not write proper equations but just now I am thinking that we can probably write the equation for probability mass within each grid cell at time t+dt after application of Fokker-planck equation and equate it with probability of grid cell at previous time t.

3. In our 2D solution of fokker-planck, the asset value is based on Z-grid or CDF lines that are different for each value of SV grid. There are no absolute values on asset grid. In a situation like this, if we have to add the probability of asset across different SV grids, we have to first convert the asset values from Z-grid(that is different for each value of SV on SV grid) to some constant grid and then add probability across the constant grid. Sometimes due to weakness in the simple algorithm or possibly due to volatility explosions, there are cells on extreme side of volatility that blow or become unstable but with extremely small probability. So I came up with a quick ad-hoc way to construct this constant grid so that whole of the interesting domain is covered. Here is the algorithm I used

minX=10;
maxX=.001;
for mm=ynStart:Mm
%Xy(mm)=sum(((1-gammaX).*y(wnStart:Nn,mm)).^(1/(1-gammaX)).*ZProb(wnStart:Nn));
Xy_limits(mm)=sum(((1-gammaX).*y(45,mm)).^(1/(1-gammaX)));

end

Xy_limits

if(Xy_limits(4)<minX)
minX=real(Xy_limits(4));
end
if(Xy_limits(33)>maxX)
maxX=real(Xy_limits(31));             %please convert both to 33 or 31.
end

minX=minX-.025;
maxX=maxX+2.50;
minX=max(.001,minX);

In the above copied code, minX is the starting value of the constant asset grid and maxX is the highest value(largest grid point) of the constant asset grid. For this I chose  "Xy_limits(mm)=sum(((1-gammaX).*y(45,mm)).^(1/(1-gammaX)));", please note that y(45,mm) means that we have chosen asset values associated with +4 SD on the volatility grid.(45 out of total 50 subdivisions corresponds to +4 on Z-grid) and then I chose
minX=real(Xy_limits(4));  meaning -4 SD on asset grid(and +4 SD on vol grid)for minimum value and similarly
maximum value is chosen equal to 31 asset grid point on asset grid (and +4 SD on vol grid). I added 2.50 on top of this maximum. This is all very ad-hoc but works reasonably for an initial program.

4.

The above function returns a constant X-grid and probability mass associated with each grid point. It takes size of grid Mm0=400 as an input. If this function goes into an infinite loop,  then reason is that there are too many NaNs in the y (which is 40X 50 value of asset on 2D grid) input to the algorithm or it has gone into zero and these cases are not covered yet. The algorithm is very basic and becomes unstable at too high volatility so if you come across an infinite loop here, you would have to abort by pressing  Ctrl+C .

Amin
Topic Author
Posts: 2894
Joined: July 14th, 2002, 3:00 am

### Re: Breakthrough in the theory of stochastic differential equations and their simulation

Friends, I expect that 2D SV PDE with correlations is solvable based on conservation of mass just like we did for 1D PDE from first principles. I have had some success and I hope that I will be able to solve the problem in a few days and would post the solution here.

Amin
Topic Author
Posts: 2894
Joined: July 14th, 2002, 3:00 am

### Re: Breakthrough in the theory of stochastic differential equations and their simulation

I have been able to solve the 2D SV PDE without correlations from principle of conservation of mass just like we did for one dimensional FPE PDE. In a day or two, I will be posting the program and also would write equations here on the forum that I derived from principle of conservation of mass to solve the PDE.

Amin
Topic Author
Posts: 2894
Joined: July 14th, 2002, 3:00 am

### Re: Breakthrough in the theory of stochastic differential equations and their simulation

I am going to share the analytics for solution of 2D stochastic volatility equations. Before I give the solution, I want to mention that in typical SV system of SDEs, there are two SDEs out of which one (usually Stochastic volatility) is independent of the other SDE and its one dimensional PDE can be solved independently. If we had a system of two SDEs which are both dependent on each other, in that case the solution though still possible would be more difficult.

I have taken both asset and SV Sdes in bessel coordinates. SV SDE can be solved independently and its PDE is given as

$\frac{\partial P(w,t)}{\partial t} \, = \, -\frac{\partial [\mu_w \, P(w,t)]}{\partial w} \, +\, .5 \, \frac{\partial^2 [{\sigma_w}^2 \, P(w,t)]}{\partial w^2} \,$

When we solve this equation with conservation of mass as I showed in a previous post, we get

$-\mu_w \, \frac{\partial Z_1}{\partial w} - .5 {\sigma_w}^2 \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2- P(Z_1)\, {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \, \big] \,$
$+ \frac{\partial w_b}{\partial t} \frac{\partial Z_1}{\partial w} =0$

In the above equation b in $w_b$ denotes that it is the boundary point that is moving so as to keep the probability conserved up till that point. We will need the above equation because when we move to two dimensional stochastic volatility PDE, the whole of above equation would be integrated in another dimension but since it equates to zero, it will also remain zero even after an extra integration and we will eliminate it from 2D equation.

our two dimensional SV and asset partial differential equation in bessel coordinates that we want to solve is
$\frac{\partial P(w,y,t)}{\partial t} \, = \, -\frac{\partial \big[\mu_y \, P(w,y,t) \big]}{\partial y} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(w,y,t) \big]}{\partial y^2} \, -\frac{\partial \big[\mu_w \, P(w,y,t) \big]}{\partial w} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_w}^2 \, P(w,y,t) \big]}{\partial w^2} \,+ \,\rho \, \frac{\partial^2 \big[ \sigma_y \, \sigma_w \, {V(w)}^{ \gamma_V} \, P(w,y,t) \big]}{\partial w \, \partial y} \,$

though correlation term is conceptually simple and when we do a two dimensional integration to preserve probability mass, this term becomes even more easy but correlation has some other issues when taking derivatives with respect of density of y with respect to Z_1 in diffusion terms. When there is no correlation our jacobian for two dimensional density when taken with respect to two independent standard nromals becomes

$\frac{\partial (Z_1,Z_2)}{\partial (w,y)}=\frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}$

since $Z_1$ and $Z_2$ are independent, we can write

$P(w,y)=P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}$

Now we set ourselves in the probability conservation set up, and we have to write a two dimensional integral from infinity(or rather zero ) to moving boundary (boundary might be a misnomer since this is not grid boundary but rather boundary of mass conserved till a certain grid point and this so called boundary is moving or alternatively the grid point is moving) points $w_b$ and $y_b$. Please note that on a two dimensional grid there are a series of these points $w_b$'s and $y_b$'s. The grid in inependent SV dimension is given by index n and the grid on dependent asset dimension is given by (n,m) and respective moving boundaries in conservation set up are indicated as  $w_b(n)$(that conserves mass till nth point on SV grid) and $y_b(n,m)$ (that conserves mass till (n,m) points on the grid) . Please note that  $y_b(n,m)$ is not a straight line and is curved in both indices.

We apply the conservation of probability mass till arbitrary nth and mth points on the grid as

$\frac{\partial \Big[\int_0^{w_b(n)} \, \int_0^{y_b(n,m)} \, P(w,y,t) \,\, dy \, dw \Big]}{\partial t} \, =0$
$=\int_0^{w_b(n)} \, \int_0^{y_b(n,m)} \, \frac{P(w,y,t)}{\partial t} \, \, dy \, dw \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(w_b(n),y,t) dy$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(w,y_b,t) dw =0$

Substituting the 2D uncorrelated PDE in the first RHS integral, we get the expression

$=\int_0^{w_b(n)} \, \int_0^{y_b(n,m)} \, \Big[\, -\frac{\partial \big[\mu_y \, P(w,y,t) \big]}{\partial y} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(w,y,t) \big]}{\partial y^2} \, -\frac{\partial \big[\mu_w \, P(w,y,t) \big]}{\partial w} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_w}^2 \, P(w,y,t) \big]}{\partial w^2} \, \Big] \, \, dy \, dw \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(w_b(n),y,t) dy$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(w,y_b,t) dw =0$

We do the integrations in the first term to get
$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, P(w,y,t) \, +\, .5 \, \frac{\partial \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(w,y,t) \big]}{\partial y} \, \Big] dw$
$\int_0^{y_b(n,m)} \Big[ -\mu_w \, P(w,y,t) \, +\, .5 \, \frac{\partial \big[{\sigma_w}^2 \, P(w,y,t) \big]}{\partial w} \, \Big] \, \, dy \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(w_b(n),y,t) dy$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(w,y_b,t) dw =0$

Now we convert our coordinates to $Z_1$ and $Z_2$ instead of w and y as

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, \frac{\partial \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial y} \, \Big] dw$
$\int_0^{y_b(n,m)} \Big[ -\mu_w \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, \frac{\partial \big[{\sigma_w}^2 \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial w} \, \Big] \, \, dy \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(Z_1(n)) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} dy$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1(n)}{\partial w} \, \frac{\partial Z_2(m)}{\partial y} dw =0$

which can be further simplified as

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$\int_0^{y_b(n,m)} \Big[ -\mu_w \, \frac{\partial Z_1}{\partial w} \, +\, .5 \, {\sigma_w}^2 \, \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2 - {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \big] \Big] \, \, dZ_2 \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, \frac{\partial Z_1}{\partial w} \, dZ_2$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0$

But we notice that second and third term together are merely the solution of 1D independent PDE integrated in an extra dimension and since the original solution of both terms for independent 1D PDE was zero, the integrated term should also be zero. (Even if it were not zero, we could have found this solution from solving the 1D independent PDE), therefore
$\int_0^{y_b(n,m)} \Big[ -\mu_w \, \frac{\partial Z_1}{\partial w} \, +\, .5 \, {\sigma_w}^2 \, \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2 - {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \big] \Big] \, \, dZ_2 \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, \frac{\partial Z_1}{\partial w} \, dZ_2=0$

Now we are left with the terms to solve as

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0$

For a given value of m, we can recursively solve the above PDE for each value of n. To make matters simple, I write the above equation as

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$+ \int_0^{w_b(n-1)} \frac{\partial y_b(n,m)}{\partial t} \, \frac{\partial Z_2}{\partial y} dZ_1$
$+ \int_{w_b(n-1)}^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, \frac{\partial Z_2}{\partial y} dZ_1 =0$

where in the above equation, I have expanded the last terms as sum of all previous terms and one current term to solve and now the equation can be solved as

$\frac{\partial y_b(n,m)}{\partial t} =\Bigg[ \int_0^{w_b(n)} \, \Big[\, \mu_y \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, \big[ Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 + {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$- \int_0^{w_b(n-1)} \frac{\partial y_b(n,m)}{\partial t} \, \frac{\partial Z_2}{\partial y} dZ_1 \Bigg] / \Big[ \int_{w_b(n-1)}^{w_b(n)} \, \frac{\partial Z_2}{\partial y} dZ_1 \Big]$

please note that $\frac{\partial y_b(n,m)}{\partial t} \Delta t$ is the change in coordinates y(n,m) so as to conserve the mass in a small time step according to given PDE dynamics.
The last equation solves the 2D PDE in the asset dimension. Adding correlation term should not be difficult but then we have to make a few changes as correlation will appear in jacobian. Also derivatives of asset with respect to $Z_1$ would have to be included something that was making my equation unstable when I tried it last time but I will try it again to see if I could make it work.
I will be posting the program for uncorrelated PDE based on these analytics tomorrow.

Amin
Topic Author
Posts: 2894
Joined: July 14th, 2002, 3:00 am

### Re: Breakthrough in the theory of stochastic differential equations and their simulation

I am going to share the analytics for solution of 2D stochastic volatility equations. Before I give the solution, I want to mention that in typical SV system of SDEs, there are two SDEs out of which one (usually Stochastic volatility) is independent of the other SDE and its one dimensional PDE can be solved independently. If we had a system of two SDEs which are both dependent on each other, in that case the solution though still possible would be more difficult.

I have taken both asset and SV Sdes in bessel coordinates. SV SDE can be solved independently and its PDE is given as

$\frac{\partial P(w,t)}{\partial t} \, = \, -\frac{\partial [\mu_w \, P(w,t)]}{\partial w} \, +\, .5 \, \frac{\partial^2 [{\sigma_w}^2 \, P(w,t)]}{\partial w^2} \,$

When we solve this equation with conservation of mass as I showed in a previous post, we get

$-\mu_w \, \frac{\partial Z_1}{\partial w} - .5 {\sigma_w}^2 \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2- P(Z_1)\, {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \, \big] \,$
$+ \frac{\partial w_b}{\partial t} \frac{\partial Z_1}{\partial w} =0$

In the above equation b in $w_b$ denotes that it is the boundary point that is moving so as to keep the probability conserved up till that point. We will need the above equation because when we move to two dimensional stochastic volatility PDE, the whole of above equation would be integrated in another dimension but since it equates to zero, it will also remain zero even after an extra integration and we will eliminate it from 2D equation.

our two dimensional SV and asset partial differential equation in bessel coordinates that we want to solve is
$\frac{\partial P(w,y,t)}{\partial t} \, = \, -\frac{\partial \big[\mu_y \, P(w,y,t) \big]}{\partial y} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(w,y,t) \big]}{\partial y^2} \, -\frac{\partial \big[\mu_w \, P(w,y,t) \big]}{\partial w} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_w}^2 \, P(w,y,t) \big]}{\partial w^2} \,+ \,\rho \, \frac{\partial^2 \big[ \sigma_y \, \sigma_w \, {V(w)}^{ \gamma_V} \, P(w,y,t) \big]}{\partial w \, \partial y} \,$

though correlation term is conceptually simple and when we do a two dimensional integration to preserve probability mass, this term becomes even more easy but correlation has some other issues when taking derivatives with respect of density of y with respect to Z_1 in diffusion terms. When there is no correlation our jacobian for two dimensional density when taken with respect to two independent standard nromals becomes

$\frac{\partial (Z_1,Z_2)}{\partial (w,y)}=\frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}$

since $Z_1$ and $Z_2$ are independent, we can write

$P(w,y)=P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}$

Now we set ourselves in the probability conservation set up, and we have to write a two dimensional integral from infinity(or rather zero ) to moving boundary (boundary might be a misnomer since this is not grid boundary but rather boundary of mass conserved till a certain grid point and this so called boundary is moving or alternatively the grid point is moving) points $w_b$ and $y_b$. Please note that on a two dimensional grid there are a series of these points $w_b$'s and $y_b$'s. The grid in inependent SV dimension is given by index n and the grid on dependent asset dimension is given by (n,m) and respective moving boundaries in conservation set up are indicated as  $w_b(n)$(that conserves mass till nth point on SV grid) and $y_b(n,m)$ (that conserves mass till (n,m) points on the grid) . Please note that  $y_b(n,m)$ is not a straight line and is curved in both indices.

We apply the conservation of probability mass till arbitrary nth and mth points on the grid as

$\frac{\partial \Big[\int_0^{w_b(n)} \, \int_0^{y_b(n,m)} \, P(w,y,t) \,\, dy \, dw \Big]}{\partial t} \, =0$
$=\int_0^{w_b(n)} \, \int_0^{y_b(n,m)} \, \frac{P(w,y,t)}{\partial t} \, \, dy \, dw \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(w_b(n),y,t) dy$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(w,y_b,t) dw =0$

Substituting the 2D uncorrelated PDE in the first RHS integral, we get the expression

$=\int_0^{w_b(n)} \, \int_0^{y_b(n,m)} \, \Big[\, -\frac{\partial \big[\mu_y \, P(w,y,t) \big]}{\partial y} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(w,y,t) \big]}{\partial y^2} \, -\frac{\partial \big[\mu_w \, P(w,y,t) \big]}{\partial w} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_w}^2 \, P(w,y,t) \big]}{\partial w^2} \, \Big] \, \, dy \, dw \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(w_b(n),y,t) dy$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(w,y_b,t) dw =0$

We do the integrations in the first term to get
$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, P(w,y,t) \, +\, .5 \, \frac{\partial \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(w,y,t) \big]}{\partial y} \, \Big] dw$
$\int_0^{y_b(n,m)} \Big[ -\mu_w \, P(w,y,t) \, +\, .5 \, \frac{\partial \big[{\sigma_w}^2 \, P(w,y,t) \big]}{\partial w} \, \Big] \, \, dy \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(w_b(n),y,t) dy$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(w,y_b,t) dw =0$

Now we convert our coordinates to $Z_1$ and $Z_2$ instead of w and y as

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, \frac{\partial \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial y} \, \Big] dw$
$\int_0^{y_b(n,m)} \Big[ -\mu_w \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, \frac{\partial \big[{\sigma_w}^2 \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial w} \, \Big] \, \, dy \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(Z_1(n)) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} dy$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1(n)}{\partial w} \, \frac{\partial Z_2(m)}{\partial y} dw =0$

which can be further simplified as

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$\int_0^{y_b(n,m)} \Big[ -\mu_w \, \frac{\partial Z_1}{\partial w} \, +\, .5 \, {\sigma_w}^2 \, \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2 - {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \big] \Big] \, \, dZ_2 \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, \frac{\partial Z_1}{\partial w} \, dZ_2$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0$

But we notice that second and third term together are merely the solution of 1D independent PDE integrated in an extra dimension and since the original solution of both terms for independent 1D PDE was zero, the integrated term should also be zero. (Even if it were not zero, we could have found this solution from solving the 1D independent PDE), therefore
$\int_0^{y_b(n,m)} \Big[ -\mu_w \, \frac{\partial Z_1}{\partial w} \, +\, .5 \, {\sigma_w}^2 \, \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2 - {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \big] \Big] \, \, dZ_2 \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, \frac{\partial Z_1}{\partial w} \, dZ_2=0$

Now we are left with the terms to solve as

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0$

For a given value of m, we can recursively solve the above PDE for each value of n. To make matters simple, I write the above equation as

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$+ \int_0^{w_b(n-1)} \frac{\partial y_b(n,m)}{\partial t} \, \frac{\partial Z_2}{\partial y} dZ_1$
$+ \int_{w_b(n-1)}^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, \frac{\partial Z_2}{\partial y} dZ_1 =0$

where in the above equation, I have expanded the last terms as sum of all previous terms and one current term to solve and now the equation can be solved as

$\frac{\partial y_b(n,m)}{\partial t} =\Bigg[ \int_0^{w_b(n)} \, \Big[\, \mu_y \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, \big[ Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 + {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$- \int_0^{w_b(n-1)} \frac{\partial y_b(n,m)}{\partial t} \, \frac{\partial Z_2}{\partial y} dZ_1 \Bigg] / \Big[ \int_{w_b(n-1)}^{w_b(n)} \, \frac{\partial Z_2}{\partial y} dZ_1 \Big]$

please note that $\frac{\partial y_b(n,m)}{\partial t} \Delta t$ is the change in coordinates y(n,m) so as to conserve the mass in a small time step according to given PDE dynamics.
The last equation solves the 2D PDE in the asset dimension. Adding correlation term should not be difficult but then we have to make a few changes as correlation will appear in jacobian. Also derivatives of asset with respect to $Z_1$ would have to be included something that was making my equation unstable when I tried it last time but I will try it again to see if I could make it work.
I will be posting the program for uncorrelated PDE based on these analytics tomorrow.
.
.
.
Sorry friends, I made a mistake. I took both $P(Z_1) \, P(Z_2)$ common out of equations which is not possible due to integrations. Sorry that I did these calculations while writing while I have not taken them common in my calculations that I implemented in matlab. Here are modified equations

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, P(Z_2) \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$\int_0^{y_b(n,m)} \Big[ -\mu_w \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, +\, .5 \, {\sigma_w}^2 \, P(Z_1) \, P(Z_2) \, \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2 - {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \big] \Big] \, \, dZ_2 \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, dZ_2$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0$

independent PDE for SV still goes to zero, leaving us with

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, P(Z_2) \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0$

We solve the above equation for a given value of m recursively over different values of n so $P(Z_2)$ can be taken common in the above equation. Please note that we can eliminate P(Z_2) which is constant for a given index m even when index n changes but we cannot eliminate $\frac{\partial Z_2}{\partial y}$ since we have curved lines and for given index m when n is different the value of $\frac{\partial Z_2}{\partial y}$ would be different. so it gives us

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, P(Z_1) \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0$

which takes us to final solution
$\frac{\partial y_b(n,m)}{\partial t} =\Bigg[ \int_0^{w_b(n)} \, P(Z_1) \, \Big[\, \mu_y \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, \big[ Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 + {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$- \int_0^{w_b(n-1)} \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, \frac{\partial Z_2}{\partial y} dZ_1 \Bigg] / \Big[ \int_{w_b(n-1)}^{w_b(n)} \, P(Z_1) \, \frac{\partial Z_2}{\partial y} dZ_1 \Big]$

Amin
Topic Author
Posts: 2894
Joined: July 14th, 2002, 3:00 am

### Re: Breakthrough in the theory of stochastic differential equations and their simulation

I am going to share the analytics for solution of 2D stochastic volatility equations. Before I give the solution, I want to mention that in typical SV system of SDEs, there are two SDEs out of which one (usually Stochastic volatility) is independent of the other SDE and its one dimensional PDE can be solved independently. If we had a system of two SDEs which are both dependent on each other, in that case the solution though still possible would be more difficult.

I have taken both asset and SV Sdes in bessel coordinates. SV SDE can be solved independently and its PDE is given as

$\frac{\partial P(w,t)}{\partial t} \, = \, -\frac{\partial [\mu_w \, P(w,t)]}{\partial w} \, +\, .5 \, \frac{\partial^2 [{\sigma_w}^2 \, P(w,t)]}{\partial w^2} \,$

When we solve this equation with conservation of mass as I showed in a previous post, we get

$-\mu_w \, \frac{\partial Z_1}{\partial w} - .5 {\sigma_w}^2 \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2- P(Z_1)\, {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \, \big] \,$
$+ \frac{\partial w_b}{\partial t} \frac{\partial Z_1}{\partial w} =0$

In the above equation b in $w_b$ denotes that it is the boundary point that is moving so as to keep the probability conserved up till that point. We will need the above equation because when we move to two dimensional stochastic volatility PDE, the whole of above equation would be integrated in another dimension but since it equates to zero, it will also remain zero even after an extra integration and we will eliminate it from 2D equation.

our two dimensional SV and asset partial differential equation in bessel coordinates that we want to solve is
$\frac{\partial P(w,y,t)}{\partial t} \, = \, -\frac{\partial \big[\mu_y \, P(w,y,t) \big]}{\partial y} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(w,y,t) \big]}{\partial y^2} \, -\frac{\partial \big[\mu_w \, P(w,y,t) \big]}{\partial w} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_w}^2 \, P(w,y,t) \big]}{\partial w^2} \,+ \,\rho \, \frac{\partial^2 \big[ \sigma_y \, \sigma_w \, {V(w)}^{ \gamma_V} \, P(w,y,t) \big]}{\partial w \, \partial y} \,$

though correlation term is conceptually simple and when we do a two dimensional integration to preserve probability mass, this term becomes even more easy but correlation has some other issues when taking derivatives with respect of density of y with respect to Z_1 in diffusion terms. When there is no correlation our jacobian for two dimensional density when taken with respect to two independent standard nromals becomes

$\frac{\partial (Z_1,Z_2)}{\partial (w,y)}=\frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}$

since $Z_1$ and $Z_2$ are independent, we can write

$P(w,y)=P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}$

Now we set ourselves in the probability conservation set up, and we have to write a two dimensional integral from infinity(or rather zero ) to moving boundary (boundary might be a misnomer since this is not grid boundary but rather boundary of mass conserved till a certain grid point and this so called boundary is moving or alternatively the grid point is moving) points $w_b$ and $y_b$. Please note that on a two dimensional grid there are a series of these points $w_b$'s and $y_b$'s. The grid in inependent SV dimension is given by index n and the grid on dependent asset dimension is given by (n,m) and respective moving boundaries in conservation set up are indicated as  $w_b(n)$(that conserves mass till nth point on SV grid) and $y_b(n,m)$ (that conserves mass till (n,m) points on the grid) . Please note that  $y_b(n,m)$ is not a straight line and is curved in both indices.

We apply the conservation of probability mass till arbitrary nth and mth points on the grid as

$\frac{\partial \Big[\int_0^{w_b(n)} \, \int_0^{y_b(n,m)} \, P(w,y,t) \,\, dy \, dw \Big]}{\partial t} \, =0$
$=\int_0^{w_b(n)} \, \int_0^{y_b(n,m)} \, \frac{P(w,y,t)}{\partial t} \, \, dy \, dw \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(w_b(n),y,t) dy$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(w,y_b,t) dw =0$

Substituting the 2D uncorrelated PDE in the first RHS integral, we get the expression

$=\int_0^{w_b(n)} \, \int_0^{y_b(n,m)} \, \Big[\, -\frac{\partial \big[\mu_y \, P(w,y,t) \big]}{\partial y} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(w,y,t) \big]}{\partial y^2} \, -\frac{\partial \big[\mu_w \, P(w,y,t) \big]}{\partial w} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_w}^2 \, P(w,y,t) \big]}{\partial w^2} \, \Big] \, \, dy \, dw \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(w_b(n),y,t) dy$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(w,y_b,t) dw =0$

We do the integrations in the first term to get
$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, P(w,y,t) \, +\, .5 \, \frac{\partial \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(w,y,t) \big]}{\partial y} \, \Big] dw$
$\int_0^{y_b(n,m)} \Big[ -\mu_w \, P(w,y,t) \, +\, .5 \, \frac{\partial \big[{\sigma_w}^2 \, P(w,y,t) \big]}{\partial w} \, \Big] \, \, dy \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(w_b(n),y,t) dy$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(w,y_b,t) dw =0$

Now we convert our coordinates to $Z_1$ and $Z_2$ instead of w and y as

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, \frac{\partial \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial y} \, \Big] dw$
$\int_0^{y_b(n,m)} \Big[ -\mu_w \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, \frac{\partial \big[{\sigma_w}^2 \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial w} \, \Big] \, \, dy \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(Z_1(n)) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} dy$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1(n)}{\partial w} \, \frac{\partial Z_2(m)}{\partial y} dw =0$

which can be further simplified as

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$\int_0^{y_b(n,m)} \Big[ -\mu_w \, \frac{\partial Z_1}{\partial w} \, +\, .5 \, {\sigma_w}^2 \, \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2 - {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \big] \Big] \, \, dZ_2 \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, \frac{\partial Z_1}{\partial w} \, dZ_2$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0$

But we notice that second and third term together are merely the solution of 1D independent PDE integrated in an extra dimension and since the original solution of both terms for independent 1D PDE was zero, the integrated term should also be zero. (Even if it were not zero, we could have found this solution from solving the 1D independent PDE), therefore
$\int_0^{y_b(n,m)} \Big[ -\mu_w \, \frac{\partial Z_1}{\partial w} \, +\, .5 \, {\sigma_w}^2 \, \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2 - {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \big] \Big] \, \, dZ_2 \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, \frac{\partial Z_1}{\partial w} \, dZ_2=0$

Now we are left with the terms to solve as

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0$

For a given value of m, we can recursively solve the above PDE for each value of n. To make matters simple, I write the above equation as

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$+ \int_0^{w_b(n-1)} \frac{\partial y_b(n,m)}{\partial t} \, \frac{\partial Z_2}{\partial y} dZ_1$
$+ \int_{w_b(n-1)}^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, \frac{\partial Z_2}{\partial y} dZ_1 =0$

where in the above equation, I have expanded the last terms as sum of all previous terms and one current term to solve and now the equation can be solved as

$\frac{\partial y_b(n,m)}{\partial t} =\Bigg[ \int_0^{w_b(n)} \, \Big[\, \mu_y \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, \big[ Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 + {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$- \int_0^{w_b(n-1)} \frac{\partial y_b(n,m)}{\partial t} \, \frac{\partial Z_2}{\partial y} dZ_1 \Bigg] / \Big[ \int_{w_b(n-1)}^{w_b(n)} \, \frac{\partial Z_2}{\partial y} dZ_1 \Big]$

please note that $\frac{\partial y_b(n,m)}{\partial t} \Delta t$ is the change in coordinates y(n,m) so as to conserve the mass in a small time step according to given PDE dynamics.
The last equation solves the 2D PDE in the asset dimension. Adding correlation term should not be difficult but then we have to make a few changes as correlation will appear in jacobian. Also derivatives of asset with respect to $Z_1$ would have to be included something that was making my equation unstable when I tried it last time but I will try it again to see if I could make it work.
I will be posting the program for uncorrelated PDE based on these analytics tomorrow.
.
.
.
Sorry friends, I made a mistake. I took both $P(Z_1) \, P(Z_2)$ common out of equations which is not possible due to integrations. Sorry that I did these calculations while writing while I have not taken them common in my calculations that I implemented in matlab. Here are modified equations

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, P(Z_2) \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$\int_0^{y_b(n,m)} \Big[ -\mu_w \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, +\, .5 \, {\sigma_w}^2 \, P(Z_1) \, P(Z_2) \, \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2 - {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \big] \Big] \, \, dZ_2 \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, dZ_2$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0$

independent PDE for SV still goes to zero, leaving us with

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, P(Z_2) \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0$

We solve the above equation for a given value of m recursively over different values of n so $P(Z_2)$ can be taken common in the above equation. Please note that we can eliminate P(Z_2) which is constant for a given index m even when index n changes but we cannot eliminate $\frac{\partial Z_2}{\partial y}$ since we have curved lines and for given index m when n is different the value of $\frac{\partial Z_2}{\partial y}$ would be different. so it gives us

$=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, P(Z_1) \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0$

which takes us to final solution
$\frac{\partial y_b(n,m)}{\partial t} =\Bigg[ \int_0^{w_b(n)} \, P(Z_1) \, \Big[\, \mu_y \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, \big[ Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 + {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1$
$- \int_0^{w_b(n-1)} \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, \frac{\partial Z_2}{\partial y} dZ_1 \Bigg] / \Big[ \int_{w_b(n-1)}^{w_b(n)} \, P(Z_1) \, \frac{\partial Z_2}{\partial y} dZ_1 \Big]$
.
.
Friends, very sorry. I feel so goofy making corrections again and again.
I wrote that when we have solved for the independent 1D PDE for volatility, it should go to zero when integrated over $dZ_2$ which is a mistake. Below is the contribution from independent volatility PDE to our 2D PDE

$\int_0^{y_b(n,m)} \Big[ -\mu_w \, \frac{\partial Z_1}{\partial w} \, +\, .5 \, {\sigma_w}^2 \, \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2 - {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \big] \Big] \, P(Z_1) \, P(Z_2) \, \, dZ_2 \,$
$+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, \frac{\partial Z_1}{\partial w} \, P(Z_1) \, P(Z_2) \, dZ_2=0$

This contribution would have gone to zero only if we integrated over $dZ_1$ but our independent PDE for volatility with individual values of index n does not equate to zero so we would have to include this effect in our calculations. Since we have already calculated $frac{\partial w_b}{\partial t}$ from independent volatility PDE this is not a problem and we can add this effect.
I just quickly made changes to my matlab implementation and I was really thrilled to see how sharply my implementation matched with monte carlo simulations. There is a marked improvement in quality of fit from my earlier implementations in which I neglected this effect. The new program though still dealing only with uncorrelated asset dynamics is better than my implementation that I posted few days ago. The correct implementation based on proper conservation of mass is better than the heuristic calculations based program that I posted a few days ago and stability is also much better.
I will post my matlab program tomorrow and hope that friends would like it.

Amin
Topic Author
Posts: 2894
Joined: July 14th, 2002, 3:00 am

### Re: Breakthrough in the theory of stochastic differential equations and their simulation

Friends, my previous post is erroneous. The contribution from independent volatility PDE solution to asset solution in joint PDE indeed goes to zero.

I am posting a new program based on analytics I posted yesterday with independent volatility PDE contribution going to zero. The program still deals with uncorrelated case only. It is a much improved program. The new 2D grid size is 32 X 40 only. The grid spacing is the same as old but I have decreased the grid boundaries to roughly Z=+4 to Z=-4  as opposed to Z=+5 to Z=-5 earlier. Larger grid had no special benefit and it was blowing all the time.
I have also written a very good special function to aggregate probability weighted CDF for the asset and this function very faithfully reproduces the assset density. Earlier function I used was very ad hoc but I went ahead and wrote this special function that I really like.

Here is the new progam.
.
function [] = SVSDEAnalyticAndMonteCarlo2DCorrWmtNew08()

% dX(t)=(a X(t)^alpha1 + b X(t)^alpha2) dt + sqrt(1-rho^2) sigma1 V(t)^gammaV X(t)^gammaX dZ1(t) + rho sigma1 V(t)^gammaV X(t)^gammaX dZ2(t)
% dV(t)=(mu1 V(t)^beta1 + mu2 V(t)^beta2)*dt + sigma0 V(t)^gamma dZ2(t)

%%%%%%%%%%%%%
X0=1.0;
alpha1=0;
alpha2=1;
kappaX=0.0;
thetaX=0.750;
a=kappaX*thetaX;
b=-kappaX;
rho=-.7*0;
sigma1=.75;
gammaV=.5;
gammaX=.65;
%%%%%%%%%%%%%
V0=.35;
thetaV=.1;
kappaV=1.50;
mu1=kappaV*thetaV;
mu2=-kappaV;
beta1=0;
beta2=1;
gamma=.995;
sigma0=.950;
%%%%%%%%%%%%

dtM=.03125;
TtM=32;

dt=dtM/16;
Tt=TtM*16;
Tt=256*2;
dtM=dt*32;
TtM=Tt/32;
T=TtM*dtM;
%paths=100000;

% dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
% Nn=50;  % 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;
%
%
% dMm=.25/1;   % Normal density subdivisions width. would change with number of subdivisions
% Mm=40;  % No of normal density subdivisions
% MmMidl=20;%One half density Subdivision left from mid of normal density(low)
% MmMidh=21;%One half density subdivision right from the mid of normal density(high)
% MmMid=4.0;

dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=40;  % No of normal density subdivisions
NnMidl=20;%One half density Subdivision left from mid of normal density(low)
NnMidh=21;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;

dMm=.25/1;   % Normal density subdivisions width. would change with number of subdivisions
Mm=32;  % No of normal density subdivisions
MmMidl=16;%One half density Subdivision left from mid of normal density(low)
MmMidh=17;%One half density subdivision right from the mid of normal density(high)
MmMid=4.0;

w(1:Nn)=V0^(1-gamma)/(1-gamma);
y(1:Nn,1:Mm)=X0^(1-gammaX)/(1-gammaX);

Z(1:Nn)=(((1:Nn)-.5)*dNn-NnMid);

%Z0(1:Mm)=(((1:Mm)-1.50)*dMm-MmMid);
Z0(1:Mm)=(((1:Mm)-.50)*dMm-MmMid);
Z
Z0
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);
ZProb(2:Nn-1)=normcdf(.5*Z(2:Nn-1)+.5*Z(3:Nn),0,1)-normcdf(.5*Z(2:Nn-1)+.5*Z(1:Nn-2),0,1);

ZCProb(1:Nn)=normcdf(Z(1:Nn));
ZProb0(1)=normcdf(Z(1),0,1);
ZProb0(2:Nn)=normcdf(Z(2:Nn),0,1)-normcdf(Z(1:Nn-1),0,1);

Z0Prob(1)=normcdf(.5*Z0(1)+.5*Z0(2),0,1)-normcdf(.5*Z0(1)+.5*Z0(2)-dMm,0,1);
Z0Prob(Mm)=normcdf(.5*Z0(Mm)+.5*Z0(Mm-1)+dMm,0,1)-normcdf(.5*Z0(Mm)+.5*Z0(Mm-1),0,1);
Z0Prob(2:Mm-1)=normcdf(.5*Z0(2:Mm-1)+.5*Z0(3:Mm),0,1)-normcdf(.5*Z0(2:Mm-1)+.5*Z0(1:Mm-2),0,1);

Z0Prob0(1)=normcdf(Z0(1),0,1);
Z0Prob0(2:Mm)=normcdf(Z0(2:Mm),0,1)-normcdf(Z0(1:Mm-1),0,1);

Z0CProb(1:Mm)=normcdf(Z0(1:Mm));

%%%%%%%%%%%%%%%%%%%%%%%%5
wnStart=1;
ynStart=1;
dwdZ(wnStart:Nn)=0.0;
tic
for tt=1:Tt

V(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
X(wnStart:Nn,ynStart:Mm)=((1-gammaX)*y(wnStart:Nn,ynStart:Mm)).^(1/(1-gammaX));

wMu0dt(wnStart:Nn)= (mu1*V(wnStart:Nn).^(beta1-gamma)+mu2*V(wnStart:Nn).^(beta2-gamma) - ...
.5*gamma*sigma0^2*V(wnStart:Nn).^(gamma-1))*dt + ...
(mu1*(beta1-gamma).*V(wnStart:Nn).^(beta1-gamma-1)+mu2*(beta2-gamma)*V(wnStart:Nn).^(beta2-gamma-1) - ...
.5*gamma*sigma0^2*(gamma-1)*V(wnStart:Nn).^(gamma-2)).* ...
(mu1*V(wnStart:Nn).^(beta1-gamma)+mu2*V(wnStart:Nn).^(beta2-gamma))*dt^2/2+ ...
(mu1*(beta1-gamma).*(beta1-gamma-1).*V(wnStart:Nn).^(beta1-gamma-2)+mu2*(beta2-gamma-1).*(beta2-gamma)*V(wnStart:Nn).^(beta2-gamma-2) - ...
.5*gamma*sigma0^2*(gamma-1)*(gamma-2).*V(wnStart:Nn).^(gamma-3)).* ...
.5.*sigma0^2.*V(wnStart:Nn).^(2*gamma)*dt^2/2;
for nn=wnStart:Nn
yMu0dt(nn,ynStart:Mm)= (a*X(nn,ynStart:Mm).^(alpha1-gammaX)+b*X(nn,ynStart:Mm).^(alpha2-gammaX) - ...
.5*gammaX*sigma1^2*V(nn).^(2*gammaV).*X(nn,ynStart:Mm).^(gammaX-1))*dt + ...
(a*(alpha1-gammaX).*X(nn,ynStart:Mm).^(alpha1-gammaX-1)+b*(alpha2-gammaX)*X(nn,ynStart:Mm).^(alpha2-gammaX-1) - ...
.5*gammaX*sigma1^2*V(nn).^(2*gammaV).*(gammaX-1)*X(nn,ynStart:Mm).^(gammaX-2)).* ...
(a*X(nn,ynStart:Mm).^(alpha1)+b*X(nn,ynStart:Mm).^(alpha2))*dt^2/2+ ...
-(.5*gammaX*sigma1^2*(2*gammaV).*V(nn).^(2*gammaV-1).*X(nn,ynStart:Mm).^(gammaX-1)).* ...
(mu1*V(nn).^(beta1)+mu2*V(nn).^(beta2))*dt^2/2 + ...
(a*(alpha1-gammaX).*(alpha1-gammaX-1).*X(nn,ynStart:Mm).^(alpha1-gammaX-2)+b*(alpha2-gammaX-1).*(alpha2-gammaX)*X(nn,ynStart:Mm).^(alpha2-gammaX-2) - ...
.5*gammaX*sigma1^2*(gammaX-1)*(gammaX-2)*V(nn).^(2*gammaV).*X(nn,ynStart:Mm).^(gammaX-3)).* ...
.5*sigma1^2*V(nn).^(2*gammaV).*X(nn,ynStart:Mm).^(2*gammaX)*dt^2/2+ ...
-(.5*gammaX*sigma1^2*(2*gammaV).*(2*gammaV-1).*V(nn).^(2*gammaV-2).*X(nn,ynStart:Mm).^(gammaX-1)).* ...
(sigma0^2*V(nn).^(2*gamma))*dt^2/2;
end
[dwdZ,d2wdZ2] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z);

for nn=wnStart:Nn
yA(ynStart:Mm)=y(nn,ynStart:Mm);
[dydZ0A,d2ydZ02A] = First2Derivatives2ndOrderEqSpacedA(ynStart,Mm,dMm,yA,Z0);
dydZ0(nn,ynStart:Mm)=dydZ0A(ynStart:Mm);
d2ydZ02(nn,ynStart:Mm)=d2ydZ02A(ynStart:Mm);
d2ydZ02(nn,Mm-1:Mm)=0.0;
end

for nn=ynStart:Mm
yA(wnStart:Nn)=y(wnStart:Nn,nn);
[dydZA,d2ydZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,yA,Z);
dydZ(wnStart:Nn,nn)=dydZA(wnStart:Nn);
d2ydZ2(wnStart:Nn,nn)=d2ydZ2A(wnStart:Nn);
%d2ydZ02(nn,Mm-1:Mm)=0.0;
end

if(tt>7)
dwb(wnStart:Nn)= wMu0dt(wnStart:Nn)+.5*sigma0^2*dt*Z(wnStart:Nn).*dwdZ(wnStart:Nn).^(-1)+ ...
.5*sigma0^2*dt*(dwdZ(wnStart:Nn).^(-2)).*d2wdZ2(wnStart:Nn);
w(wnStart:Nn)=w(wnStart:Nn)+ dwb(wnStart:Nn);
dwbdt(wnStart:Nn)=dwb(wnStart:Nn)/dt;
elseif((tt<=7)&&(tt>1))
dwb(wnStart:Nn)= wMu0dt(wnStart:Nn)+.5*sigma0^2*dt*Z(wnStart:Nn).*dwdZ(wnStart:Nn).^(-1);

w(wnStart:Nn)=w(wnStart:Nn)+ dwb(wnStart:Nn);
dwbdt(wnStart:Nn)=dwb(wnStart:Nn)/dt;
elseif((tt==1))
dwb(wnStart:Nn)= wMu0dt(wnStart:Nn)+sigma0*sqrt(dt).*Z(wnStart:Nn);
w(wnStart:Nn)=w(wnStart:Nn)+ dwb(wnStart:Nn);
dwbdt(wnStart:Nn)=dwb(wnStart:Nn)/dt;
end

dwdZPrev(wnStart:Nn)=dwdZ(wnStart:Nn);
d2wdZ2Prev(wnStart:Nn)=d2wdZ2(wnStart:Nn);

[dwdZNew,d2wdZ2New] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z);

%    Ratio(wnStart:Nn,ynStart:Mm)=(dwdZPrev(wnStart:Nn)*dNn+.25*d2wdZ2Prev(wnStart:Nn)*dNn^2)/(dwdZNew(wnStart:Nn)*dNn+.25*d2wdZ2New(wnStart:Nn)*dNn^2);

tt

if(tt>=3)
for nn=wnStart:Nn
if(nn>wnStart)

TMuy(nn,ynStart:Mm)=TMuy(nn-1,ynStart:Mm)+yMu0dt(nn,ynStart:Mm).*normpdf(Z0(ynStart:Mm)).*(.5*dydZ0(nn,ynStart:Mm)+.5*dydZ0(nn-1,ynStart:Mm)).^(-1).*ZProb0(nn);
TDiffy(nn,ynStart:Mm)=TDiffy(nn-1,ynStart:Mm)+.5*sigma1^2* V(nn).^(2*gammaV).*  ...
(-sqrt(1-rho^2).^0.*Z0(ynStart:Mm).*dydZ0(nn,ynStart:Mm).^(-1) ...
-0*rho*-Z(nn).*dydZ(nn,ynStart:Mm).^(-1)- ...
(dydZ0(nn,ynStart:Mm).^(-2)).*d2ydZ02(nn,ynStart:Mm)).*normpdf(Z0(ynStart:Mm)).*(.5*dydZ0(nn,ynStart:Mm)+.5*dydZ0(nn-1,ynStart:Mm)).^(-1).*ZProb0(nn);
Tpydw(nn,ynStart:Mm)=Tpydw(nn-1,ynStart:Mm)+normpdf(Z0(ynStart:Mm)).*(.5*dydZ0(nn,ynStart:Mm)+.5*dydZ0(nn-1,ynStart:Mm)).^(-1).*ZProb0(nn);
else
TMuy(nn,ynStart:Mm)=yMu0dt(nn,ynStart:Mm).*normpdf(Z0(ynStart:Mm)).*dydZ0(nn,ynStart:Mm).^(-1).*ZProb0(nn);
TDiffy(nn,ynStart:Mm)=.5*sigma1^2* V(nn).^(2*gammaV).* ...
(-sqrt(1-rho^2)^0.*Z0(ynStart:Mm).*dydZ0(nn,ynStart:Mm).^(-1) ...
-0*rho*Z(nn).*dydZ(nn,ynStart:Mm).^(-1)- ...
(dydZ0(nn,ynStart:Mm).^(-2)).*d2ydZ02(nn,ynStart:Mm)).*normpdf(Z0(ynStart:Mm)).*dydZ0(nn,ynStart:Mm).^(-1).*ZProb0(nn);
Tpydw(nn,ynStart:Mm)=normpdf(Z0(ynStart:Mm)).*dydZ0(nn,ynStart:Mm).^(-1).*ZProb0(nn);

end

Tcor(nn,ynStart:Mm)=rho*sigma0*sigma1*V(nn).^gammaV.*normpdf(Z0(ynStart:Mm)).*dydZ0(nn,ynStart:Mm).^(-1).* ...
normpdf(Z(nn)).*dwdZ(nn).^(-1);
Ratio(nn,ynStart:Mm)=(dwdZPrev(nn)./dwdZNew(nn));
end
for nn=wnStart:Nn
if(nn==wnStart)
dyb(nn,ynStart:Mm)=-(-TMuy(nn,ynStart:Mm) + TDiffy(nn,ynStart:Mm)*dt + ...
Tcor(nn,ynStart:Mm)*dt)./Tpydw(nn,ynStart:Mm);
Tyb(nn,ynStart:Mm)=dyb(nn,ynStart:Mm).*normpdf(Z0(ynStart:Mm)).*dydZ0(nn,ynStart:Mm).^(-1).*ZProb0(nn);
else
dyb(nn,ynStart:Mm)=-((-TMuy(nn,ynStart:Mm) + TDiffy(nn,ynStart:Mm)*dt + ...
Tcor(nn,ynStart:Mm)*dt +   Tyb(nn-1,ynStart:Mm)))./ ...
( normpdf(Z0(ynStart:Mm)).*(.5*dydZ0(nn,ynStart:Mm)+.5*dydZ0(nn-1,ynStart:Mm)).^(-1).*ZProb0(nn));
Tyb(nn,ynStart:Mm)=Tyb(nn-1,ynStart:Mm)+dyb(nn,ynStart:Mm).*normpdf(Z0(ynStart:Mm)).*(.5*dydZ0(nn,ynStart:Mm)+.5*dydZ0(nn-1,ynStart:Mm)).^(-1).*ZProb0(nn);
end

y(nn,ynStart:Mm)=y(nn,ynStart:Mm)+dyb(nn,ynStart:Mm).*Ratio(nn,ynStart:Mm).^0;%

end

elseif ((tt<=2)&&(tt>1))
for nn=wnStart:Nn
Ratio(nn,ynStart:Mm)=(dwdZPrev(nn)./dwdZNew(nn));
yv(nn,ynStart:Mm)=.5*sigma1^2.*V(nn).^(2*gammaV).*dt.* ...
(Z0(ynStart:Mm).*dydZ0(nn,ynStart:Mm).^(-1)).*Ratio(nn,ynStart:Mm);
end

y(wnStart:Nn,ynStart:Mm)=y(wnStart:Nn,ynStart:Mm)+yMu0dt(wnStart:Nn,ynStart:Mm)+yv(wnStart:Nn,ynStart:Mm);
elseif(tt==1)
for mm=1:Mm
yv(wnStart:Nn,mm)=sigma1*V(:).^(gammaV).*sqrt(dt)*Z0(mm);%- ...
end

y(wnStart:Nn,ynStart:Mm)=y(wnStart:Nn,ynStart:Mm)+yMu0dt(wnStart:Nn,ynStart:Mm)+yv(wnStart:Nn,ynStart:Mm);
end

end
Vw(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
DfVw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1

DfVw(nn) = (Vw(nn + 1) - Vw(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end

pVw(1:Nn)=0;
for nn = wnStart:Nn-1

pVw(nn) = (normpdf(Z(nn),0, 1))/abs(DfVw(nn));
end

minX=10;
maxX=.001;
for mm=ynStart:Mm
%Xy(mm)=sum(((1-gammaX).*y(wnStart:Nn,mm)).^(1/(1-gammaX)).*ZProb(wnStart:Nn));
Xy_limits(mm)=sum(((1-gammaX).*y(38,mm)).^(1/(1-gammaX)));

end
w
Xy_limits
if(Xy_limits(2)<minX)
minX=real(Xy_limits(2));
end
if(Xy_limits(30)>maxX)
maxX=real(Xy_limits(30));
end

minX=minX-.05;
maxX=maxX+1.50;
minX=max(.001,minX);
minX
maxX
Mm0=150;
dX=(maxX-minX)/(Mm0-1);

[X,XCdf] = ConstructCDFOnNewGrid02(y,wnStart,Nn,ynStart,Mm,Mm0,ZProb,Z0,gammaX,minX,maxX)

Zx(2:Mm0-1)=norminv(XCdf(2:Mm0-1));%is calculated at boundaries of cells
%ZX(2:Mm0)=(Zb(2:Mm0)+Zb(1:Mm0-1))/2.0;%is calculated at center of grid cells.
%Zx(2:Mm0-1)=Zb(2:Mm0-1);%is calculated at center of grid cells.
[Zx(1)] = InterpolateOrderN6(4,X(1),X(2),X(3),X(4),X(5),0,0,Zx(2),Zx(3),Zx(4),Zx(5),0,0);
dXdZx(3:Mm0-2)=(X(4:Mm0-1)-X(2:Mm0-3))./(Zx(4:Mm0-1)-Zx(2:Mm0-3))
[dXdZx(2)] = InterpolateOrderN6(4,Zx(2),Zx(3),Zx(4),Zx(5),Zx(6),0,0,dXdZx(3),dXdZx(4),dXdZx(5),dXdZx(6),0,0);
[dXdZx(1)] = InterpolateOrderN6(4,Zx(1),Zx(3),Zx(4),Zx(5),Zx(6),0,0,dXdZx(3),dXdZx(4),dXdZx(5),dXdZx(6),0,0);

Zb(1:Mm0-3)=(Zx(1:Mm0-3)+Zx(2:Mm0-2))/2;

ZbProb(1)=normcdf(Zb(1));
ZbProb(2:Mm0-3)=normcdf(Zb(2:Mm0-3))-normcdf(Zb(1:Mm0-4));

for mm=1:Mm0-10
pX(mm) = (normpdf(Zx(mm),0, 1))/abs(dXdZx(mm));
end

toc
%Zx

%pX
pX(isnan(pX)==1)=0.0;
pX(isinf(pX)==1)=0.0;
%pX
%XProb
XX=X;
plot(XX(1:Mm0-10),XCdf(1:Mm0-10),'r');
str=input('Look at graph');
plot(XX(2:Mm0-10),pX(2:Mm0-10),'r');
str=input('Look at graph');

ItoHermiteSVMean=sum(Vw(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1))
TrueSVMean=thetaV+(V0-thetaV)*exp(-kappaV*dt*Tt)

ItoHermiteAssetMeanX2=sum(X(2:Mm0-3).*ZbProb(2:Mm0-3))
TrueAssetMean=thetaX+(X0-thetaX)*exp(-kappaX*dt*Tt)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%rng(29079137, 'twister')
rng(15898837, 'twister')
paths=100000;
V(1:paths)=V0;  %Original process monte carlo.
X=0.0;
X(1:paths)=X0;
Random1(1:paths)=0;
Random2(1:paths)=0;

for ttM=1:TtM
Random1=randn(size(Random1));
Random2=randn(size(Random2));

X(1:paths)=X(1:paths)+ ...
(a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dtM + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) * sqrt(dtM) + ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*sqrt(dtM) + ...
(a*alpha1* X(1:paths).^(alpha1-1)+b*alpha2* X(1:paths).^(alpha2-1)).* ...
(((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dtM^2/2)+ ...
(sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) *(1-1/sqrt(3)).*dtM^1.5+ ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*(1-1/sqrt(3)).*dtM^1.5))+ ...
.5*(a*alpha1*(alpha1-1)* X(1:paths).^(alpha1-2)+b*alpha2*(alpha2-1).* X(1:paths).^(alpha2-2)).* ...
( sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX)) *dtM^2/2 + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random1(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*(Random1(1:paths).^2-1) * dtM/2 + ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths)*dtM/2)+ ...
.5*sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
(sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random1(1:paths).*1/sqrt(3).*dtM^1.5)+ ...
sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random1(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
sigma0*V(1:paths).^gamma.*Random1(1:paths).*Random2(1:paths)*dtM/2)+ ...
.5*sqrt(1-rho^2)* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
(sigma0^2*V(1:paths).^(2*gamma).*Random1(1:paths)*1/sqrt(3)*dtM^1.5)+ ...
sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random1(1:paths)*1/sqrt(3)*dtM^1.5+ ...
rho* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random2(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths) * dtM/2 + ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*(Random2(1:paths).^2-1)*dtM/2)+ ...
.5*rho* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
(sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random2(1:paths).*1/sqrt(3).*dtM^1.5)+ ...
rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
sigma0*V(1:paths).^gamma.*(Random2(1:paths).^2-1)*dtM/2)+ ...
.5*rho* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
sigma0^2.*V(1:paths).^(2*gamma).*Random2(1:paths) * 1/sqrt(3).* dtM^1.5+ ...
rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random2(1:paths)*1/sqrt(3)*dtM^1.5;

V(1:paths)=V(1:paths)+ ...
(mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dtM + ...
sigma0*V(1:paths).^gamma .*Random2(1:paths)*sqrt(dtM) + ...
(mu1.*beta1*V(1:paths).^(beta1-1) + mu2.*beta2.*V(1:paths).^(beta2-1)).* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dtM^2/2 + ...
sigma0*V(1:paths).^gamma .*Random2(1:paths)*(1-1/sqrt(3))*dtM^1.5) + ...
.5*(mu1.*beta1.*(beta1-1).*V(1:paths).^(beta1-2) + mu2.*beta2.*(beta2-1).*V(1:paths).^(beta2-2)).* ...
sigma0^2.*V(1:paths).^(2*gamma).*dtM^2/2 + ...
sigma0*gamma*V(1:paths).^(gamma-1) .* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths).*1/sqrt(3)*dtM^1.5 + ...
sigma0.*V(1:paths).^gamma .*(Random2(1:paths).^2-1)*dtM/2) + ...
.5*sigma0*gamma*(gamma-1).*V(1:paths).^(gamma-2) .* ...
sigma0^2.*V(1:paths).^(2*gamma) .*Random2(1:paths).*1/sqrt(3)*dtM^1.5;

end

SVolMeanAnalytic=thetaV+(V0-thetaV)*exp(-kappaV*dt*Tt)
SVolMeanMC=sum(V(1:paths))/paths
AssetMeanAnalytic=thetaX+exp(-kappaX*T)*(X0-thetaX)
AssetMeanMC=sum(X(1:paths))/paths

MaxCutOff=30;
%NoOfBins=round(600*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too
NoOfBins=500;
[XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti_NEW(X,paths,NoOfBins,MaxCutOff );
plot(IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g',XX(1:Mm0-10),pX(1:Mm0-10),'r');
%plot(IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g',XX(3:80),pX(3:80),'r');
%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(sprintf('X0 = %.4f,thetaX=%.3f,kappaX=%.2f,gammaX=%.3f,gammaV=%.3f,rho=%.2f,sigma1=%.2f,T=%.2f,dt=%.5f,AMean=%.4f,McMean=%.4f,TMean=%.4f', X0,thetaX,kappaX,gammaX,gammaV,rho,sigma1,T,dt,ItoHermiteAssetMeanX2,AssetMeanMC,AssetMeanAnalytic));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

%legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')

str=input('Green line is density of Asset SDE.');

MaxCutOff=30;
%NoOfBins=round(600*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too
NoOfBins=1000;
[VDensity,IndexOutV,IndexMaxV] = MakeDensityFromSimulation_Infiniti_NEW(V,paths,NoOfBins,MaxCutOff );
plot(IndexOutV(1:IndexMaxV),VDensity(1:IndexMaxV),'g',Vw(wnStart+1:Nn-1),pVw(wnStart+1:Nn-1),'r');
%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');
%plot(v(wnStart+1:Nn-1),pv(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');

title(sprintf('V0 = %.4f,thetaV=%.3f,kappaV=%.2f,gamma=%.3f,sigma0=%.2f,T=%.2f,dt=%.5f,SVolMean=%.4f,SVolMeanAnalytic=%.4f', V0,thetaV,kappaV,gamma,sigma0,T,dt,SVolMeanMC,SVolMeanAnalytic));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

%legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')

str=input('red line is density of Asset SDE.');

end


.
.
Here is the new sub-function
.
function [X,XCdf] = ConstructCDFOnNewGrid02(y,wnStart,Nn,ynStart,Mm,Mm0,ZProb,Z1,gammaX,minX,maxX)

Xyy(wnStart:Nn,ynStart:Mm)=((1-gammaX).*y(wnStart:Nn,ynStart:Mm)).^(1/(1-gammaX));

% maxX
% minX
% str=input('Look at numbers');

X(1)=minX;
X(Mm0)=maxX;
dX=(X(Mm0)-X(1))/(Mm0-1);
X(1:Mm0)=X(1)+(0:Mm0-1)*dX;
%Xa(1:Mm0)=X(1:Mm0)-dX/2;
%Xb(1:Mm0)=X(1:Mm0)+dX/2;
%Xa(1)=X(1)-dX/20;
%if(Xa(1)<0)
%    Xa(1)=.0000001;
%end
XCdf(1:Mm0)=0.0;

ZProb(wnStart+2:Nn-2)=ZProb(wnStart+2:Nn-2)/sum(ZProb(wnStart+2:Nn-2));

for nn=wnStart+2:Nn-2
Xy(ynStart:Mm)=real(Xyy(nn,ynStart:Mm));
%str=input('Look at Xy');
mm=2;
mm0=1;
while(mm0<Mm0-2)
if((mm==2))
while((X(mm0)<=Xy(mm)) &&(mm0<Mm0))
Zx(mm0) =InterpolateOrderN6(4,X(mm0),Xy(mm),Xy(mm+1),Xy(mm+2),Xy(mm+3),0,0,Z1(mm),Z1(mm+1),Z1(mm+2),Z1(mm+3),0,0);
XCdf(mm0)=XCdf(mm0)+normcdf(real(Zx(mm0))).*ZProb(nn);
mm0=mm0+1;
end
mm=mm+1;
elseif( (mm==3))
while((X(mm0)<=Xy(mm))&&(mm0<Mm0))
Zx(mm0) =InterpolateOrderN6(5,X(mm0),Xy(mm-1),Xy(mm),Xy(mm+1),Xy(mm+2),Xy(mm+3),0,Z1(mm-1),Z1(mm),Z1(mm+1),Z1(mm+2),Z1(mm+3),0);
if((mm0>1)&&(Zx(mm0)<Zx(mm0-1)))
Zx(mm0)=Zx(mm0-1);
end
XCdf(mm0)=XCdf(mm0)+normcdf(real(Zx(mm0))).*ZProb(nn);
mm0=mm0+1;
end
mm=mm+1;
elseif( (mm==4))
while((X(mm0)<=Xy(mm))&&(mm0<Mm0))
Zx(mm0) =InterpolateOrderN6(5,X(mm0),Xy(mm-2),Xy(mm-1),Xy(mm),Xy(mm+1),Xy(mm+2),0,Z1(mm-2),Z1(mm-1),Z1(mm),Z1(mm+1),Z1(mm+2),0);
if((mm0>1)&&(Zx(mm0)<Zx(mm0-1)))
Zx(mm0)=Zx(mm0-1);
end
XCdf(mm0)=XCdf(mm0)+normcdf(real(Zx(mm0))).*ZProb(nn);
mm0=mm0+1;
end
mm=mm+1;
elseif( (mm>4) && (mm<=Mm-3))
while((X(mm0)<=Xy(mm))&&(mm0<Mm0))
Zx(mm0) =InterpolateOrderN6(6,X(mm0),Xy(mm-3),Xy(mm-2),Xy(mm-1),Xy(mm),Xy(mm+1),Xy(mm+2),Z1(mm-3),Z1(mm-2),Z1(mm-1),Z1(mm),Z1(mm+1),Z1(mm+2));
if((mm0>1)&&(Zx(mm0)<Zx(mm0-1)))
Zx(mm0)=Zx(mm0-1);
end
XCdf(mm0)=XCdf(mm0)+normcdf(real(Zx(mm0))).*ZProb(nn);
mm0=mm0+1;
end
mm=mm+1;
elseif(mm==Mm-2)
while((X(mm0)<=Xy(mm))&&(mm0<Mm0))
Zx(mm0) =InterpolateOrderN6(5,X(mm0),Xy(mm-3),Xy(mm-2),Xy(mm-1),Xy(mm),Xy(mm+1),0,Z1(mm-3),Z1(mm-2),Z1(mm-1),Z1(mm),Z1(mm+1),0);
XCdf(mm0)=XCdf(mm0)+normcdf(real(Zx(mm0))).*ZProb(nn);
if((mm0>1)&&(Zx(mm0)<Zx(mm0-1)))
Zx(mm0)=Zx(mm0-1);
end
mm0=mm0+1;
end
mm=mm+1;
elseif(mm==Mm-1)
while((X(mm0)<=Xy(mm))&&(mm0<Mm0))
Zx(mm0) =InterpolateOrderN6(4,X(mm0),Xy(mm-3),Xy(mm-2),Xy(mm-1),Xy(mm),0,0,Z1(mm-3),Z1(mm-2),Z1(mm-1),Z1(mm),0,0);
if((mm0>1)&&(Zx(mm0)<Zx(mm0-1)))
Zx(mm0)=Zx(mm0-1);
end
XCdf(mm0)=XCdf(mm0)+normcdf(real(Zx(mm0))).*ZProb(nn);
mm0=mm0+1;
end
while((X(mm0)>Xy(mm))&&(mm0<Mm0))
Zx(mm0) =InterpolateOrderN6(4,X(mm0),Xy(mm-3),Xy(mm-2),Xy(mm-1),Xy(mm),0,0,Z1(mm-3),Z1(mm-2),Z1(mm-1),Z1(mm),0,0);
if((mm0>1)&&(Zx(mm0)<Zx(mm0-1)))
Zx(mm0)=Zx(mm0-1);
end
XCdf(mm0)=XCdf(mm0)+normcdf(real(Zx(mm0))).*ZProb(nn);
mm0=mm0+1;
end
end
end

end

end

.
You will also need another minor function for interpolation.
.
function [y0] = InterpolateOrderN6(N,x0,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6)

if(N==6)
y0=(x0-x2).*(x0-x3).*(x0-x4).*(x0-x5).*(x0-x6)./((x1-x2).*(x1-x3).*(x1-x4).*(x1-x5).*(x1-x6)).*y1+ ...
(x0-x1).*(x0-x3).*(x0-x4).*(x0-x5).*(x0-x6)./((x2-x1).*(x2-x3).*(x2-x4).*(x2-x5).*(x2-x6)).*y2+ ...
(x0-x1).*(x0-x2).*(x0-x4).*(x0-x5).*(x0-x6)./((x3-x1).*(x3-x2).*(x3-x4).*(x3-x5).*(x3-x6)).*y3+ ...
(x0-x1).*(x0-x2).*(x0-x3).*(x0-x5).*(x0-x6)./((x4-x1).*(x4-x2).*(x4-x3).*(x4-x5).*(x4-x6)).*y4+ ...
(x0-x1).*(x0-x2).*(x0-x3).*(x0-x4).*(x0-x6)./((x5-x1).*(x5-x2).*(x5-x3).*(x5-x4).*(x5-x6)).*y5+ ...
(x0-x1).*(x0-x2).*(x0-x3).*(x0-x4).*(x0-x5)./((x6-x1).*(x6-x2).*(x6-x3).*(x6-x4).*(x6-x5)).*y6;
end

if(N==5)
y0=(x0-x2).*(x0-x3).*(x0-x4).*(x0-x5)./((x1-x2).*(x1-x3).*(x1-x4).*(x1-x5)).*y1+ ...
(x0-x1).*(x0-x3).*(x0-x4).*(x0-x5)./((x2-x1).*(x2-x3).*(x2-x4).*(x2-x5)).*y2+ ...
(x0-x1).*(x0-x2).*(x0-x4).*(x0-x5)./((x3-x1).*(x3-x2).*(x3-x4).*(x3-x5)).*y3+ ...
(x0-x1).*(x0-x2).*(x0-x3).*(x0-x5)./((x4-x1).*(x4-x2).*(x4-x3).*(x4-x5)).*y4+ ...
(x0-x1).*(x0-x2).*(x0-x3).*(x0-x4)./((x5-x1).*(x5-x2).*(x5-x3).*(x5-x4)).*y5;
end
if(N==4)
y0=(x0-x2).*(x0-x3).*(x0-x4)./((x1-x2).*(x1-x3).*(x1-x4)).*y1+ ...
(x0-x1).*(x0-x3).*(x0-x4)./((x2-x1).*(x2-x3).*(x2-x4)).*y2+ ...
(x0-x1).*(x0-x2).*(x0-x4)./((x3-x1).*(x3-x2).*(x3-x4)).*y3+ ...
(x0-x1).*(x0-x2).*(x0-x3)./((x4-x1).*(x4-x2).*(x4-x3)).*y4;
end
if(N==3)
y0=(x0-x2).*(x0-x3)./((x1-x2).*(x1-x3)).*y1+ ...
(x0-x1).*(x0-x3)./((x2-x1).*(x2-x3)).*y2+ ...
(x0-x1).*(x0-x2)./((x3-x1).*(x3-x2)).*y3;
end
if(N==2)
y0=(x0-x2)./((x1-x2)).*y1+ ...
(x0-x1)./((x2-x1)).*y2;
end

end
Here is the output of the program

ItoHermiteSVMean =
0.155813360436568
TrueSVMean =
0.155782540037107
ItoHermiteAssetMeanX2 =
0.999871301686896
TrueAssetMean =
1
SVolMeanAnalytic =
0.155782540037107
SVolMeanMC =
0.156372408825739
AssetMeanAnalytic =
1
AssetMeanMC =
1.000169980402764
IndexMax =
501
Green line is density of Asset SDE.

Here are the graphs that you will see when you run this program.