Sorry Friends, the algorithm I posted yesterday to calculate variances had an extra unnecessary nested loop. I have removed the redundant loop and now the algorithm is way faster and you can far easily use it for calculation of covariances on longer time horizons.

Here is the new modified efficient program in place of program posted yesterday.

```
function [] = SDETransProb08WmtGrid200b2()
%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
%Besse1 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.
%When you run this program, please make sure that you look at graph zoomed
%at main body of the graph. Sometimes all the SDE is properly simulated but
%the tail explodes for some reason so you would know. These are minor
%issues that will be fixed in next versions of the program.
%I have put the terminal date at a very small value. Since calculations are
%done for all possible pairs of covariances on a very large number of time
%steps, increasing time to a large value is going to take quite a while.
%Start with a small date. You might want to modify this algorithm to
%calculate covariance pairs certain time spacing larger than minimum time spacing.
T=1/4;
dt=.125/16/2/2; % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=T*128*2*2; % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
T=Tt*dt;
OrderA=4; %
dtM=.125/2;%Monte carlo time interval size dtM.
TtM=T*2*8;%Monte carlo number of simulation intervals.
dNn=.05; % Normal density subdivisions width. would change with number of subdivisions
Nn=200; % No of normal density subdivisions
NnMidl=100;%One half density Subdivision left from mid of normal density(low)
NnMidh=101;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;
Nn2=280;
dNn2=.05;
Nn2Midl=140;
Nn2Midh=141;
Nn2Mid=4.0;
x0=1.00; % starting value of SDE
beta1=0.0;
beta2=1.00; % Second drift term power.
gamma=.95;%50; % volatility power.
kappa=0.0;%.950; %mean reversion parameter.
theta=1.0;%mean reversion target
sigma0=.850;%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
w(1:Nn)=x0^(1-gamma)/(1-gamma);
%Z(1:Nn)=(((1:Nn)-5.5)*dNn-NnMid);
%Z(1:Nn)=(((1:Nn)-20.5)*dNn-NnMid);
Z(1:Nn)=(((1:Nn)-20.5)*dNn-NnMid);
%Z(1)=-5.118453;
%Z(Nn)=5.118453;
Z2(1:Nn2)=(((1:Nn2)-60.5)*dNn2-Nn2Mid);
Z
Z2
str=input('Look at Zs')
%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);
%Above calculate probability mass in each probability subdivision.
ZProb(1)=normcdf(-4.95);
ZProb(Nn)=1-normcdf(4.95);
ZProb(2)=normcdf(-4.9)-normcdf(-4.95);
ZProb(Nn-1)=normcdf(4.95)-normcdf(4.9);
ZProb(3:Nn-2)=normcdf(.5*Z(3:Nn-2)+.5*Z(4:Nn-1),0,1)-normcdf(.5*Z(3:Nn-2)+.5*Z(2:Nn-3),0,1);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Current program has 50 grid points. I have divided every grid subdivision
%into three further subdivisions with probabilities P0,Pb,Pa.
wnStart=1;%
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
%The above yy(wnStart:Nn)=x0;
Pprev(1:Nn)=ZProb(1:Nn);%probability mass at starting node.
Pnew(1:Nn2)=0.0;
EZwnew(1:Nn2)=0.0;
EZwprev(1:Nn)=0.0;
Ewnew(1:Nn2)=0;
Ewprev(1:Nn)=0;
Ew1new(1:Nn2)=0;
Ew1prev(1:Nn)=0;
CovarZt(1:Tt,1:Tt)=0.0;
Covarwt(1:Tt,1:Tt)=0.0;
%ECovZnew0(1:Tt,1:Tt,1:Nn2)=0.0;
ECovZprev0(1:Tt,1:Nn)=0.0;
%ECovwnew0(1:Tt,1:Tt,1:Nn2)=0.0;
ECovwprev0(1:Tt,1:Nn)=0.0;
VarZ(1:Tt)=0;
Varw(1:Tt)=0;
ValZprev(1:Nn)=0;
Valwprev(1:Nn)=0;
ValZnew(1:Nn2)=0;
Valwnew(1:Nn2)=0;
ValZnow(1:Nn)=0;
Valwnow(1:Nn)=0;
Zwt0(1:Nn)=0;
wt0(1:Nn)=0;
tic
load('C:\\Users\\Ahsan\\Documents\\MATLAB\\TrProbs0O4b.mat','pt0','Zt0');
load('C:\\Users\\Ahsan\\Documents\\MATLAB\\TrProbsO4b.mat','ptt','Ztt','ptA','ptB');
%load('C:\\Users\\Ahsan\\Documents\\MATLAB\\FPEwI.mat','wI');
%load('C:\\Users\\Ahsan\\Documents\\MATLAB\\FPEdwZI.mat','dwZI');
for tt=1:Tt
t2=tt*dt;
t1=(tt-1)*dt;
if(tt==1)
[wMu0dt,c1,c22] = CalculateDriftAndVolA8Trans(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
Pnew(1:Nn2)=pt0(1:Nn2);
EZwnew(1:Nn2)=pt0(1:Nn2).*Zt0(1:Nn2);
Ewnew(1:Nn2)=(x0^(1-gamma)/(1-gamma)+wMu0dt(1)+c1(1).*Zt0(1:Nn2)+c22(1).*(Zt0(1:Nn2).^2-1)).*pt0(1:Nn2);
Ew1new(1:Nn2)=(x0^(1-gamma)/(1-gamma)+wMu0dt(1)+c1(1).*Zt0(1:Nn2)+c22(1).*(Zt0(1:Nn2).^2-1)).*pt0(1:Nn2);
end
if(tt>1)
%Convert the expected values in each cell to standard Values associated with the grid cell after removing the
%integrated probability in the cell.
w(1:Nn)=Ewprev(1:Nn)./Pprev(1:Nn);
w1(1:Nn)=Ew1prev(1:Nn)./Pprev(1:Nn);
[wMu0dt,wMu1dt,dwMu0dtdw,d2wMu0dtdw2,c1] = CalculateDriftAndVolA404(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
[w1Mu0dt,w1Mu1dt,dw1Mu0dtdw1,d2w1Mu0dtdw2,cw1] = CalculateDriftAndVolA404(w1,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
pt(1:Nn,1:Nn2)=ptt(tt,1:Nn,1:Nn2);
Zt(1:Nn,1:Nn2)=Ztt(tt,1:Nn,1:Nn2);
%#Girsanov
%For solution of drifted SDEs with Girsanov,uncomment the next three lines.
% wmStart=1;
% Mm=Nn2;
% [GirsanovExp] = CalculateGirsanovExponential(wnStart,Nn,wmStart,Mm,w,Zt,dt,mu1,mu2,beta1,beta2,gamma,sigma0);
for mm=1:Nn
%#Girsanov
%Uncomment the next line to see that expected value of Girsanov exponential
%should be unity.
% sum(GirsanovExp(mm,1:Nn2).*pt(mm,1:Nn2))
A0=ptA(mm,tt);
B0=ptB(mm,tt);
%A0 and B0 are starting and ending boundaries for non-zero
%elements in larger target array. These are boundaries for
%transition probabilities originating from cell mm on first grid
% at time tt. Out of the 280 target grid
%cell, only a few ten are non-zero. These non-zero elements
%boundaries(A0 and B0) were calculated in a different module and saved on
%disk just when we calculated the transition probabilities.
%Later we retrieved them from disk to use in this program.
Pnew(A0:B0)=Pnew(A0:B0)+Pprev(mm).*pt(mm,A0:B0);
EZwnew(A0:B0)=EZwnew(A0:B0)+EZwprev(mm).*pt(mm,A0:B0)+ ...
Zt(mm,A0:B0).*pt(mm,A0:B0).*Pprev(mm);
%Below is the transition probabilities simulation
%with improvised term. These improvised terms make an improvement
%over the base model which is given in next paragraph of code.
Ewnew(A0:B0)=Ewnew(A0:B0)+Ewprev(mm).*pt(mm,A0:B0)+ ...
(wMu0dt(mm)+c1(mm).*Zt(mm,A0:B0)).* ...
pt(mm,A0:B0).*Pprev(mm);
%#Girsanov
%For Girsanov, uncomment the next block and comment the previous block.
% Ewnew(A0:B0)=Ewnew(A0:B0)+Ewprev(mm).*pt(mm,A0:B0)+ ...
% GirsanovExp(mm,A0:B0).*sigma0.*sqrt(dt).*Zt(mm,A0:B0).* ...
% pt(mm,A0:B0).*Pprev(mm);
Ewnew(Ewnew<0)=0;
Ewnew(isnan(Ewnew)==1)=0;
%Below is the base transition probabilities simulation
%with same terms as in monte carlo simulation. cw1 term is analogue
%of c1 which is coefficient of Zt in our bessel form SDE.
Ew1new(A0:B0)=Ew1new(A0:B0)+Ew1prev(mm).*pt(mm,A0:B0)+ ...
(w1Mu0dt(mm)+cw1(mm).*Zt(mm,A0:B0)).* ...
pt(mm,A0:B0).*Pprev(mm);
Ew1new(Ew1new<0)=0;
Ew1new(isnan(Ew1new)==1)=0;
end
wa(1:Nn-1)=w(1:Nn-1);
wa(Nn)=w(Nn);
wb(1:Nn-1)=w(2:Nn);
wb(Nn)=w(Nn-1)+w(Nn-1)-w(Nn-2);
w(wa(:)>wb(:))=wb(wa(:)>wb(:));%Be careful;might not universally hold;
%if(w(Nn)>(2*w(Nn-1)-w(Nn-2)))
% w(Nn)=(2*w(Nn-1)-w(Nn-2));
%end
end
Pprev(1)=sum(Pnew(1:41));
Pprev(2:Nn-1)=Pnew(42:40+Nn-1);
Pprev(Nn)=sum(Pnew(40+Nn:80+Nn));%All the probability mass in last 41
%subdivisions on target grid is assigned to end boundary point on
%the grid that will become smaller first grid for next simulation step.
EZwprev(1)=sum(EZwnew(1:41));
EZwprev(2:Nn-1)=EZwnew(42:40+Nn-1);
EZwprev(Nn)=sum(EZwnew(40+Nn:80+Nn));
Zwt(tt,1:Nn)=EZwprev(1:Nn)./Pprev(1:Nn);
MuZ(tt)=sum(EZwprev(1:Nn));
Ewprev(1)=sum(Ewnew(1:41));
Ewprev(2:Nn-1)=Ewnew(42:40+Nn-1);
Ewprev(Nn)=sum(Ewnew(40+Nn:80+Nn));
wt(tt,1:Nn)=Ewprev(1:Nn)./Pprev(1:Nn);
Muw(tt)=sum(Ewprev(1:Nn));
Ew1prev(1)=sum(Ew1new(1:41));
Ew1prev(2:Nn-1)=Ew1new(42:40+Nn-1);
Ew1prev(Nn)=sum(Ew1new(40+Nn:80+Nn));
Pnew(1:Nn2)=0;
EZwnew(1:Nn2)=0.0;
Ewnew(1:Nn2)=0.0;
Ew1new(1:Nn2)=0.0;
%#Covariances Calculation.
%Below, we calculate covariances between the each date with all of the
%later following dates
if(tt==1)
VarZ(tt)=sum((Zwt(tt,1:Nn)-MuZ(tt)).^2.*ZProb(1:Nn));
Varw(tt)=sum((wt(tt,1:Nn)-Muw(tt)).^2.*ZProb(1:Nn));
end
if(tt>1)
VarZ(tt)=sum((Zwt(tt,1:Nn)-MuZ(tt)).^2.*ZProb(1:Nn));
Varw(tt)=sum((wt(tt,1:Nn)-Muw(tt)).^2.*ZProb(1:Nn));
%Below loop with variable tt1 indicates the starting time for various
%covariances
for tt1=1:tt-1
tt1
tt
if(tt1<tt-1)
%If the starting time for covariance calculations is smaller than one date
%ago, it means that we do not have to initialize the algorithm and we are
%in the middle of the calculations.
ValZprev(1:Nn)=ECovZprev0(tt1,1:Nn);
Valwprev(1:Nn)=ECovwprev0(tt1,1:Nn);
end
if(tt1==tt-1)
%If the starting time for covariance calculations is equal to one date
%ago, it means that we haveto initialize the algorithm and we are
%starting the calculations. Here we assign the starting values for covariance
%calculation.
Zwt0(1:Nn)=Zwt(tt-1,1:Nn);
ValZprev(1:Nn)=(Zwt0(1:Nn)-MuZ(tt-1)).*ZProb(1:Nn);
wt0(1:Nn)=wt(tt-1,1:Nn);
Valwprev(1:Nn)=(wt0(1:Nn)-Muw(tt-1)).*ZProb(1:Nn);
end
ValZnew(1:Nn2)=0;
Valwnew(1:Nn2)=0;
%Below we calculate the running string for the transition probability increment
%between first covariance date and last covariance date.
for mm=1:Nn
A0=ptA(mm,tt);
B0=ptB(mm,tt);
pt(mm,1:Nn2)=ptt(tt,mm,1:Nn2);
ValZnew(A0:B0)=ValZnew(A0:B0)+ ...
ValZprev(mm).*pt(mm,A0:B0);
Valwnew(A0:B0)=Valwnew(A0:B0)+ ...
Valwprev(mm).*pt(mm,A0:B0);
end
%ECovZnew0(tt1,tt2,1:Nn2)=ValZnew(1:Nn2);
%ECovwnew0(tt1,tt2,1:Nn2)=Valwnew(1:Nn2);
ECovZprev0(tt1,1)=sum(ValZnew(1:41));
ECovZprev0(tt1,2:Nn-1)=ValZnew(42:40+Nn-1);
ECovZprev0(tt1,Nn)=sum(ValZnew(40+Nn:80+Nn));
ECovwprev0(tt1,1)=sum(Valwnew(1:41));
ECovwprev0(tt1,2:Nn-1)=Valwnew(42:40+Nn-1);
ECovwprev0(tt1,Nn)=sum(Valwnew(40+Nn:80+Nn));
%Below when second covariance pair date is equal to present date, we do
%final(ending) calculations for covariance calculation with all first pairs from
%all of the earlier dates.
ValZnow(1:Nn)=Zwt(tt,1:Nn)-MuZ(tt);
Valwnow(1:Nn)=wt(tt,1:Nn)-Muw(tt);
ValZprev(1:Nn)=ECovZprev0(tt1,1:Nn);
Valwprev(1:Nn)=ECovwprev0(tt1,1:Nn);
CovarZt(tt1,tt)=sum(ValZprev(1:Nn).*ValZnow(1:Nn));
Covarwt(tt1,tt)=sum(Valwprev(1:Nn).*Valwnow(1:Nn));
end
end
% %#Covariances Calculation.
% %Below, we calculate covariances between the each date with all of the
% %later following dates
%
% if(tt==1)
% VarZ(tt)=sum((Zwt(tt,1:Nn)-MuZ(tt)).^2.*ZProb(1:Nn));
% Varw(tt)=sum((wt(tt,1:Nn)-Muw(tt)).^2.*ZProb(1:Nn));
% end
% if(tt>1)
% VarZ(tt)=sum((Zwt(tt,1:Nn)-MuZ(tt)).^2.*ZProb(1:Nn));
% Varw(tt)=sum((wt(tt,1:Nn)-Muw(tt)).^2.*ZProb(1:Nn));
%
% %Below loop with variable tt1 indicates the starting time for various
% %covariances
% %Below loop variable tt2 indicates the algorithm calculations on dates
% %following the strting time.
% for tt1=1:tt-1
% for tt2=tt:Tt
% tt1
% tt2
% if(tt1<tt-1)
% %If the starting time for covariance calculations is smaller than one date
% %ago, it means that we do not have to initialize the algorithm and we are
% %in the middle of the calculations.
% ValZprev(1:Nn)=ECovZprev0(tt1,tt2,1:Nn);
% Valwprev(1:Nn)=ECovwprev0(tt1,tt2,1:Nn);
% end
% if(tt1==tt-1)
% %If the starting time for covariance calculations is equal to one date
% %ago, it means that we haveto initialize the algorithm and we are
% %starting the calculations. Here we assign the starting values for covariance
% %calculation.
% Zwt0(1:Nn)=Zwt(tt-1,1:Nn);
% ValZprev(1:Nn)=(Zwt0(1:Nn)-MuZ(tt-1)).*ZProb(1:Nn);
% wt0(1:Nn)=wt(tt-1,1:Nn);
% Valwprev(1:Nn)=(wt0(1:Nn)-Muw(tt-1)).*ZProb(1:Nn);
% end
% ValZnew(1:Nn2)=0;
% Valwnew(1:Nn2)=0;
% %Below we calculate the running string for the transition probability increment
% %between first covariance date and last covariance date.
% for mm=1:Nn
% A0=ptA(mm,tt);
% B0=ptB(mm,tt);
% pt(mm,1:Nn2)=ptt(tt,mm,1:Nn2);
% ValZnew(A0:B0)=ValZnew(A0:B0)+ ...
% ValZprev(mm).*pt(mm,A0:B0);
% Valwnew(A0:B0)=Valwnew(A0:B0)+ ...
% Valwprev(mm).*pt(mm,A0:B0);
% end
% %ECovZnew0(tt1,tt2,1:Nn2)=ValZnew(1:Nn2);
% %ECovwnew0(tt1,tt2,1:Nn2)=Valwnew(1:Nn2);
% ECovZprev0(tt1,tt2,1)=sum(ValZnew(1:41));
% ECovZprev0(tt1,tt2,2:Nn-1)=ValZnew(42:40+Nn-1);
% ECovZprev0(tt1,tt2,Nn)=sum(ValZnew(40+Nn:80+Nn));
%
% ECovwprev0(tt1,tt2,1)=sum(Valwnew(1:41));
% ECovwprev0(tt1,tt2,2:Nn-1)=Valwnew(42:40+Nn-1);
% ECovwprev0(tt1,tt2,Nn)=sum(Valwnew(40+Nn:80+Nn));
% end
% %Below when second covariance pair date is equal to present date, we do
% %final(ending) calculations for covariance calculation with all first pairs from
% %all of the earlier dates.
% ValZnow(1:Nn)=Zwt(tt,1:Nn)-MuZ(tt);
% Valwnow(1:Nn)=wt(tt,1:Nn)-Muw(tt);
% ValZprev(1:Nn)=ECovZprev0(tt1,tt,1:Nn);
% Valwprev(1:Nn)=ECovwprev0(tt1,tt,1:Nn);
% CovarZt(tt1,tt)=sum(ValZprev(1:Nn).*ValZnow(1:Nn));
% Covarwt(tt1,tt)=sum(Valwprev(1:Nn).*Valwnow(1:Nn));
%
% end
% end
end
%At the end, we want to see the results for our covariance calulations
for tt=1:Tt-1
%Below, we see the variance of standard increments normal on a particular date
VarZ(tt)
%Now, we want to see the covariance of all the following dates with present date
%and for normal noise this should be equal to the variance at starting
%date.
CovarZt(tt,tt+1:Tt)
tt
str=input('Look at CovarZt');
%Below, we see the variance of a general SDE(but this current program works well
%with driftless noises only) on a particular date
Varw(tt)
%Now, we want to see the covariance of all the following dates with present date
%and for CEV type driftless SDEs, this should be equal to the variance at starting
%date.
Covarwt(tt,tt+1:Tt)
tt
str=input('Look at Covarwt');
end
plot((1:Nn),ZProb(1:Nn),'g',(1:Nn),Pprev(1:Nn),'r');
%plot((1:Nn),ZProb(1:Nn),'g',(1:Nn),Pprev(1:Nn),'r');
%str=input('Look at distributions');
Pnew=Pprev;
%Convert the expected values in each cell to standard Values associated with the grid cell after removing the
%integrated probability in the cell.
w(wnStart:Nn)=Ewprev(wnStart:Nn)./Pprev(wnStart:Nn);
w1(wnStart:Nn)=Ew1prev(wnStart:Nn)./Pprev(wnStart:Nn);
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
yy1(wnStart:Nn)=((1-gamma)*w1(wnStart:Nn)).^(1/(1-gamma));
%below D's (the names of variables starting with D) are
%change of probability derivatives.
y_w(1:Nn)=0;
y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
y_w1(1:Nn)=0;
y_w1(wnStart:Nn) = ((1-gamma)*w1(wnStart:Nn)).^(1/(1-gamma));
Dfy_w(wnStart:Nn)=0;
Dfy_w1(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
Dfy_w(nn) = (y_w(nn + 1) - y_w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
Dfy_w1(nn) = (y_w1(nn + 1) - y_w1(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
py_w(1:Nn)=0;
py_w1(1:Nn)=0;
for nn = wnStart:Nn-1
py_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
py_w1(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w1(nn));
end
toc
ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-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
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.
MCVar=sum((YY(:)-MCMean).^2)/paths
disp('Original process average from our simulation');
ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates
ItoHermiteVar=sum((y_w(wnStart+1:Nn-1)-ItoHermiteMean).^2.*ZProb(wnStart+1:Nn-1))
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(y_w1(wnStart+1:Nn-1),py_w1(wnStart+1:Nn-1),'b',y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
plot(y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',y_w1(wnStart+1:Nn-1),py_w1(wnStart+1:Nn-1),'b',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
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({'Improvised terms Model Density','Monte Carlo Density'},'Location','northeast')
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
end
```

.

Here is the last output of the program.

true Mean only applicble to standard SV mean reverting type models otherwise disregard