Friends, I have caught malaria and have very high fever. I will post the complete program after I recover. Here is the tentative program and I wanted to change a lot of things in it to make it easy to understand but now I am posting as such. I will post a reworked version once I recover.

.

.

```
function [] = FPERevisitedTransProb08ABwNew04DownloadedAB1()
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=mu1 x(t)^beta1 dt + mu2 x(t)^beta2 dt +sigma x(t)^gamma dz(t)
%I have not directly simulated the SDE but simulated the transformed
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.
dt=.125/16/2/2; % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=128*2*2*2; % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
T=Tt*dt;
OrderA=4; %
OrderM=4; %
dtM=.0625/2;
TtM=T/dtM;
dNn=.2/1; % Normal density subdivisions width. would change with number of subdivisions
Nn=41; % No of normal density subdivisions
NnMidl=18;%One half density Subdivision left from mid of normal density(low)
NnMidh=19;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;
NnT=46;
%NnD=(NnT-Nn)/2;
NnD=5;
x0=.25; % starting value of SDE
beta1=0.0;
beta2=1.0; % Second drift term power.
gamma=.75;%50; % volatility power.
kappa=1.0;%.950; %mean reversion parameter.
theta=.06;%mean reversion target
sigma0=.9;%Volatility value
%you can specify any general mu1 and mu2 and beta1 and beta2.
mu1=1*theta*kappa; %first drift coefficient.
mu2=-1*kappa; % Second drift coefficient.
%mu1=0;
%mu2=0;
alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;%This is for expansion of integrals for calculation of drift
%and volatility coefficients
yy(1:Nn)=x0;
w(1:Nn)=x0^(1-gamma)/(1-gamma);
%Z(1:Nn)=(((1:Nn)-5.5)*dNn-NnMid);
Z(1:Nn)=(((1:Nn)+1.5)*dNn-NnMid);
ZT(1:NnT)=(((1:NnT)-3.5)*dNn-NnMid);
ZT
Z
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);
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);
ZProbT(1)=normcdf(.5*ZT(1)+.5*ZT(2),0,1)-normcdf(.5*ZT(1)+.5*ZT(2)-dNn,0,1);
ZProbT(NnT)=normcdf(.5*ZT(NnT)+.5*ZT(NnT-1)+dNn,0,1)-normcdf(.5*ZT(NnT)+.5*ZT(NnT-1),0,1);
ZProbT(2:NnT-1)=normcdf(.5*ZT(2:NnT-1)+.5*ZT(3:NnT),0,1)-normcdf(.5*ZT(2:NnT-1)+.5*ZT(1:NnT-2),0,1);
%Above calculate probability mass in each probability subdivision.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
sigma11(1:OrderA+1)=0;
mu11(1:OrderA+1)=0;
mu22(1:OrderA+1)=0;
sigma22(1:OrderA+1)=0;
% index 1 correponds to zero level since matlab indexing starts at one.
sigma11(1)=1;
mu11(1)=1;
mu22(1)=1;
sigma22(1)=1;
for k=1:(OrderA+1)
if sigma0~=0
sigma11(k)=sigma0^(k-1);
end
if mu1 ~= 0
mu11(k)=mu1^(k-1);
end
if mu2 ~= 0
mu22(k)=mu2^(k-1);
end
if sigma0~=0
sigma22(k)=sigma0^(2*(k-1));
end
end
%Ft(1:TtM+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0; %General time powers on hermite polynomials
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers on coefficients of hermite polynomials.
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers for bessel transformed coordinates.
%YCoeff0 and YCoeff are coefficents for original coordinates monte carlo.
%YqCoeff0 and YqCoeff are bessel/lamperti version monte carlo.
YCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
YqCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
%Pre-compute the time and power exponent values in small multi-dimensional arrays
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma); %expand y^alpha where alpha=1;
YqCoeff = ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);%expand y^alpha1 where alpha1=(1-gamma)
YqCoeff=YqCoeff/(1-gamma); %Transformed coordinates coefficients have to be
%further divided by (1-gamma)
for k = 0 : (OrderA)
for m = 0:k
l4 = k - m + 1;
for n = 0 : m
l3 = m - n + 1;
for j = 0:n
l2 = n - j + 1;
l1 = j + 1;
%Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
- (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
- (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
YCoeff0(l1,l2,l3,l4) =YCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
YqCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
end
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wnStart=1;%
tic
Zt1(wnStart:Nn)=0.0;
Zt2(wnStart:Nn)=0.0;
for tt=1:Tt
%[wMu0dt,dwMu0dtdw,c1] = CalculateDriftAndVolA4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
%[wMu0dt,c1,c2] = CalculateDriftAndVolA404AB(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
[wMu0dt,c1,c2,dwMu0dtdw] = CalculateDriftAndVolA404AB(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
[wMid] = InterpolateOrderN8(8,0,Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidl),Z(NnMidh),Z(NnMidh+1),Z(NnMidh+2),Z(NnMidh+3),w(NnMidl-3),w(NnMidl-2),w(NnMidl-1),w(NnMidl),w(NnMidh),w(NnMidh+1),w(NnMidh+2),w(NnMidh+3));
Zt1(wnStart:Nn)=w(wnStart:Nn)-wMid;
[dwdZ,d2wdZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt1,Z);
[dwdZ,d2wdZ2,d3wdZ3] = First3Derivatives4thOrderEqSpaced(wnStart,Nn,dNn,Zt1,Z);
tt
C0(wnStart:Nn)=Zt1(wnStart:Nn)-dwdZ(wnStart:Nn).*Z(wnStart:Nn);
%c1(wnStart:Nn)=sigma0*sqrt(dt);
Zt2(wnStart:Nn)=C0(wnStart:Nn)+sign(dwdZ(wnStart:Nn)+c1(wnStart:Nn)).* ...
abs(sqrt(sign(dwdZ(wnStart:Nn)).*(dwdZ(wnStart:Nn)).^2+ ...
sign(c1(wnStart:Nn)).*c1(wnStart:Nn).^2)).*Z(wnStart:Nn);
[dwdZ,d2wdZ2,d3wdZ3] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt2,Z);
[dw2dZ,d2w2dZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt2,Z);
[dw2dZ,d2w2dZ2A,d3wdZ3] = First3Derivatives4thOrderEqSpaced(wnStart,Nn,dNn,Zt2,Z);
if(tt>10)
% dwdt(wnStart:Nn)=(wMu0dt(wnStart:Nn)+Zt2(wnStart:Nn)-Zt1(wnStart:Nn)+ ...
% .5*sigma0^2*(dw2dZ(wnStart:Nn).^(-2)).*d2w2dZ2A(wnStart:Nn).*dt)/dt;
else
% dwdt(wnStart:Nn)=(wMu0dt(wnStart:Nn)+Zt2(wnStart:Nn)-Zt1(wnStart:Nn))/dt;%+ ...
end
wnStart
if(tt>10)
dZdw(wnStart:Nn)=1./dw2dZ(wnStart:Nn);
d2Zdw2(wnStart:Nn)=-d2w2dZ2A(wnStart:Nn).*dZdw(wnStart:Nn).^3;
d3Zdw3(wnStart:Nn)=-d3wdZ3(wnStart:Nn).*dZdw(wnStart:Nn).^4+3*d2w2dZ2A(wnStart:Nn).^2.*dZdw(wnStart:Nn).^5;
%w(wnStart:Nn)=w(wnStart:Nn)+ ...
% wMu0dt(wnStart:Nn)+ ...
% 0*dwMu0dtdw(wnStart:Nn)./(-Z(wnStart:Nn).*dZdw(wnStart:Nn).^2+d2Zdw2(wnStart:Nn)) - ...
% +.5*sigma0^2*dt*( (Z(wnStart:Nn).^2-1).*dZdw(wnStart:Nn).^3- ...
% 0* 3*Z(wnStart:Nn).*dZdw(wnStart:Nn).*d2Zdw2(wnStart:Nn))./(-Z(wnStart:Nn).*dZdw(wnStart:Nn).^2+d2Zdw2(wnStart:Nn));
w(wnStart:Nn)=w(wnStart:Nn)+ ...
wMu0dt(wnStart:Nn)+1*(Zt2(wnStart:Nn)-Zt1(wnStart:Nn))+ ...
1*.5*sigma0^2*(dw2dZ(wnStart:Nn).^(-2)).*d2w2dZ2A(wnStart:Nn).*dt;
else
w(wnStart:Nn)=w(wnStart:Nn)+ ...
wMu0dt(wnStart:Nn)+Zt2(wnStart:Nn)-Zt1(wnStart:Nn);
end
% %yy(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
% [dwdZ,d2wdZ2A,d3wdZ3] = First3Derivatives4thOrderEqSpaced(wnStart,Nn,dNn,w,Z);
%FirstCoordinate
% for nn=NnMidl-11:-1:wnStart+1
% if((w(nn)-w(nn-1))>(w(nn+1)-w(nn)))
%
% w(nn) = InterpolateOrderN8(6,Z(nn),Z(nn+1),Z(nn+2),Z(nn+3),Z(nn+4),Z(nn+5),Z(nn+6),Z(nn+7),Z(nn+8),w(nn+1),w(nn+2),w(nn+3),w(nn+4),w(nn+5),w(nn+6),w(nn+7),w(nn+8));
%
% end
% end
% if(w(wnStart+1)-w(wnStart)>w(wnStart+2)-w(wnStart+1))
% nn=wnStart;
% w(wnStart) = InterpolateOrderN8(6,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),Z(wnStart+8),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7),w(wnStart+8));
% w(nn)=w(nn+1)-(w(nn+2)-w(nn+1));
% end
[wE] = InterpolateOrderN6(6,Z(Nn)+dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));
for nn=wnStart:Nn-1
if(w(nn)>w(nn+1))
w(nn)=0;
end
end
% w1(wnStart:Nn-1)=w(wnStart:Nn-1);
% w1(Nn)=w(Nn);
% w2(wnStart:Nn-1)=w(wnStart+1:Nn);
% w2(Nn)=wE;
% w(w1(wnStart:Nn)>w2(:))=0;%Be careful;might not universally hold;
% %Change 3:7/25/2020: I have improved zero correction in above.
%for nn=1:Nn
% if(w2(nn)>w1(nn))
w(w<0)=0.0;
StartChangeFlag=0;
nnChange=0;
wnStartPrev=wnStart;
for nn=wnStart:Nn
if(w(nn)<=0)
wnStart=nn+1;
StartChangeFlag=1;
nnChange=nnChange+1;
end
end
%nn=wnStart;
%while(w(nn)<=0)
% wnStart=nn+1;
% nn=nn+1;
%end
InterpolatePrevFlag=1;
wnChange=wnStart-wnStartPrev;
Zs=0
ws=0;
for mm=1:nnChange
Zs(mm)=Z(wnStart)-mm*dNn;
ws(mm)=0.0
end
if(nnChange>=1)
ws = InterpolateOrderN8(8,Zs,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
ContinueFlag=1;
end
for mm=1:nnChange
if((ws(mm)>0)&&(ws(mm)<w(wnStart))&&(ContinueFlag==1))
wnStart=wnStart-1;
w(wnStart)=ws(mm);
ContinueFlag=1;
else
ContinueFlag=0;
end
end
% if(InterpolatePrevFlag==1)
% InterpolatePrevFlag=0;
% %wStart1 = InterpolateOrderN8(8,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
% NnStart=nnChange
% wStart1 = InterpolateOrderN8(8,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
% if(wStart1>0)
% wnStart=wnStart-1;
% w(wnStart)=wStart1;
% InterpolatePrevFlag=1;
% end
% if(wnStart==wnStartPrev)
% % InterpolatePrevFlag=0;
% end
% end
% end
% if(nnChange>=1)
% wStart1 = InterpolateOrderN8(8,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
% if(wStart1>0)
% wnStart=wnStart-1
% w(wnStart)=wStart1;
% InterpolatePrevFlag=1;
% end
% end
% if((nnChange>=2)&&( InterpolatePrevFlag==1))
% wStart1 = InterpolateOrderN8(8,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
% InterpolatePrevFlag=0;
% if(wStart1>0)
% wnStart=wnStart-1
% w(wnStart)=wStart1;
% InterpolatePrevFlag=1;
% end
% end
% if((nnChange>=3)&&( InterpolatePrevFlag==1))
% wStart1 = InterpolateOrderN8(8,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
% InterpolatePrevFlag=0;
% if(wStart1>0)
% wnStart=wnStart-1
% w(wnStart)=wStart1;
% InterpolatePrevFlag=1;
% end
% end
%
w
wnStart
tt
end
wT(wnStart+NnD:Nn+NnD)=w(wnStart:Nn);
Zs=0
ws=0;
for mm=1:wnStart+NnD-1
Zs(mm)=Z(wnStart)-mm*dNn;
ws(mm)=0.0
end
wnStartT=wnStart+NnD;
%NnT=Nn+2*NnD;
ws = InterpolateOrderN8(8,Zs,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
ContinueFlag=1;
wnStart0=wnStart+NnD-1;
for mm=1:wnStart0
if((ws(mm)>0)&&(ws(mm)<wT(wnStartT))&&(ContinueFlag==1))
wnStartT=wnStartT-1;
wT(wnStartT)=ws(mm);
ContinueFlag=1;
else
ContinueFlag=0;
end
end
% Zs=0;
% ws=0;
% for mm=1:NnD
% Zs(mm)=Z(Nn)+mm*dNn;
% ws(mm)=0.0
% end
%
% ws = InterpolateOrderN8(2,Zs,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
%
%
% %for mm=1:NnD
%
% wT(Nn+NnD+1:Nn+2*NnD)=ws(1:NnD);
% %end
% for mm=1:wnStart0-1
% if((ws(mm)>0)&&(ws(mm)<w(wnStart))&&(ContinueFlag==1))
% wnStart=wnStart-1;
% w(wnStart)=ws(mm);
% ContinueFlag=1;
% else
% ContinueFlag=0;
% end
% end
wnStartT
%str=input('Look at numbers')
yyT(wnStartT:NnT)=((1-gamma).*wT(wnStartT:NnT)).^(1/(1-gamma));
DfyyT(wnStartT:NnT)=0;
for nn=wnStartT+1:NnT-1
DfyyT(nn) = (yyT(nn + 1) - yyT(nn - 1))/(ZT(nn + 1) - ZT(nn - 1));
%Change of variable derivative for densities
end
pyyT(1:Nn)=0;
for nn = wnStartT:NnT-1
pyyT(nn) = (normpdf(ZT(nn),0, 1))/abs(DfyyT(nn));
end
if(tt>=23)
%plot(yy(wnStart:Nn),pyy(wnStart:Nn),'r');
%yy
%str=input('Look at the output graph');
end
toc
ItoHermiteMean=sum(yyT(wnStartT+1:NnT-1).*ZProbT(wnStartT+1:NnT-1)) %Original process average from coordinates
disp('true Mean only applicable to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
theta1=1;
rng(29079137, 'twister')
paths=200000;
YY(1:paths)=x0; %Original process monte carlo.
Random1(1:paths)=0;
for tt=1:TtM
Random1=randn(size(Random1));
HermiteP1(1,1:paths)=1;
HermiteP1(2,1:paths)=Random1(1:paths);
HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;
YY(1:paths)=YY(1:paths) + ...
(YCoeff0(1,1,2,1).*YY(1:paths).^Fp(1,1,2,1)+ ...
YCoeff0(1,2,1,1).*YY(1:paths).^Fp(1,2,1,1)+ ...
YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1))*dtM + ...
(YCoeff0(1,1,3,1).*YY(1:paths).^Fp(1,1,3,1)+ ...
YCoeff0(1,2,2,1).*YY(1:paths).^Fp(1,2,2,1)+ ...
YCoeff0(2,1,2,1).*YY(1:paths).^Fp(2,1,2,1)+ ...
YCoeff0(1,3,1,1).*YY(1:paths).^Fp(1,3,1,1)+ ...
YCoeff0(2,2,1,1).*YY(1:paths).^Fp(2,2,1,1)+ ...
YCoeff0(3,1,1,1).*YY(1:paths).^Fp(3,1,1,1))*dtM^2 + ...
((YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
(YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+ ...
YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5) .*HermiteP1(2,1:paths) + ...
((YCoeff0(1,1,1,3).*YY(1:paths).^Fp(1,1,1,3) *dtM) + ...
(YCoeff0(1,1,2,3).*YY(1:paths).^Fp(1,1,2,3)+ ...
YCoeff0(1,2,1,3).*YY(1:paths).^Fp(1,2,1,3)+ ...
YCoeff0(2,1,1,3).*YY(1:paths).^Fp(2,1,1,3)).*dtM^2).*HermiteP1(3,1:paths) + ...
((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )).*HermiteP1(4,1:paths) + ...
(YCoeff0(1,1,1,5).*YY(1:paths).^Fp(1,1,1,5)*dtM^2.0).*HermiteP1(5,1:paths);
end
YY(YY<0)=0;
disp('Original process average from monte carlo');
MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
disp('Original process average from our simulation');
ItoHermiteMean=sum(yyT(wnStartT+1:NnT-1).*ZProbT(wnStartT+1:NnT-1)) %Original process average from coordinates
disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
MaxCutOff=30;
NoOfBins=round(300*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
plot(yyT(wnStartT+1:NnT-1),pyyT(wnStartT+1:NnT-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');
title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,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 SDE from Ito-Hermite method, green is monte carlo.');
end
```

.

.

```
function [wMu0dt,c1,c2,dwMu0dtdw] = CalculateDriftAndVolA404AB(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt)
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
Fp2=Fp1/(1-gamma);
wMu0dt(wnStart:Nn)=(YqCoeff0(1,1,2,1).*yy(wnStart:Nn).^Fp1(1,1,2,1)+ ...
YqCoeff0(1,2,1,1).*yy(wnStart:Nn).^Fp1(1,2,1,1)+ ...
YqCoeff0(2,1,1,1).*yy(wnStart:Nn).^Fp1(2,1,1,1))*dt + ...
(YqCoeff0(1,1,3,1).*yy(wnStart:Nn).^Fp1(1,1,3,1)+ ...
YqCoeff0(1,2,2,1).*yy(wnStart:Nn).^Fp1(1,2,2,1)+ ...
YqCoeff0(2,1,2,1).*yy(wnStart:Nn).^Fp1(2,1,2,1)+ ...
YqCoeff0(1,3,1,1).*yy(wnStart:Nn).^Fp1(1,3,1,1)+ ...
YqCoeff0(2,2,1,1).*yy(wnStart:Nn).^Fp1(2,2,1,1)+ ...
YqCoeff0(3,1,1,1).*yy(wnStart:Nn).^Fp1(3,1,1,1))*dt^2 + ...
(YqCoeff0(1,1,4,1).*yy(wnStart:Nn).^Fp1(1,1,4,1)+ ...
YqCoeff0(1,2,3,1).*yy(wnStart:Nn).^Fp1(1,2,3,1)+ ...
YqCoeff0(2,1,3,1).*yy(wnStart:Nn).^Fp1(2,1,3,1)+ ...
YqCoeff0(1,3,2,1).*yy(wnStart:Nn).^Fp1(1,3,2,1)+ ...
YqCoeff0(2,2,2,1).*yy(wnStart:Nn).^Fp1(2,2,2,1)+ ...
YqCoeff0(3,1,2,1).*yy(wnStart:Nn).^Fp1(3,1,2,1)+ ...
YqCoeff0(1,4,1,1).*yy(wnStart:Nn).^Fp1(1,4,1,1)+ ...
YqCoeff0(2,3,1,1).*yy(wnStart:Nn).^Fp1(2,3,1,1)+ ...
YqCoeff0(3,2,1,1).*yy(wnStart:Nn).^Fp1(3,2,1,1)+ ...
YqCoeff0(4,1,1,1).*yy(wnStart:Nn).^Fp1(4,1,1,1))*dt^3+ ...
(YqCoeff0(1,1,5,1).*yy(wnStart:Nn).^Fp1(1,1,5,1)+ ...
YqCoeff0(1,2,4,1).*yy(wnStart:Nn).^Fp1(1,2,4,1)+ ...
YqCoeff0(2,1,4,1).*yy(wnStart:Nn).^Fp1(2,1,4,1)+ ...
YqCoeff0(1,3,3,1).*yy(wnStart:Nn).^Fp1(1,3,3,1)+ ...
YqCoeff0(2,2,3,1).*yy(wnStart:Nn).^Fp1(2,2,3,1)+ ...
YqCoeff0(3,1,3,1).*yy(wnStart:Nn).^Fp1(3,1,3,1)+ ...
YqCoeff0(1,4,2,1).*yy(wnStart:Nn).^Fp1(1,4,2,1)+ ...
YqCoeff0(2,3,2,1).*yy(wnStart:Nn).^Fp1(2,3,2,1)+ ...
YqCoeff0(3,2,2,1).*yy(wnStart:Nn).^Fp1(3,2,2,1)+ ...
YqCoeff0(4,1,2,1).*yy(wnStart:Nn).^Fp1(4,1,2,1)+ ...
YqCoeff0(1,5,1,1).*yy(wnStart:Nn).^Fp1(1,5,1,1)+ ...
YqCoeff0(2,4,1,1).*yy(wnStart:Nn).^Fp1(2,4,1,1)+ ...
YqCoeff0(3,3,1,1).*yy(wnStart:Nn).^Fp1(3,3,1,1)+ ...
YqCoeff0(4,2,1,1).*yy(wnStart:Nn).^Fp1(4,2,1,1)+ ...
YqCoeff0(5,1,1,1).*yy(wnStart:Nn).^Fp1(5,1,1,1))*dt^4;
dwMu0dtdw(wnStart:Nn)=(YqCoeff0(1,1,2,1).*Fp1(1,1,2,1).*((1-gamma)*w(wnStart:Nn)).^(Fp2(1,1,2,1)-1)+ ...
YqCoeff0(1,2,1,1).*Fp1(1,2,1,1).*((1-gamma)*w(wnStart:Nn)).^(Fp2(1,2,1,1)-1)+ ...
YqCoeff0(2,1,1,1).*Fp1(2,1,1,1).*((1-gamma)*w(wnStart:Nn)).^(Fp2(2,1,1,1)-1))*dt ;
%wMu0dtOdt(wnStart:Nn)=(YqCoeff0(1,1,2,1).*yy(wnStart:Nn).^Fp1(1,1,2,1)*0+ ...
% YqCoeff0(1,2,1,1).*yy(wnStart:Nn).^Fp1(1,2,1,1)+ ...
% YqCoeff0(2,1,1,1).*yy(wnStart:Nn).^Fp1(2,1,1,1))*dt ;
c1(wnStart:Nn)=((YqCoeff0(1,1,1,2).*yy(wnStart:Nn).^Fp1(1,1,1,2).*sqrt(dt))+ ...
(YqCoeff0(1,1,2,2).*yy(wnStart:Nn).^Fp1(1,1,2,2)+YqCoeff0(1,2,1,2).*yy(wnStart:Nn).^Fp1(1,2,1,2)+ ...
YqCoeff0(2,1,1,2).*yy(wnStart:Nn).^Fp1(2,1,1,2)).*dt^1.5+ ...
(YqCoeff0(1,1,3,2).*yy(wnStart:Nn).^Fp1(1,1,3,2)+YqCoeff0(1,2,2,2).*yy(wnStart:Nn).^Fp1(1,2,2,2)+ ...
YqCoeff0(2,1,2,2).*yy(wnStart:Nn).^Fp1(2,1,2,2)+YqCoeff0(1,3,1,2).*yy(wnStart:Nn).^Fp1(1,3,1,2)+ ...
YqCoeff0(2,2,1,2).*yy(wnStart:Nn).^Fp1(2,2,1,2)+YqCoeff0(3,1,1,2).*yy(wnStart:Nn).^Fp1(3,1,1,2)).*dt^2.5+ ...
(YqCoeff0(1,1,4,2).*yy(wnStart:Nn).^Fp1(1,1,4,2)+YqCoeff0(1,2,3,2).*yy(wnStart:Nn).^Fp1(1,2,3,2)+ ...
YqCoeff0(2,1,3,2).*yy(wnStart:Nn).^Fp1(2,1,3,2)+YqCoeff0(1,3,2,2).*yy(wnStart:Nn).^Fp1(1,3,2,2)+ ...
YqCoeff0(2,2,2,2).*yy(wnStart:Nn).^Fp1(2,2,2,2)+ YqCoeff0(3,1,2,2).*yy(wnStart:Nn).^Fp1(3,1,2,2)+ ...
YqCoeff0(1,4,1,2).*yy(wnStart:Nn).^Fp1(1,4,1,2)+YqCoeff0(2,3,1,2).*yy(wnStart:Nn).^Fp1(2,3,1,2)+ ...
YqCoeff0(3,2,1,2).*yy(wnStart:Nn).^Fp1(3,2,1,2)+YqCoeff0(4,1,1,2).*yy(wnStart:Nn).^Fp1(4,1,1,2)).*dt^3.5);
c2(wnStart:Nn)=(YqCoeff0(1,1,2,3).*yy(wnStart:Nn).^Fp1(1,1,2,3)+YqCoeff0(1,2,1,3).*yy(wnStart:Nn).^Fp1(1,2,1,3)+ ...
YqCoeff0(2,1,1,3).*yy(wnStart:Nn).^Fp1(2,1,1,3)).*dt^2;%+ ...
(YqCoeff0(1,1,3,3).*yy(wnStart:Nn).^Fp1(1,1,3,3)+YqCoeff0(1,2,2,3).*yy(wnStart:Nn).^Fp1(1,2,2,3)+ ...
YqCoeff0(2,1,2,3).*yy(wnStart:Nn).^Fp1(2,1,2,3) + YqCoeff0(1,3,1,3).*yy(wnStart:Nn).^Fp1(1,3,1,3)+ ...
YqCoeff0(2,2,1,3).*yy(wnStart:Nn).^Fp1(2,2,1,3)+YqCoeff0(3,1,1,3).*yy(wnStart:Nn).^Fp1(3,1,1,3)).*dt^3;
% (YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
% (YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
% YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5+ ...
% (YCoeff0(1,1,3,2).*YY(1:paths).^Fp(1,1,3,2)+YCoeff0(1,2,2,2).*YY(1:paths).^Fp(1,2,2,2)+ ...
% YCoeff0(2,1,2,2).*YY(1:paths).^Fp(2,1,2,2)+YCoeff0(1,3,1,2).*YY(1:paths).^Fp(1,3,1,2)+ ...
% YCoeff0(2,2,1,2).*YY(1:paths).^Fp(2,2,1,2)+YCoeff0(3,1,1,2).*YY(1:paths).^Fp(3,1,1,2)).*dtM^2.5+ ...
% (YCoeff0(1,1,4,2).*YY(1:paths).^Fp(1,1,4,2)+YCoeff0(1,2,3,2).*YY(1:paths).^Fp(1,2,3,2)+ ...
% YCoeff0(2,1,3,2).*YY(1:paths).^Fp(2,1,3,2)+YCoeff0(1,3,2,2).*YY(1:paths).^Fp(1,3,2,2)+ ...
% YCoeff0(2,2,2,2).*YY(1:paths).^Fp(2,2,2,2)+ YCoeff0(3,1,2,2).*YY(1:paths).^Fp(3,1,2,2)+ ...
% YCoeff0(1,4,1,2).*YY(1:paths).^Fp(1,4,1,2)+YCoeff0(2,3,1,2).*YY(1:paths).^Fp(2,3,1,2)+ ...
% YCoeff0(3,2,1,2).*YY(1:paths).^Fp(3,2,1,2)+YCoeff0(4,1,1,2).*YY(1:paths).^Fp(4,1,1,2)).*dtM^3.5) .*HermiteP1(2,1:paths)
end
```

.

.

.

```
function [dwdZ,d2wdZ2,d3wdZ3] = First3Derivatives4thOrderEqSpaced(wnStart,Nn,dNn,w,Z)
% [wS] = InterpolateOrderN6(4,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5));
% [wE] = InterpolateOrderN6(4,Z(Nn)+dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));
% [wS] = InterpolateOrderN8(7,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
% [wE] = InterpolateOrderN8(7,Z(Nn)+dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),Z(Nn-6),Z(Nn-7),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5),w(Nn-6),w(Nn-7));
[wS] = InterpolateOrderN6(6,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5));
[wE] = InterpolateOrderN6(6,Z(Nn)+dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));
[wS2] = InterpolateOrderN6(6,Z(wnStart)-2*dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5));
[wE2] = InterpolateOrderN6(6,Z(Nn)+2*dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));
% Zi=[Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),Z(wnStart+8),Z(wnStart+9)];
% wi=[w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7),w(wnStart+8),w(wnStart+9)];
% % Zi=[Z(Nn-6),Z(Nn-5),Z(Nn-4),Z(Nn-3),Z(Nn-2),Z(Nn-1),Z(Nn)];
% % wi=[w(Nn-6),w(Nn-5),w(Nn-4),w(Nn-3),w(Nn-2),w(Nn-1),w(Nn)];
%
% Zi=[Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6)];
% wi=[w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6)];
% Zq=Z(wnStart)-dNn;
% wS=spline(Zi,wi,Zq);
%
% Zq=Z(wnStart)-2*dNn;
% wS2=spline(Zi,wi,Zq);
%
%
%
% % Zi=[Z(Nn-9),Z(Nn-8),Z(Nn-7),Z(Nn-6),Z(Nn-5),Z(Nn-4),Z(Nn-3),Z(Nn-2),Z(Nn-1),Z(Nn)];
% % wi=[w(Nn-9),w(Nn-8),w(Nn-7),w(Nn-6),w(Nn-5),w(Nn-4),w(Nn-3),w(Nn-2),w(Nn-1),w(Nn)];
% Zi=[Z(Nn-6),Z(Nn-5),Z(Nn-4),Z(Nn-3),Z(Nn-2),Z(Nn-1),Z(Nn)];
% wi=[w(Nn-6),w(Nn-5),w(Nn-4),w(Nn-3),w(Nn-2),w(Nn-1),w(Nn)];
% Zq=Z(Nn)+dNn;
% wE=spline(Zi,wi,Zq);
%
%
% Zq=Z(Nn)+2*dNn;
% wE2=spline(Zi,wi,Zq);
ZS=Z(wnStart)-dNn;
ZE=Z(Nn)+dNn;
% d2wdZ2(wnStart)=(Z(wnStart)*wS-2*w(wnStart)*Z(wnStart)+w(wnStart+1).*Z(wnStart))/(dNn.^2);
% d2wdZ2(wnStart+1:Nn-1)=(w(wnStart+1:Nn-1).*Z(wnStart+1:Nn-1)-2*w(wnStart+1:Nn-1).*Z(wnStart+1:Nn-1)+1*w(wnStart+2:Nn).*Z(wnStart+1:Nn-1))/(dNn.^2);
% d2wdZ2(Nn)=(wE.*Z(Nn)-2*w(Nn).*Z(Nn)+w(Nn-1).*Z(Nn))/(dNn.^2);
dwdZ(wnStart+2:Nn-2)=(1*w(wnStart:Nn-4)-8*w(wnStart+1:Nn-3)+0*w(wnStart+2:Nn-2)+8*w(wnStart+3:Nn-1)-1*w(wnStart+4:Nn))/(12*dNn);
dwdZ(wnStart+1)=(-3*w(wnStart)-10*w(wnStart+1)+18*w(wnStart+2)-6*w(wnStart+3)+1*w(wnStart+4))/(12*dNn);
dwdZ(wnStart)=(-25/12*w(wnStart)+4*w(wnStart+1)-3*w(wnStart+2)+4/3*w(wnStart+3)-1/4*w(wnStart+4))/dNn;%/(12*dNn);
dwdZ(Nn-1)=(-1*w(Nn-4)+6*w(Nn-3)-18*w(Nn-2)+10*w(Nn-1)+3*w(Nn))/(12*dNn);
dwdZ(Nn)=(1/4*w(Nn-4)-4/3*w(Nn-3)+3*w(Nn-2)-4*w(Nn-1)+25/12*w(Nn))/dNn;%/(12*dNn);
%−25/12 4 −3 4/3 −1/4
%(-3*f[i-1]-10*f[i+0]+18*f[i+1]-6*f[i+2]+1*f[i+3])/(12*1.0*h**1)
% (-1*f[i-3]+6*f[i-2]-18*f[i-1]+10*f[i+0]+3*f[i+1])/(12*1.0*h**1)
d2wdZ2(wnStart+2:Nn-2)=(-1*w(wnStart:Nn-4)+16*w(wnStart+1:Nn-3)-30*w(wnStart+2:Nn-2)+16*w(wnStart+3:Nn-1)-1*w(wnStart+4:Nn))/(12*dNn^2);
d2wdZ2(wnStart+1)=(10*w(wnStart)-15*w(wnStart+1)-4*w(wnStart+2)+14*w(wnStart+3)-6*w(wnStart+4)+1*w(wnStart+5))/(12*dNn^2);
d2wdZ2(wnStart)=(15/4*w(wnStart)-77/6*w(wnStart+1)+107/6*w(wnStart+2)-13*w(wnStart+3)+61/12*w(wnStart+4)-5/6*w(wnStart+5))/(dNn^2);
d2wdZ2(Nn-1)=(10*w(Nn)-15*w(Nn-1)-4*w(Nn-2)+14*w(Nn-3)-6*w(Nn-4)+1*w(Nn-5))/(12*dNn^2);
%d2wdZ2(Nn)=(-1*w(Nn-2)+16*w(Nn-1)-30*w(Nn)+16*wE-1*wE2)/(12*dNn^2);
d2wdZ2(Nn)=(15/4*w(Nn)-77/6*w(Nn-1)+107/6*w(Nn-2)-13*w(Nn-3)+61/12*w(Nn-4)-5/6*w(Nn-5))/(dNn^2);
% (1*f[i-4]-6*f[i-3]+14*f[i-2]-4*f[i-1]-15*f[i+0]+10*f[i+1])/(12*1.0*h**2)
% (10*f[i-1]-15*f[i+0]-4*f[i+1]+14*f[i+2]-6*f[i+3]+1*f[i+4])/(12*1.0*h**2)
% 15/4 −77/6 107/6 −13 61/12 −5/6
%
%
% d2wdZ2(wnStart)=(wS-2*w(wnStart)+w(wnStart+1))/(dNn.^2);
% d2wdZ2(wnStart+1:Nn-1)=(w(wnStart:Nn-2)-2*w(wnStart+1:Nn-1)+1*w(wnStart+2:Nn))/(dNn.^2);
% d2wdZ2(Nn)=(wE-2*w(Nn)+w(Nn-1))/(dNn.^2);
%
%
% dwdZ(wnStart)=(-1*wS+1*w(wnStart+1))/(2*dNn);
% dwdZ(wnStart+1:Nn-1)=(-1*w(wnStart:Nn-2)+1*w(wnStart+2:Nn))/(2*dNn);
% dwdZ(Nn)=(wE-w(Nn-1))/(2*dNn);
%
d3wdZ3(wnStart+2:Nn-2)=(-1*w(wnStart:Nn-4)+2*w(wnStart+1:Nn-3)+0*w(wnStart+2:Nn-2)-2*w(wnStart+3:Nn-1)+w(wnStart+4:Nn))/(2*dNn^3);
d3wdZ3(wnStart+1)=(-1*wS+2*w(wnStart)+0*w(wnStart+1)-2*w(wnStart+2)+w(wnStart+3))/(2*dNn^3);
d3wdZ3(wnStart)=(-1*wS2+2*wS+0*w(wnStart)-2*w(wnStart+1)+w(wnStart+2))/(2*dNn^3);
d3wdZ3(Nn-1)=(-1*w(Nn-3)+2*w(Nn-2)+0*w(Nn-1)-2*w(Nn)+wE)/(2*dNn^3);
d3wdZ3(Nn)=(-1*w(Nn-2)+2*w(Nn-1)+0*w(Nn)-2*wE+wE2)/(2*dNn^3);
end
```