SERVING THE QUANTITATIVE FINANCE COMMUNITY

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

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

I want friends to protest against my increasing persecution. For past 3-4 days, mind control agencies have been trying to viciously drug me. I have seen several incidents of outside people coming into bakeries after I arrive there and taking position with bakery staff. They even drugged food once inside my locked car after I had gone inside a store and had parked my car outside the store but I got to know about it and threw the food. Mind control agencies have totally gone rabid. Another grave issue I noticed today was that my LinkedIn connections all of a sudden deceased from 2988 to 2936. Just a day or two ago, I noticed that my connections were 2988 and I was thinking that I would get to 3000 connections in a few weeks. They have removed a large number of good people from my Linkedin connections. I have not myself removed a single connection. I truly value good and respected people in my Linkedin especially the ones in Academia. I will request those people who notice that they have been removed to please connect with me again. I deeply apologize to people they have removed from my linkedin.
For the context of my persecution, please see
My Letter to HRW About My Human Rights Abuse : https://forum.wilmott.com/viewtopic.php?f=15&t=102298
Muslims, mind control, Muslim-Jewish relationship and American defensehttps://forum.wilmott.com/viewtopic.php?f=15&t=101453
My Open Letter to Prime Minister of Pakistan, Imran Khan, to Intervene In Order to Protect the Interests of Pakistan  https://forum.wilmott.com/viewtopic.php?f=15&t=94796&start=1125#p858119

I was extremely lucky to have noticed that mind control agencies have removed more than fifty people from my linkedin connections within a span of two days and mentioned it here. Machinatory people would have been telling those good people who have struggled to protect me that I have removed them from my linkedin connections.Please do not believe in a single word of these machinatory people and their schemes. It was machinations all along for past twenty three years. And I will make an apology to  all the good people they have removed from my connections and who had made efforts to protect me.
I want to take this opportunity to tell the good people who made efforts to protect me that I am indebted to them all my life for what they did for me and this is something so meaningful for me to think that I could be free in a few months and I can never repay the debt good people have on me.

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

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

First of all I am very sorry to all the friends for such a long delay in the new completely worked out version of the program. I got an antipsychotic injection in the first week of this month and I really could not work well and it was difficult to concentrate. Anyway, I am posting a new interim version program. It just has a few minor changes and also has the bad integral in previous version done right. I am posting it since I believe that the main algorithm will see only a few minor changes and will remain the same as in this program. And the program still has problems and I believe that reasons behind the problems are as follows:
1. We need to take a smaller step size for transition probabilities calculations. Currently I am using Freq=8 which means that one transition probabilities step size is equal to eight steps in the fokker-planck equation. Freq=4.0 and Freq=2.0 give far better results(if you can get it to work) but on these smaller steps, transition probabilities become unstable and create problems. So major issue is to have very good transition probabilities on smaller step levels. I plan to use my mean-correction algorithm that changes the transition probabilities slightly while making sure they sum exactly to one and remain stable. Later I also want to introduce mean correction for transition probabilities of SDEs where mean is not analytically known.

2. I also want to try whether using a Z-grid of 100 grid points improves accuracy as opposed to grid of 50 grid points that we are currently using.

3. For high kappa and for fast moving densities, I will try if adding higher order iterated integrals would add to accuracy.

There are other such minor issues that affect the performance of the algorithm and once solved, we would be able to get a good result even with the same main algorithm we have used for path integrals.

Here is the new program with a few minor changes. I have also posted the dependency function that is used for calculation of variance integrals.
.
.
function [] = FPERevisitedTransProb08()

%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%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*4/4*2*2;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
%Tt=8*60;
T=Tt*dt;
OrderA=4;  %
OrderM=4;  %
%dtM=.125/8;%Monte carlo time interval size dtM.
%TtM=8*8;%Monte carlo number of simulation intervals.
dtM=dt*4*2;
TtM=Tt/4/2;

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;

%theta=mu1/(-mu2);
%kappa=-mu2;

x0=.250;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.95;%50;   % volatility power.
kappa=1.0;%.950;   %mean reversion parameter.
theta=.10;%mean reversion target
sigma0=.70;%Volatility value
omega1=1.0; %power on yy(t) for fy and fyPI.
%This is yy(t)^omega1 is the function that is summed in path integral
%I have tested powers between omega1=.5 and omega1=2;
theta1=1;%This is the weight on arithmetic sum in the path integral.
%keep it equal to one otherwise results will be erroneous.
%This will be activated in future and made time dependent.

%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
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);
%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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ss=0;  %ss is index for transition probability calculation time steps.
Freq=8.0;%transition probabilities are calculated Freq time intervals apart.
wnStart=1;%
d2wdZ2(1:Nn)=0;
d2wdZ2A(1:Nn)=0;
dwdZ(1:Nn)=0;
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:Nn)=0.0;
Efprev(1:Nn)=theta1*yy(1:Nn).^omega1.*ZProb(1:Nn);%value of function at
%starting nodes.
Efnew(1:Nn)=0;
Eyyprev(1:Nn)=yy(1:Nn).*ZProb(1:Nn);%Value of SDE variable in original
%coordinates at starting nodes.
Eyynew(1:Nn)=0;
Efdtprev(1:Nn)=0;%Value of path-integral at starting nodes is zero.
Efdtnew(1:Nn)=0;

tic

for tt=1:Tt
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
[wMu0dt,dwMu0dtdw,c1] = CalculateDriftAndVolA4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn) ;% ...
%dw(wnStart:Nn)=sigma0.*sqrt(dt).*Z(wnStart:Nn) ;
dw2(wnStart:Nn)=dw(wnStart:Nn).^2;
w(isnan(w)==1)=0;
wMu0dt(isnan(wMu0dt)==1)=0;
%wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%     [dwdZ,d2wdZ2A,d3wdZ3A] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z);
%     d2wdZ2(wnStart:Nn)=0.0;
%     d3wdZ3(wnStart:Nn)=0.0;
%     if (tt>36)
%         d3wdZ3(wnStart:Nn) = (d3wdZ3A(wnStart:Nn));
%         d2wdZ2(wnStart:Nn)=d2wdZ2A(wnStart:Nn);
%         d2wdZ2(1:15)=d2wdZ2(1:15).*exp(-12*kappa^2*(1-gamma).*dNn*(15-(1:15)));
%         d3wdZ3(1:15)=d3wdZ3(1:15).*exp(-12*kappa^2*(1-gamma).*dNn*(15-(1:15)));
%
%     end

%tt
%plot((wnStart:Nn),d2wdZ2A(wnStart:Nn),'g',(wnStart:Nn),d2wdZ2(wnStart:Nn),'r');
%str=input('Look at d2wdZ2A(wnStart:Nn)');

%If the calculation of the program blows, it means that above calculated
%derivative is unstable or bad. Plese replace with your favourite method of
%calculating this second derivative.

dZdw(wnStart:Nn)=1.0./dwdZ(wnStart:Nn);
if(tt==1)
wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));

w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+ ...
sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)));
end
if(tt>1)

C22(wnStart:Nn)=.0125/pi*25.0;
C33(wnStart:Nn)=.0125/pi*2;
C44(wnStart:Nn)=.0125/pi*25.0;

C22(wnStart:Nn)=.0125/pi*25.0*1.5;
C33(wnStart:Nn)=.0125/pi*2*1.0;
C44(wnStart:Nn)=.0125/pi*25.0*1.5;

TermA0(wnStart:Nn)=-C22(wnStart:Nn).*3.*Z(wnStart:Nn).^1.*dZdw(wnStart:Nn).*d2wdZ2(wnStart:Nn)*sigma0^2.*dt + ...
+C44(wnStart:Nn).*(-d3wdZ3(wnStart:Nn).*dZdw(wnStart:Nn)+3*d2wdZ2(wnStart:Nn).^2.*dZdw(wnStart:Nn).^2)*sigma0^2.*dt+ ...
-C33(wnStart:Nn).*dwMu0dtdw(wnStart:Nn).*dwdZ(wnStart:Nn).^2;

TermA(wnStart:Nn)= sqrt(abs(TermA0(wnStart:Nn)));
SignTermA(wnStart:Nn)=sign(TermA0(wnStart:Nn));

w(wnStart:Nn)=w(wnStart:Nn)+wMu0dt(wnStart:Nn)+(-(Z(wnStart:Nn).*dwdZ(wnStart:Nn))+d2wdZ2(wnStart:Nn))+ ...
sign((Z(wnStart:Nn).*dwdZ(wnStart:Nn))-d2wdZ2(wnStart:Nn)+dw(wnStart:Nn)+SignTermA(wnStart:Nn).*TermA(wnStart:Nn)).* ...
sqrt(abs(sign((Z(wnStart:Nn).*dwdZ(wnStart:Nn))-d2wdZ2(wnStart:Nn)).* ...
((Z(wnStart:Nn).*dwdZ(wnStart:Nn))-d2wdZ2(wnStart:Nn)).^2+ ...
sign(+dw(wnStart:Nn)).*dw2(wnStart:Nn)+ ...
SignTermA(wnStart:Nn).*TermA(wnStart:Nn).^2));% + ...
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Below are instability correction equations.
%Interpolation of three points in instability region with lagrange
%polynomials.

w(NnMidl) = InterpolateOrderN8(8,Z(NnMidl),Z(NnMidl-4),Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidh+2),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),w(NnMidl-4),w(NnMidl-3),w(NnMidl-2),w(NnMidl-1),w(NnMidh+2),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5));
w(NnMidh) = InterpolateOrderN8(8,Z(NnMidh),Z(NnMidl-4),Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidh+2),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),w(NnMidl-4),w(NnMidl-3),w(NnMidl-2),w(NnMidl-1),w(NnMidh+2),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5));
w(NnMidh+1) = InterpolateOrderN8(8,Z(NnMidh+1),Z(NnMidl-4),Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidh+2),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),w(NnMidl-4),w(NnMidl-3),w(NnMidl-2),w(NnMidl-1),w(NnMidh+2),w(NnMidh+3),w(NnMidh+4),w(NnMidh+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));

w1(1:Nn-1)=w(1:Nn-1);
w1(Nn)=w(Nn);
w2(1:Nn-1)=w(2:Nn);
w2(Nn)=wE;
w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
%Change 3:7/25/2020: I have improved zero correction in above.
w(w<0)=0.0;
for nn=1:Nn
if(w(nn)<=0)
wnStart=nn+1;
end
end

%         %%1st order mean correction. We know the mean in original
%         %%coordinates so I shift the density into original coordinates,
%         %%apply the mean correction and then transform again into Lamperti
%         %%coordinates. Algebra solves two equations in two unknowns.
%         %%Equations are Sum_0^N{(Y_w(wnStart:Nn)-Y0)*W0}=1 and the second
%         %%equation is Sum_0^N{Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0}=u0
%         %%Two unknows are Y0 and W0. u0 is known mean.
%   tt;
u0=theta+(x0-theta)*exp(-kappa*(tt*dt)); %analytic mean of the density
%
%         %If you are not using stochastic volatility, replace above with
%         %true mean otherwise results would become garbage.
%
Y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
%
Yn2Pn=sum(Y_w(wnStart:Nn).*Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
YnPn=sum(Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
Pn=1.0;%%Sum(ZProb(1:Nn))
Y0=(Yn2Pn-u0*YnPn)/(YnPn-u0*Pn);

Y0Pn=Y0*Pn;
W0=1/(YnPn-Y0Pn);
YCorrected(wnStart:Nn)=Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0;
wCorrected(wnStart:Nn)=(YCorrected(wnStart:Nn).^(1-gamma))./(1-gamma);
w(wnStart:Nn)=wCorrected(wnStart:Nn);
%
%I have disabled the mean correction. The above mean correction is only valid for
%standard mean reverting stochastic volatility type SDEs. To enable mean
%correction please uncomment the above block.

[dwdZ,d2wdZ2A,d3wdZ3A] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z);
d2wdZ2(wnStart:Nn)=0.0;
d3wdZ3(wnStart:Nn)=0.0;
if (tt>36)
d3wdZ3(wnStart:Nn) = (d3wdZ3A(wnStart:Nn));
d2wdZ2(wnStart:Nn)=d2wdZ2A(wnStart:Nn);
d2wdZ2(1:15)=d2wdZ2(1:15).*exp(-12*kappa^2*(1-gamma).*dNn*(15-(1:15)));
d3wdZ3(1:15)=d3wdZ3(1:15).*exp(-12*kappa^2*(1-gamma).*dNn*(15-(1:15)));

end

%Start of functions and path integrals calculations
%For transition probabilities calculations, the time grid is Freq intervals apart.
% so time interval ds for transition probability calculation is
% ds=dt*Freq; and t1A is start time for transition probabilities interval
% while t2A is end time of transition probabilities calculation interval.
% so t2A=t1A+ds; W1 is the bessel coordinates value w at t1A and W2 is
% equal to w at t2A. W1=w(t1A) and W2=w(t2A). W1 and W2 are used for
% transition probabilities calculations.

if(tt==1)

W1(1:Nn)=x0^(1-gamma)/(1-gamma);

ds=dt*Freq;
t1A=0;
t2A=ds;
%below calculates integrals for drift and volatility for transition
%probabilities.
[wMu0dt2,dwMu0dt2,dw02]=CalculateDriftAndVol(mu1,mu2,beta1,beta2,gamma,sigma0,W1,1,Nn,ds);
%below calculates integrals for drift and volatility for yy(t).
%yy(t) is the SDE variable in original coordinates as opposed to
%bessel coordinates.
[yyMu0dt,dyyz0,dyy2z0]=CalculateFunctionDriftAndVol(1.0,1.0,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds);
%below calculates integrals for drift and volatility for
%yy(t)^omega1 which is used in Ito change of variable calculations
%on the lattice. This is not needed for path integrals.

[fMu0dt,dfz0,df2z0]=CalculateFunctionDriftAndVol(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds);

%below calculates integrals for drift and volatility for
%yy(s)^omega1 ds. which is used in path integrals calculations
%on the lattice.
%[fdtMu0dt,dfdtz01,dfdt2z02]=CalculateFdtDriftAndVol08(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds,t1A,t2A);

[fdtMu0dt,dfdtz01,dfdt2z02]=CalculateFdtDriftAndVol08(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds,t1A,t2A);

end

if(rem(tt,Freq)==0)

%When this above condition is true it means that we are at end of
%particular transition probability interval and we have end of
%interval W2(s)=w(t) while we will use the value W1 that was
%previously saved Freq dt time intervals earlier or ds=Freq*dt
%time ago.

ss=ss+1; %transition loop variable and is one at the end of first
%transition probability interval.
%Div0=Div0+.5*ds;

W2(1:Nn)=w(1:Nn);
ds=dt*Freq;
dW2(2:Nn-1)=(W2(3:Nn)-W2(1:Nn-2))/2.0;%For calculation of width of
%interval for integrated probability.
dW2(1)=dW2(2);%rough approximation, will change
dW2(Nn)=dW2(Nn-1);%rough approximation, will change
yy1(1:Nn)=((1-gamma)*W1(1:Nn)).^(1/(1-gamma));%Value in origianl
%coordinates at the start of transition probabilities interval.
yy2(1:Nn)=((1-gamma)*W2(1:Nn)).^(1/(1-gamma));%Value in origianl
%coordinates at the end of transition probabilities interval.
fyy1(1:Nn)=theta1.*yy1(1:Nn).^omega1;%function value at the start
%both for Ito calculation of function and path integrals.
fyy2(1:Nn)=theta1.*yy2(1:Nn).^omega1;%function value at the end
%both for Ito calculation of function and path integrals.

for mm=1:Nn
%Below is the gaussian increment between two grid points W1(mm)
% at time t1A and W2(wnStart:Nn) at time t2A, backed
%out from the original stochastic differential equation drift
%and volatility. This gaussian increment is used for
%probability calculation between corresponding grid points and
%also for calculation of stochastic integrals between the
%corresponding grid points. Please note that this gaussian
%increment remains the same in bessel coordinates and the
%original coordinates
Zt(mm,wnStart:Nn)=(W2(wnStart:Nn)-W1(mm)-1*wMu0dt2(mm))/dw02(mm);
%Below is calculation of transition probability between two
%grid points W1(mm)  at time t1A and W2(wnStart:Nn) at time t2A
%using the gaussian increment preeviously calculated.
pt(mm,wnStart:Nn)=exp(-.5* Zt(mm,wnStart:Nn).^2)./sqrt(2*pi*dw02(mm).^2).*dW2(wnStart:Nn);
pt(isnan(pt(mm,:))==1)=0;
pt(isinf(pt(mm,:))==1)=0;
%Below is the counterpart of original SDE that is used for the
%evolution of expcted value of yy(t). It is totally independent
%for its own sake and is not needed for any function or
%path integral calculations.
dyyGenerator(mm,wnStart:Nn)=yyMu0dt(mm)+dyyz0(mm).*Zt(mm,wnStart:Nn)+dyy2z0(mm).*(Zt(mm,wnStart:Nn).^2-1);

%Below is the Ito-generator for fy=theta1*yy(t).^omega1 which
%is used for the evolution of expected value of fy and this is
%totally independent for its own sake and is not needed in any
%other calculations.
dfGenerator(mm,wnStart:Nn)=fMu0dt(mm)+dfz0(mm).*Zt(mm,wnStart:Nn)+df2z0(mm).*(Zt(mm,wnStart:Nn).^2-1);

%Below is the the relevant generator that is used to subtract
%noise/variance from the path integral so as to match monte
%carlo density. Please note that this uses the gaussian
%increment between two grid points at different times that we
%calculated at start.
dfdtGenerator(mm,wnStart:Nn)=-(dfdtz01(mm)).*Zt(mm,wnStart:Nn)-1*dfdt2z02(mm).*(Zt(mm,wnStart:Nn).^2-1);
%Below is the evolution equation for evolution of probability
%density. Ideally Pprev(mm)=Pnew(mm)=ZProb(mm) that we started
%with since probability in each subdivision does not change but
%these Pnew and ZProb can be different at both extreme ends and
%sometimes results would be better when we would use Pnew to
%calculate the original variables from expectations by division
%of probability since the expectation of original variables
%used the same transition probabilities as were used for Pnew.
Pnew(wnStart:Nn)=Pnew(wnStart:Nn)+Pprev(mm).*pt(mm,wnStart:Nn);
%Below is evolution of expectation of yy(t) in original
%coordinates in our set up. This is not needed for anything but
%is a useful exercise
Eyynew(wnStart:Nn)=Eyynew(wnStart:Nn)+Eyyprev(mm).*pt(mm,wnStart:Nn)+ ...
dyyGenerator(mm,wnStart:Nn).*pt(mm,wnStart:Nn).*Pprev(mm);
%Below is evolution of expectation of yy(t)^omega1 in original
%coordinates in our set up. This is not needed for anything but
%is a useful exercise. This is counterpart of Ito change of
%variable
Efnew(wnStart:Nn)=Efnew(wnStart:Nn)+Efprev(mm).*pt(mm,wnStart:Nn)+ ...
dfGenerator(mm,wnStart:Nn).*pt(mm,wnStart:Nn).*Pprev(mm);
%Below is the calculation of path integrals.
if((ss==1))
%Below calculate the function value and save it in the
%appropriate grid point at the end of transition increment.

Efdtnew(wnStart:Nn)=Efdtnew(wnStart:Nn)+ ...
(ds*.5*yy1(mm).^omega1+ds*.5*yy2(wnStart:Nn).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm);
%Below subtract the volatility needed. This uses both the
%gaussian increment and transition probabilities between
%grid points at different time.
Efdtnew(wnStart:Nn)=Efdtnew(wnStart:Nn)+ ...
1.0/pi^2*dfdtGenerator(mm,wnStart:Nn).*pt(mm,wnStart:Nn).*ZProb(mm);
elseif ((ss==Tt/Freq))
%Evolve the time integrated expected value in the same grid
%point till end.
Efdtnew(mm)=Efdtnew(mm)+Efdtprev(mm);
%Below calculate the function value and save it in the
%appropriate grid point.
%Below subtract the volatility needed. This uses both the
%gaussian increment and transition probabilities between
%grid points at different time.
Efdtnew(wnStart:Nn)=Efdtnew(wnStart:Nn)+ ...
(ds*.5*yy1(mm).^omega1+ds*.5*yy2(wnStart:Nn).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm) + ...
1.0/pi^2*dfdtGenerator(mm,wnStart:Nn).*pt(mm,wnStart:Nn).*ZProb(mm);
%Efdtnew(wnStart:Nn)=Efdtnew(wnStart:Nn)+ds*(.5*yy2(wnStart:Nn).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm)+ ...
%    .50/Div0*dfdtGenerator(mm,wnStart:Nn).*pt(mm,wnStart:Nn).*ZProb(mm);

else
%Evolve the time integrated expected value in the same grid
%point till end. It will continue to evolve in the same grid
%point till end.

Efdtnew(mm)=Efdtnew(mm)+Efdtprev(mm);
%Below calculate the function value and save it in the
%appropriate grid point.
%Below subtract the volatility needed. This uses both the
%gaussian increment and transition probabilities between
%grid points at different time.

Efdtnew(wnStart:Nn)=Efdtnew(wnStart:Nn)+ ...
(ds*.5*yy1(mm).^omega1+ds*.5*yy2(wnStart:Nn).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm) + ...
1.0/pi^2*dfdtGenerator(mm,wnStart:Nn).*pt(mm,wnStart:Nn).*ZProb(mm);
end

end
Pprev(1:Nn)=Pnew(1:Nn);
Efprev(1:Nn)=Efnew(1:Nn);
Eyyprev(1:Nn)=Eyynew(1:Nn);
Efdtprev(1:Nn)=Efdtnew(1:Nn);
Pnew(1:Nn)=0;
Efnew(1:Nn)=0.0;
Eyynew(1:Nn)=0.0;
Efdtnew(1:Nn)=0.0;

%Below calculate the start and end time t1A and t2A for the next
%transition probability time intervals.
t1A=(ss)*ds;
t2A=(ss+1)*ds;
%Below set the starting value W1 for next time interval equal to end
%value of the current time interval W2 and calculate various
%integrals. All these integrals only need the starting value W1
W1(1:Nn)=W2(1:Nn);
%[fdtMu0dt,dfdtz01,dfdt2z02]=CalculateFdtDriftAndVol08(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds,t1A,t2A);
[wMu0dt2,dwMu0dt2,dw02]=CalculateDriftAndVol(mu1,mu2,beta1,beta2,gamma,sigma0,W1,1,Nn,ds);
[yyMu0dt,dyyz0,dyy2z0]=CalculateFunctionDriftAndVol(1.0,1.0,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds);
[fMu0dt,dfz0,df2z0]=CalculateFunctionDriftAndVol(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds);
[fdtMu0dt,dfdtz01,dfdt2z02]=CalculateFdtDriftAndVol08(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds,t1A,t2A);

end

end
Pnew=Pprev;
Efnew=Efprev;
Eyynew=Eyyprev;
Efdtnew=Efdtprev;

yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
%Below calculate value of yy obtained from expected value in transition
%probabilities framework after division of expected value with
%integrated probability. This value is called yyTr
yyTr(wnStart:Nn)=Eyynew(wnStart:Nn)./Pnew(wnStart:Nn)  ;

%plot((wnStart:Nn),yy(wnStart:Nn),'g',(wnStart:Nn),yyTr(wnStart:Nn),'b');

%str=input('Look at graph 01');

%Below calculate value of yy(T)^omega1 obtained from its expected value in
%transition probabilities framework after division of expected value with
%integrated probability. This value is called ffy

ffy(wnStart:Nn)=(Efnew(wnStart:Nn))./Pnew(wnStart:Nn);

%plot((wnStart:Nn),theta1.*yy(wnStart:Nn).^omega1,'g',(wnStart:Nn),ffy(wnStart:Nn),'b');
%str=input('Look at graph 02');

%Convert the expected values in each cell to standard Values
%associated with the grid cell after removing the
%integrated probability in the cell.Above for fy ito process. below
%for arithmetic sum path integral process.

%Below we repeat the process for path integral density
fydt(wnStart:Nn)=Efdtnew(wnStart:Nn)./ZProb(wnStart:Nn);

%below D's (the names of variables starting with D) are
%change of probability derivatives.
%Dfy_w(wmStart:wmEnd)=0;
Dffy(1:Nn)=0;
Dfydt(1:Nn)=0;
DyyTr(wnStart:Nn)=0.0;
for mm=wnStart+1:Nn-1
Dffy(mm) = (ffy(mm + 1) - ffy(mm - 1))/(Z(mm + 1) - Z(mm - 1));
Dfydt(mm) = (fydt(mm + 1) - fydt(mm - 1))/(Z(mm + 1) - Z(mm - 1));
DyyTr(mm) = (yyTr(mm + 1) - yyTr(mm - 1))/(Z(mm + 1) - Z(mm - 1));
end

%below pfy and pfydt are respective density amplitudes.
pfy(1:Nn)=0;
pfydt(1:Nn)=0;
pyyTr(1:Nn)=0.0;
for mm = wnStart+1:Nn-1
pfy(mm) = (normpdf(Z(mm),0, 1))/abs(Dffy(mm));
pfydt(mm) = (normpdf(Z(mm),0, 1))/abs(Dfydt(mm));
pyyTr(mm) = Pnew(mm)/abs(yyTr(mm+1)-yyTr(mm-1))*2;
end

y_w(1:Nn)=0;
y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
Dfy_w(wnStart:Nn)=0;
Dfw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
Dfy_w(nn) = (y_w(nn + 1) - y_w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
Dfw(nn) = (w(nn + 1) - w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
py_w(1:Nn)=0;
pw(1:Nn)=0;
for nn = wnStart:Nn-1
py_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
pw(nn) = (normpdf(Z(nn),0, 1))/abs(Dfw(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

theta1=1;
rng(29079137, 'twister')
paths=200000;
YY(1:paths)=x0;  %Original process monte carlo.
fYYdt(1:paths)=0.0;
Random1(1:paths)=0;
YYMean(1:TtM)=0;
for tt=1:TtM
t1M=(tt-1)*dtM;
t2M=tt*dtM;
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;

YYPrev(1:paths)=YY(1:paths);

fYYdt(1:paths)=fYYdt(1:paths)+.5*dtM*YY(1:paths).^omega1;

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);

%fYYdt(1:paths)=fYYdt(1:paths)+ t2M.*YY(1:paths).^omega1-t1M.*YYPrev(1:paths).^omega1- ...
%5(omega1.*YYPrev(1:paths).^(omega1-1).*(mu1*YYPrev(1:paths).^beta1+ mu2*YYPrev(1:paths).^beta2) + ...
%    .5*omega1.*(omega1-1).*YYPrev(1:paths).^(omega1-2).*sigma0^2.*YYPrev(1:paths).^(2*gamma))*(t2M^2-t1M^2)/2.0 - ...
%    omega1.*YYPrev(1:paths).^(omega1-1).*(sigma0.*YYPrev(1:paths).^gamma).*sqrt((t2M^3.0-t1M^3.0)/3).*HermiteP1(2,1:paths);

fYYdt(1:paths)=fYYdt(1:paths)+.5*dtM*YY(1:paths).^omega1;

%Uncomment for fourth order monte carlo

%
%   YY(1:paths)=YY(1:paths) + ...
%       (YCoeff0(1,1,2,1).*YY(1:paths).^Fp(1,1,2,1)+YCoeff0(1,2,1,1).*YY(1:paths).^Fp(1,2,1,1)+ ...
%       YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1))*dtM + ...
%       (YCoeff0(1,1,3,1).*YY(1:paths).^Fp(1,1,3,1)+YCoeff0(1,2,2,1).*YY(1:paths).^Fp(1,2,2,1)+ ...
%       YCoeff0(2,1,2,1).*YY(1:paths).^Fp(2,1,2,1)+YCoeff0(1,3,1,1).*YY(1:paths).^Fp(1,3,1,1)+ ...
%       YCoeff0(2,2,1,1).*YY(1:paths).^Fp(2,2,1,1)+YCoeff0(3,1,1,1).*YY(1:paths).^Fp(3,1,1,1))*dtM^2 + ...
%       (YCoeff0(1,1,4,1).*YY(1:paths).^Fp(1,1,4,1)+YCoeff0(1,2,3,1).*YY(1:paths).^Fp(1,2,3,1)+ ...
%       YCoeff0(2,1,3,1).*YY(1:paths).^Fp(2,1,3,1)+YCoeff0(1,3,2,1).*YY(1:paths).^Fp(1,3,2,1)+ ...
%       YCoeff0(2,2,2,1).*YY(1:paths).^Fp(2,2,2,1)+YCoeff0(3,1,2,1).*YY(1:paths).^Fp(3,1,2,1)+ ...
%       YCoeff0(1,4,1,1).*YY(1:paths).^Fp(1,4,1,1)+YCoeff0(2,3,1,1).*YY(1:paths).^Fp(2,3,1,1)+ ...
%       YCoeff0(3,2,1,1).*YY(1:paths).^Fp(3,2,1,1)+YCoeff0(4,1,1,1).*YY(1:paths).^Fp(4,1,1,1))*dtM^3 + ...
%        (YCoeff0(1,1,5,1).*YY(1:paths).^Fp(1,1,5,1)+YCoeff0(1,2,4,1).*YY(1:paths).^Fp(1,2,4,1)+ ...
%     YCoeff0(2,1,4,1).*YY(1:paths).^Fp(2,1,4,1)+YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1)+ ...
%      YCoeff0(2,2,3,1).*YY(1:paths).^Fp(2,2,3,1)+ ...
%      YCoeff0(3,1,3,1).*YY(1:paths).^Fp(3,1,3,1)+YCoeff0(1,4,2,1).*YY(1:paths).^Fp(1,4,2,1)+ ...
%      YCoeff0(2,3,2,1).*YY(1:paths).^Fp(2,3,2,1)+YCoeff0(3,2,2,1).*YY(1:paths).^Fp(3,2,2,1)+ ...
%      YCoeff0(4,1,2,1).*YY(1:paths).^Fp(4,1,2,1)+YCoeff0(1,5,1,1).*YY(1:paths).^Fp(1,5,1,1)+ ...
%      YCoeff0(2,4,1,1).*YY(1:paths).^Fp(2,4,1,1)+ ...
%      YCoeff0(3,3,1,1).*YY(1:paths).^Fp(3,3,1,1)+YCoeff0(4,2,1,1).*YY(1:paths).^Fp(4,2,1,1)+ ...
%      YCoeff0(5,1,1,1).*YY(1:paths).^Fp(5,1,1,1))*dtM^4+ ...
%       ((YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
%     (YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
%     YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5+ ...
%     (YCoeff0(1,1,3,2).*YY(1:paths).^Fp(1,1,3,2)+YCoeff0(1,2,2,2).*YY(1:paths).^Fp(1,2,2,2)+ ...
%     YCoeff0(2,1,2,2).*YY(1:paths).^Fp(2,1,2,2)+YCoeff0(1,3,1,2).*YY(1:paths).^Fp(1,3,1,2)+ ...
%     YCoeff0(2,2,1,2).*YY(1:paths).^Fp(2,2,1,2)+YCoeff0(3,1,1,2).*YY(1:paths).^Fp(3,1,1,2)).*dtM^2.5+ ...
%     (YCoeff0(1,1,4,2).*YY(1:paths).^Fp(1,1,4,2)+YCoeff0(1,2,3,2).*YY(1:paths).^Fp(1,2,3,2)+ ...
%     YCoeff0(2,1,3,2).*YY(1:paths).^Fp(2,1,3,2)+YCoeff0(1,3,2,2).*YY(1:paths).^Fp(1,3,2,2)+ ...
%     YCoeff0(2,2,2,2).*YY(1:paths).^Fp(2,2,2,2)+ YCoeff0(3,1,2,2).*YY(1:paths).^Fp(3,1,2,2)+ ...
%     YCoeff0(1,4,1,2).*YY(1:paths).^Fp(1,4,1,2)+YCoeff0(2,3,1,2).*YY(1:paths).^Fp(2,3,1,2)+ ...
%     YCoeff0(3,2,1,2).*YY(1:paths).^Fp(3,2,1,2)+YCoeff0(4,1,1,2).*YY(1:paths).^Fp(4,1,1,2)).*dtM^3.5) .*HermiteP1(2,1:paths) + ...
%     ((YCoeff0(1,1,1,3).*YY(1:paths).^Fp(1,1,1,3) *dtM) + ...
%     (YCoeff0(1,1,2,3).*YY(1:paths).^Fp(1,1,2,3)+YCoeff0(1,2,1,3).*YY(1:paths).^Fp(1,2,1,3)+ ...
%     YCoeff0(2,1,1,3).*YY(1:paths).^Fp(2,1,1,3)).*dtM^2+ ...
%     (YCoeff0(1,1,3,3).*YY(1:paths).^Fp(1,1,3,3)+YCoeff0(1,2,2,3).*YY(1:paths).^Fp(1,2,2,3)+ ...
%     YCoeff0(2,1,2,3).*YY(1:paths).^Fp(2,1,2,3) + YCoeff0(1,3,1,3).*YY(1:paths).^Fp(1,3,1,3)+ ...
%     YCoeff0(2,2,1,3).*YY(1:paths).^Fp(2,2,1,3)+YCoeff0(3,1,1,3).*YY(1:paths).^Fp(3,1,1,3)).*dtM^3).*HermiteP1(3,1:paths) + ...
%     ((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )+ ...
%     (YCoeff0(1,1,2,4).*YY(1:paths).^Fp(1,1,2,4)+YCoeff0(1,2,1,4).*YY(1:paths).^Fp(1,2,1,4)+ ...
%     YCoeff0(2,1,1,4).*YY(1:paths).^Fp(2,1,1,4))*dtM^2.5).*HermiteP1(4,1:paths) + ...
%     (YCoeff0(1,1,1,5).*YY(1:paths).^Fp(1,1,1,5)*dtM^2.0).*HermiteP1(5,1:paths);
%

%    YYMean(tt)=YYMean(tt)+sum(YY(1:paths))/paths;

end
YY(YY<0)=0;
disp('Original process average from monte carlo');
MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
disp('Original process average from our simulation');
ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates

disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average

fYYdtm=sum(fYYdt(:))/paths   %path integral average from monte carlo
fydtm= sum(fydt(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %path integral from transition

fYYdtvar=sum((fYYdt(:)-fYYdtm).^2)/paths %path integral variance from monte carlo
fydtvar= sum((fydt(wnStart+1:Nn-1)-fydtm).^2.*ZProb(wnStart+1:Nn-1)) %path integral variance from analytical work.

BinSize=.0075*1/2*2;%Please change this bin size accordingly. When data range is too small decrease it.
%When density range is too large increase it. Please change it accordingly
%for a good graph. This is important.
MaxCutOff=40;
[fYYdtDensity,IndexOutfYYdt,IndexMaxfYYdt] = MakeDensityFromSimulation_Infiniti(fYYdt,paths,BinSize,MaxCutOff);
plot(fydt(wnStart+1:Nn-1),pfydt(wnStart+1:Nn-1),'r',IndexOutfYYdt(1:IndexMaxfYYdt),fYYdtDensity(1:IndexMaxfYYdt),'g');
title(sprintf('Path Integral Density, x0 = %.2f,theta=%.2f,kappa=%.2f,gamma=%.2f,sigma=%.2f,T=%.2f,omega1=%.2f', x0,theta,kappa,gamma,sigma0,T,omega1));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

legend({'Path Integral Analytic Density','Path Integral Monte Carlo Density'},'Location','northeast')

str=input('red line is arithmetic dt integral density from transition probability framework , green is monte carlo.');

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_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',yyTr(wnStart+1:Nn-1),pyyTr(wnStart+1:Nn-1),'b');
%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));

%title('text',sprintf('x0 = %f', x0),sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T))
%xlabel(a)
%title('Graphical Comparison of Ito-Hermite Method and Monte Carlo Method Density')
legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')

str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');

%plot((wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r');
end
.
Here is the updated function you will also need to run this program.
.
function [fMu0dt,dfz01,df2z02]=CalculateFdtDriftAndVol08(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,w,wnStart,Nn,dt,t1,t2)

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

%f(wnStart:Nn)=theta1*yy(wnStart:Nn).^(omega1);
df(wnStart:Nn)=omega1*theta1*yy(wnStart:Nn).^(omega1-1);
DD(wnStart:Nn)=mu1*yy(wnStart:Nn).^beta1+mu2*yy(wnStart:Nn).^beta2;
d2f(wnStart:Nn)=theta1*omega1*(omega1-1).*yy(wnStart:Nn).^(omega1-2);

QQ(wnStart:Nn)=sigma0^2*yy(wnStart:Nn).^(2*gamma);
d_dfDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*yy(wnStart:Nn).^(omega1+beta1-2)+ ...
mu2*omega1*theta1*(omega1+beta2-1).*yy(wnStart:Nn).^(omega1+beta2-2);
d2_dfDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*(omega1+beta1-2).*yy(wnStart:Nn).^(omega1+beta1-3)+ ...
mu2*omega1*theta1*(omega1+beta2-1).*(omega1+beta2-2).*yy(wnStart:Nn).^(omega1+beta2-3);

d_d2fQQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*yy(wnStart:Nn).^(omega1+2*gamma-3);

d2_d2fQQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*(omega1+2*gamma-3).*yy(wnStart:Nn).^(omega1+2*gamma-4);

VV(wnStart:Nn)=sigma0.*yy(wnStart:Nn).^gamma;

d_dfVV(wnStart:Nn)= ...
sigma0.*omega1.*theta1.*(omega1+gamma-1).*yy(wnStart:Nn).^(omega1+gamma-2);

d2_dfVV(wnStart:Nn)= ...
sigma0.*omega1.*theta1.*(omega1+gamma-1).*(omega1+gamma-2).*yy(wnStart:Nn).^(omega1+gamma-3);

%    fMu0dt1(wnStart:Nn)=f(wnStart:Nn)*dt;%+ ...

fMu0dt(wnStart:Nn)=(df(wnStart:Nn).*DD(wnStart:Nn) + ...
.5*d2f(wnStart:Nn).*QQ(wnStart:Nn))*(t2^2-t1^2)/2.0+ ...
(d_dfDD(wnStart:Nn).*DD(wnStart:Nn) + ...
.5*d_d2fQQ(wnStart:Nn).*DD(wnStart:Nn) + ...
.5*d2_dfDD(wnStart:Nn).*QQ(wnStart:Nn) + ...
.25*d2_d2fQQ(wnStart:Nn).*QQ(wnStart:Nn)).*((t2^3-t1^3)/3-(t1*t2^2-t1^3)/2.0);

dfz01(1:Nn)=df(wnStart:Nn).*VV(wnStart:Nn).*(sqrt(t2^3.0-t1^3.0)/sqrt(3))+ ...
1* (d_dfVV(wnStart:Nn).*DD(wnStart:Nn)+ ...
.5*d2_dfVV(wnStart:Nn).*QQ(wnStart:Nn)).*(sqrt(t2^5/5-t1^5/5-t1^2*t2^3/3+t1^5/3)) + ...
1* (d_dfDD(wnStart:Nn).*VV(wnStart:Nn) + ...
.5*d_d2fQQ(wnStart:Nn).*VV(wnStart:Nn)).*sqrt((t2^5.0-t1^5.0)/4*(1-1/(5))-(t2^4*t1-t1^5)/4);

df2z02(1:Nn)=1*d_dfVV(wnStart:Nn).*VV(wnStart:Nn).*sqrt((t2^4/4-t1^4/4-t1*t2^3/3+t1^4/3))/sqrt(2);

end
.
And here is the output of this program
MCMean =
0.154865203933743
Original process average from our simulation
ItoHermiteMean =
0.155180606064941
true Mean only applicble to standard SV mean reverting type models otherwise disregard
TrueMean =
0.155181916175716
fYYdtm =
0.194741089398704
fydtm =
0.194805045257125
fYYdtvar =
0.004987629740974
fydtvar =
0.005387533097219
IndexMax =
144

and here is the output graph for this program.

I have started to make the changes that I mentioned at the start of this post and hope to be able to upload a new program with most issues resolved by the end of this week.

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

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

Friends, I just made this quick change in calculation of transition probabilities and used error function for calculation of transition probabilities and got better results. I am posting the new program. This is still intermediate and will undergo many changes but I decided to post it right away.
Here is the code.
.
function [] = FPERevisitedTransProb08()

%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%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*4/4*2*2;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
%Tt=8*60;
T=Tt*dt;
OrderA=4;  %
OrderM=4;  %
%dtM=.125/8;%Monte carlo time interval size dtM.
%TtM=8*8;%Monte carlo number of simulation intervals.
dtM=dt*4*2;
TtM=Tt/4/2;

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

%theta=mu1/(-mu2);
%kappa=-mu2;

x0=.250;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.95;%50;   % volatility power.
kappa=2.0;%.950;   %mean reversion parameter.
theta=.10;%mean reversion target
sigma0=.70;%Volatility value
omega1=1.0; %power on yy(t) for fy and fyPI.
%This is yy(t)^omega1 is the function that is summed in path integral
%I have tested powers between omega1=.5 and omega1=2;
theta1=1;%This is the weight on arithmetic sum in the path integral.
%keep it equal to one otherwise results will be erroneous.
%This will be activated in future and made time dependent.

%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)-10.5)*dNn-NnMid);

Z
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);
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.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ss=0;  %ss is index for transition probability calculation time steps.
Freq=2.0;%transition probabilities are calculated Freq time intervals apart.
wnStart=1;%
d2wdZ2(1:Nn)=0;
d2wdZ2A(1:Nn)=0;
dwdZ(1:Nn)=0;
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:Nn)=0.0;
Efprev(1:Nn)=theta1*yy(1:Nn).^omega1.*ZProb(1:Nn);%value of function at
%starting nodes.
Efnew(1:Nn)=0;
Eyyprev(1:Nn)=yy(1:Nn).*ZProb(1:Nn);%Value of SDE variable in original
%coordinates at starting nodes.
Eyynew(1:Nn)=0;
Efdtprev(1:Nn)=0;%Value of path-integral at starting nodes is zero.
Efdtnew(1:Nn)=0;

tic

for tt=1:Tt
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
[wMu0dt,dwMu0dtdw,c1] = CalculateDriftAndVolA4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn) ;% ...
%dw(wnStart:Nn)=sigma0.*sqrt(dt).*Z(wnStart:Nn) ;
dw2(wnStart:Nn)=dw(wnStart:Nn).^2;
w(isnan(w)==1)=0;
wMu0dt(isnan(wMu0dt)==1)=0;
%wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%     [dwdZ,d2wdZ2A,d3wdZ3A] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z);
%     d2wdZ2(wnStart:Nn)=0.0;
%     d3wdZ3(wnStart:Nn)=0.0;
%     if (tt>36)
%         d3wdZ3(wnStart:Nn) = (d3wdZ3A(wnStart:Nn));
%         d2wdZ2(wnStart:Nn)=d2wdZ2A(wnStart:Nn);
%         d2wdZ2(1:15)=d2wdZ2(1:15).*exp(-12*kappa^2*(1-gamma).*dNn*(15-(1:15)));
%         d3wdZ3(1:15)=d3wdZ3(1:15).*exp(-12*kappa^2*(1-gamma).*dNn*(15-(1:15)));
%
%     end

%tt
%plot((wnStart:Nn),d2wdZ2A(wnStart:Nn),'g',(wnStart:Nn),d2wdZ2(wnStart:Nn),'r');
%str=input('Look at d2wdZ2A(wnStart:Nn)');

%If the calculation of the program blows, it means that above calculated
%derivative is unstable or bad. Plese replace with your favourite method of
%calculating this second derivative.

dZdw(wnStart:Nn)=1.0./dwdZ(wnStart:Nn);
if(tt==1)
wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));

w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+ ...
sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)));
end
if(tt>1)

C22(wnStart:Nn)=.0125/pi*25.0;
C33(wnStart:Nn)=.0125/pi*2;
C44(wnStart:Nn)=.0125/pi*25.0;

C22(wnStart:Nn)=.0125/pi*25.0*1.5;
C33(wnStart:Nn)=.0125/pi*2*1.0;
C44(wnStart:Nn)=.0125/pi*25.0*1.5;

TermA0(wnStart:Nn)=-C22(wnStart:Nn).*3.*Z(wnStart:Nn).^1.*dZdw(wnStart:Nn).*d2wdZ2(wnStart:Nn)*sigma0^2.*dt + ...
+C44(wnStart:Nn).*(-d3wdZ3(wnStart:Nn).*dZdw(wnStart:Nn)+3*d2wdZ2(wnStart:Nn).^2.*dZdw(wnStart:Nn).^2)*sigma0^2.*dt+ ...
-C33(wnStart:Nn).*dwMu0dtdw(wnStart:Nn).*dwdZ(wnStart:Nn).^2;

TermA(wnStart:Nn)= sqrt(abs(TermA0(wnStart:Nn)));
SignTermA(wnStart:Nn)=sign(TermA0(wnStart:Nn));

w(wnStart:Nn)=w(wnStart:Nn)+wMu0dt(wnStart:Nn)+(-(Z(wnStart:Nn).*dwdZ(wnStart:Nn))+d2wdZ2(wnStart:Nn))+ ...
sign((Z(wnStart:Nn).*dwdZ(wnStart:Nn))-d2wdZ2(wnStart:Nn)+dw(wnStart:Nn)+SignTermA(wnStart:Nn).*TermA(wnStart:Nn)).* ...
sqrt(abs(sign((Z(wnStart:Nn).*dwdZ(wnStart:Nn))-d2wdZ2(wnStart:Nn)).* ...
((Z(wnStart:Nn).*dwdZ(wnStart:Nn))-d2wdZ2(wnStart:Nn)).^2+ ...
sign(+dw(wnStart:Nn)).*dw2(wnStart:Nn)+ ...
SignTermA(wnStart:Nn).*TermA(wnStart:Nn).^2));% + ...
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Below are instability correction equations.
%Interpolation of three points in instability region with lagrange
%polynomials.

w(NnMidl) = InterpolateOrderN8(8,Z(NnMidl),Z(NnMidl-4),Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidh+2),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),w(NnMidl-4),w(NnMidl-3),w(NnMidl-2),w(NnMidl-1),w(NnMidh+2),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5));
w(NnMidh) = InterpolateOrderN8(8,Z(NnMidh),Z(NnMidl-4),Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidh+2),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),w(NnMidl-4),w(NnMidl-3),w(NnMidl-2),w(NnMidl-1),w(NnMidh+2),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5));
w(NnMidh+1) = InterpolateOrderN8(8,Z(NnMidh+1),Z(NnMidl-4),Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidh+2),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),w(NnMidl-4),w(NnMidl-3),w(NnMidl-2),w(NnMidl-1),w(NnMidh+2),w(NnMidh+3),w(NnMidh+4),w(NnMidh+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));

w1(1:Nn-1)=w(1:Nn-1);
w1(Nn)=w(Nn);
w2(1:Nn-1)=w(2:Nn);
w2(Nn)=wE;
w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
%Change 3:7/25/2020: I have improved zero correction in above.
w(w<0)=0.0;
for nn=1:Nn
if(w(nn)<=0)
wnStart=nn+1;
end
end

%         %%1st order mean correction. We know the mean in original
%         %%coordinates so I shift the density into original coordinates,
%         %%apply the mean correction and then transform again into Lamperti
%         %%coordinates. Algebra solves two equations in two unknowns.
%         %%Equations are Sum_0^N{(Y_w(wnStart:Nn)-Y0)*W0}=1 and the second
%         %%equation is Sum_0^N{Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0}=u0
%         %%Two unknows are Y0 and W0. u0 is known mean.
%   tt;
u0=theta+(x0-theta)*exp(-kappa*(tt*dt)); %analytic mean of the density
%
%         %If you are not using stochastic volatility, replace above with
%         %true mean otherwise results would become garbage.
%
Y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
%
Yn2Pn=sum(Y_w(wnStart:Nn).*Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
YnPn=sum(Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
Pn=1.0;%%Sum(ZProb(1:Nn))
Y0=(Yn2Pn-u0*YnPn)/(YnPn-u0*Pn);

Y0Pn=Y0*Pn;
W0=1/(YnPn-Y0Pn);
YCorrected(wnStart:Nn)=Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0;
wCorrected(wnStart:Nn)=(YCorrected(wnStart:Nn).^(1-gamma))./(1-gamma);
w(wnStart:Nn)=wCorrected(wnStart:Nn);
%
%I have disabled the mean correction. The above mean correction is only valid for
%standard mean reverting stochastic volatility type SDEs. To enable mean
%correction please uncomment the above block.

[dwdZ,d2wdZ2A,d3wdZ3A] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z);
d2wdZ2(wnStart:Nn)=0.0;
d3wdZ3(wnStart:Nn)=0.0;
if (tt>36)
d3wdZ3(wnStart:Nn) = (d3wdZ3A(wnStart:Nn));
d2wdZ2(wnStart:Nn)=d2wdZ2A(wnStart:Nn);
d2wdZ2(1:15)=d2wdZ2(1:15).*exp(-12*kappa^2*(1-gamma).*dNn*(15-(1:15)));
d3wdZ3(1:15)=d3wdZ3(1:15).*exp(-12*kappa^2*(1-gamma).*dNn*(15-(1:15)));

end

%Start of functions and path integrals calculations
%For transition probabilities calculations, the time grid is Freq intervals apart.
% so time interval ds for transition probability calculation is
% ds=dt*Freq; and t1A is start time for transition probabilities interval
% while t2A is end time of transition probabilities calculation interval.
% so t2A=t1A+ds; W1 is the bessel coordinates value w at t1A and W2 is
% equal to w at t2A. W1=w(t1A) and W2=w(t2A). W1 and W2 are used for
% transition probabilities calculations.

if(tt==1)

W1(1:Nn)=x0^(1-gamma)/(1-gamma);

ds=dt*Freq;
t1A=0;
t2A=ds;
%below calculates integrals for drift and volatility for transition
%probabilities.
[wMu0dt2,dwMu0dt2,dw02]=CalculateDriftAndVol(mu1,mu2,beta1,beta2,gamma,sigma0,W1,1,Nn,ds);
%below calculates integrals for drift and volatility for yy(t).
%yy(t) is the SDE variable in original coordinates as opposed to
%bessel coordinates.
[yyMu0dt,dyyz0,dyy2z0]=CalculateFunctionDriftAndVol(1.0,1.0,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds);
%below calculates integrals for drift and volatility for
%yy(t)^omega1 which is used in Ito change of variable calculations
%on the lattice. This is not needed for path integrals.

[fMu0dt,dfz0,df2z0]=CalculateFunctionDriftAndVol(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds);

%below calculates integrals for drift and volatility for
%yy(s)^omega1 ds. which is used in path integrals calculations
%on the lattice.
%[fdtMu0dt,dfdtz01,dfdt2z02]=CalculateFdtDriftAndVol08(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds,t1A,t2A);

[fdtMu0dt,dfdtz01,dfdt2z02]=CalculateFdtDriftAndVol08(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds,t1A,t2A);

end

if(rem(tt,Freq)==0)

%When this above condition is true it means that we are at end of
%particular transition probability interval and we have end of
%interval W2(s)=w(t) while we will use the value W1 that was
%previously saved Freq dt time intervals earlier or ds=Freq*dt
%time ago.

ss=ss+1; %transition loop variable and is one at the end of first
%transition probability interval.
%Div0=Div0+.5*ds;

W2(1:Nn)=w(1:Nn);
ds=dt*Freq;
%(1*f[i-2]-8*f[i-1]+0*f[i+0]+8*f[i+1]-1*f[i+2])/(12*1.0*h**1)
dW2(3:Nn-2)=(W2(1:Nn-4)-8*W2(2:Nn-3)+8*W2(4:Nn-1)-W2(5:Nn))/12.0;
dW2(2)=(W2(3)-W2(1))/2.0;
dW2(Nn-1)=(W2(Nn)-W2(Nn-2))/2.0;

%    dW2(2:Nn-1)=(W2(3:Nn)-W2(1:Nn-2))/2.0;%For calculation of width of
%interval for integrated probability.
%      dW2(1)=dW2(2);%rough approximation, will change
%      dW2(Nn)=dW2(Nn-1);%rough approximation, will change
dW2(1)=ZProb(1)./ZProb(2).*dW2(2);
dW2(Nn)=ZProb(Nn)./ZProb(Nn-1).*dW2(Nn-1);

yy1(1:Nn)=((1-gamma)*W1(1:Nn)).^(1/(1-gamma));%Value in origianl
%coordinates at the start of transition probabilities interval.
yy2(1:Nn)=((1-gamma)*W2(1:Nn)).^(1/(1-gamma));%Value in origianl
%coordinates at the end of transition probabilities interval.
fyy1(1:Nn)=theta1.*yy1(1:Nn).^omega1;%function value at the start
%both for Ito calculation of function and path integrals.
fyy2(1:Nn)=theta1.*yy2(1:Nn).^omega1;%function value at the end
%both for Ito calculation of function and path integrals.

[W2GridStarts,W2GridEnds] = CalculateGridStartsAndEnds(W2,Z,wnStart,Nn,dNn);

for mm=1:Nn
%Below is the gaussian increment between two grid points W1(mm)
% at time t1A and W2(wnStart:Nn) at time t2A, backed
%out from the original stochastic differential equation drift
%and volatility. This gaussian increment is used for
%probability calculation between corresponding grid points and
%also for calculation of stochastic integrals between the
%corresponding grid points. Please note that this gaussian
%increment remains the same in bessel coordinates and the
%original coordinates
Zt(mm,wnStart:Nn)=(W2(wnStart:Nn)-W1(mm)-1*wMu0dt2(mm))/dw02(mm);
%Below is calculation of transition probability between two
%grid points W1(mm)  at time t1A and W2(wnStart:Nn) at time t2A
%using the gaussian increment preeviously calculated.
%            pt(mm,wnStart:Nn)=exp(-.5* Zt(mm,wnStart:Nn).^2)./sqrt(2*pi*dw02(mm).^2).*dW2(wnStart:Nn);%-exp(-.5* Zt(mm,wnStart:Nn).^2)./sqrt(2*pi*dw02(mm).^2)./dw02(mm).* Zt(mm,wnStart:Nn).*dW2(wnStart:Nn).^2/8;

pt(mm,wnStart:Nn)=0.5* erf((0.707107.*W1(mm)+0.707107*wMu0dt2(mm) - 0.707107* W2GridStarts(wnStart:Nn))/dw02(mm)) - ...
0.5* erf((0.707107* W1(mm) + 0.707107* wMu0dt2(mm) - 0.707107 *W2GridEnds(wnStart:Nn))/dw02(mm));

pt(isnan(pt(mm,:))==1)=0;
pt(isinf(pt(mm,:))==1)=0;

% % % %If true mean is not known in original coordinates just comment
% % %        %the if loop below. Currently, it uses a mean-reverting mean SV models.
% % % %         if( (b1(mm)<wE0-.05*sigma0*sqrt(dt))&& (Pn(mm)>0) )
%                Pn(mm)=sum(pt(mm,1:Nn));
%                Pn(mm)=1.0;
%                pt0(1:Nn)=pt(mm,1:Nn);
% % % %             if(mm==6)
% % % %                 pt0(1:c2N(mm))=pt(mm,1:c2N(mm));
% % % %                 Pn(mm)=Pn(mm);
% % % %             end
% % % %             if(mm<6)
% % % %                 pt0(1:c2N(mm))=pt(mm,1:c2N(mm));
% % % %                 Pn(mm)=Pn(mm);
% % % %             end
%             TrueMean=theta+(yy(mm)-theta)*exp(-kappa*(ds));
%             ptCorrected0 = CorrectProbAndMean(TrueMean,1,Nn,yy,pt0,Pn(mm));
%   pt(mm,1:Nn)=ptCorrected0(1:Nn);
% % % %
% % % %             %In this if loop, the transition probabilities are corrected so
% % % %             %that sum of next stage target probabilities sum to one and
% % % %             %mean equals the process mean at next stage given its evolution
% % % %             %from point w(mm).
% % % %             %If true mean is not known in original coordinates just comment
% % % %             %this if loop.
% % % %
% % % %
% % % %      Pn(mm)
% % % %      sum(pt(mm,c1N(mm):c2N(mm)))
% % % % %     P2n(mm)
% % % %     sum(yy(c1N(mm):c2N(mm)).*pt(mm,c1N(mm):c2N(mm)))
% % % %
% % % % %    yy(mm)
% % % %     theta+(yy(mm)-theta)*exp(-kappa*dt)
% % % %     mm
% % % %   %  str=input('Look at values') ;
% % % %

%Below is the counterpart of original SDE that is used for the
%evolution of expcted value of yy(t). It is totally independent
%for its own sake and is not needed for any function or
%path integral calculations.
dyyGenerator(mm,wnStart:Nn)=yyMu0dt(mm)+dyyz0(mm).*Zt(mm,wnStart:Nn)+dyy2z0(mm).*(Zt(mm,wnStart:Nn).^2-1);

%Below is the Ito-generator for fy=theta1*yy(t).^omega1 which
%is used for the evolution of expected value of fy and this is
%totally independent for its own sake and is not needed in any
%other calculations.
dfGenerator(mm,wnStart:Nn)=fMu0dt(mm)+dfz0(mm).*Zt(mm,wnStart:Nn)+df2z0(mm).*(Zt(mm,wnStart:Nn).^2-1);

%Below is the the relevant generator that is used to subtract
%noise/variance from the path integral so as to match monte
%carlo density. Please note that this uses the gaussian
%increment between two grid points at different times that we
%calculated at start.
dfdtGenerator(mm,wnStart:Nn)=-(dfdtz01(mm)).*Zt(mm,wnStart:Nn)-1*dfdt2z02(mm).*(Zt(mm,wnStart:Nn).^2-1);
%Below is the evolution equation for evolution of probability
%density. Ideally Pprev(mm)=Pnew(mm)=ZProb(mm) that we started
%with since probability in each subdivision does not change but
%these Pnew and ZProb can be different at both extreme ends and
%sometimes results would be better when we would use Pnew to
%calculate the original variables from expectations by division
%of probability since the expectation of original variables
%used the same transition probabilities as were used for Pnew.
Pnew(wnStart:Nn)=Pnew(wnStart:Nn)+Pprev(mm).*pt(mm,wnStart:Nn);
%Below is evolution of expectation of yy(t) in original
%coordinates in our set up. This is not needed for anything but
%is a useful exercise
Eyynew(wnStart:Nn)=Eyynew(wnStart:Nn)+Eyyprev(mm).*pt(mm,wnStart:Nn)+ ...
dyyGenerator(mm,wnStart:Nn).*pt(mm,wnStart:Nn).*Pprev(mm);
%Below is evolution of expectation of yy(t)^omega1 in original
%coordinates in our set up. This is not needed for anything but
%is a useful exercise. This is counterpart of Ito change of
%variable
Efnew(wnStart:Nn)=Efnew(wnStart:Nn)+Efprev(mm).*pt(mm,wnStart:Nn)+ ...
dfGenerator(mm,wnStart:Nn).*pt(mm,wnStart:Nn).*Pprev(mm);
%Below is the calculation of path integrals.
if((ss==1))
%Below calculate the function value and save it in the
%appropriate grid point at the end of transition increment.

Efdtnew(wnStart:Nn)=Efdtnew(wnStart:Nn)+ ...
(ds*.5*yy1(mm).^omega1+ds*.5*yy2(wnStart:Nn).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm);
%Below subtract the volatility needed. This uses both the
%gaussian increment and transition probabilities between
%grid points at different time.
Efdtnew(wnStart:Nn)=Efdtnew(wnStart:Nn)+ ...
1.0/pi^2*dfdtGenerator(mm,wnStart:Nn).*pt(mm,wnStart:Nn).*ZProb(mm);
elseif ((ss==Tt/Freq))
%Evolve the time integrated expected value in the same grid
%point till end.
Efdtnew(mm)=Efdtnew(mm)+Efdtprev(mm);
%Below calculate the function value and save it in the
%appropriate grid point.
%Below subtract the volatility needed. This uses both the
%gaussian increment and transition probabilities between
%grid points at different time.
Efdtnew(wnStart:Nn)=Efdtnew(wnStart:Nn)+ ...
(ds*.5*yy1(mm).^omega1+ds*.5*yy2(wnStart:Nn).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm) + ...
1.0/pi^2*dfdtGenerator(mm,wnStart:Nn).*pt(mm,wnStart:Nn).*ZProb(mm);
%Efdtnew(wnStart:Nn)=Efdtnew(wnStart:Nn)+ds*(.5*yy2(wnStart:Nn).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm)+ ...
%    .50/Div0*dfdtGenerator(mm,wnStart:Nn).*pt(mm,wnStart:Nn).*ZProb(mm);

else
%Evolve the time integrated expected value in the same grid
%point till end. It will continue to evolve in the same grid
%point till end.

Efdtnew(mm)=Efdtnew(mm)+Efdtprev(mm);
%Below calculate the function value and save it in the
%appropriate grid point.
%Below subtract the volatility needed. This uses both the
%gaussian increment and transition probabilities between
%grid points at different time.

Efdtnew(wnStart:Nn)=Efdtnew(wnStart:Nn)+ ...
(ds*.5*yy1(mm).^omega1+ds*.5*yy2(wnStart:Nn).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm) + ...
1.0/pi^2*dfdtGenerator(mm,wnStart:Nn).*pt(mm,wnStart:Nn).*ZProb(mm);
end

end
Pprev(1:Nn)=Pnew(1:Nn);
Efprev(1:Nn)=Efnew(1:Nn);
Eyyprev(1:Nn)=Eyynew(1:Nn);
Efdtprev(1:Nn)=Efdtnew(1:Nn);
Pnew(1:Nn)=0;
Efnew(1:Nn)=0.0;
Eyynew(1:Nn)=0.0;
Efdtnew(1:Nn)=0.0;

%Below calculate the start and end time t1A and t2A for the next
%transition probability time intervals.
t1A=(ss)*ds;
t2A=(ss+1)*ds;
%Below set the starting value W1 for next time interval equal to end
%value of the current time interval W2 and calculate various
%integrals. All these integrals only need the starting value W1
W1(1:Nn)=W2(1:Nn);
%[fdtMu0dt,dfdtz01,dfdt2z02]=CalculateFdtDriftAndVol08(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds,t1A,t2A);
[wMu0dt2,dwMu0dt2,dw02]=CalculateDriftAndVol(mu1,mu2,beta1,beta2,gamma,sigma0,W1,1,Nn,ds);
[yyMu0dt,dyyz0,dyy2z0]=CalculateFunctionDriftAndVol(1.0,1.0,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds);
[fMu0dt,dfz0,df2z0]=CalculateFunctionDriftAndVol(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds);
[fdtMu0dt,dfdtz01,dfdt2z02]=CalculateFdtDriftAndVol08(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds,t1A,t2A);

end

end
Pnew=Pprev;
Efnew=Efprev;
Eyynew=Eyyprev;
Efdtnew=Efdtprev;

yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
%Below calculate value of yy obtained from expected value in transition
%probabilities framework after division of expected value with
%integrated probability. This value is called yyTr
yyTr(wnStart:Nn)=Eyynew(wnStart:Nn)./Pnew(wnStart:Nn)  ;

%plot((wnStart:Nn),yy(wnStart:Nn),'g',(wnStart:Nn),yyTr(wnStart:Nn),'b');

%str=input('Look at graph 01');

%Below calculate value of yy(T)^omega1 obtained from its expected value in
%transition probabilities framework after division of expected value with
%integrated probability. This value is called ffy

ffy(wnStart:Nn)=(Efnew(wnStart:Nn))./Pnew(wnStart:Nn);

%plot((wnStart:Nn),theta1.*yy(wnStart:Nn).^omega1,'g',(wnStart:Nn),ffy(wnStart:Nn),'b');
%str=input('Look at graph 02');

%Convert the expected values in each cell to standard Values
%associated with the grid cell after removing the
%integrated probability in the cell.Above for fy ito process. below
%for arithmetic sum path integral process.

%Below we repeat the process for path integral density
fydt(wnStart:Nn)=Efdtnew(wnStart:Nn)./ZProb(wnStart:Nn);

%below D's (the names of variables starting with D) are
%change of probability derivatives.
%Dfy_w(wmStart:wmEnd)=0;
Dffy(1:Nn)=0;
Dfydt(1:Nn)=0;
DyyTr(wnStart:Nn)=0.0;
for mm=wnStart+1:Nn-1
Dffy(mm) = (ffy(mm + 1) - ffy(mm - 1))/(Z(mm + 1) - Z(mm - 1));
Dfydt(mm) = (fydt(mm + 1) - fydt(mm - 1))/(Z(mm + 1) - Z(mm - 1));
DyyTr(mm) = (yyTr(mm + 1) - yyTr(mm - 1))/(Z(mm + 1) - Z(mm - 1));
end

%below pfy and pfydt are respective density amplitudes.
pfy(1:Nn)=0;
pfydt(1:Nn)=0;
pyyTr(1:Nn)=0.0;
for mm = wnStart+1:Nn-1
pfy(mm) = (normpdf(Z(mm),0, 1))/abs(Dffy(mm));
pfydt(mm) = (normpdf(Z(mm),0, 1))/abs(Dfydt(mm));
pyyTr(mm) = Pnew(mm)/abs(yyTr(mm+1)-yyTr(mm-1))*2;
end

y_w(1:Nn)=0;
y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
Dfy_w(wnStart:Nn)=0;
Dfw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
Dfy_w(nn) = (y_w(nn + 1) - y_w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
Dfw(nn) = (w(nn + 1) - w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
py_w(1:Nn)=0;
pw(1:Nn)=0;
for nn = wnStart:Nn-1
py_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
pw(nn) = (normpdf(Z(nn),0, 1))/abs(Dfw(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

theta1=1;
rng(29079137, 'twister')
paths=200000;
YY(1:paths)=x0;  %Original process monte carlo.
fYYdt(1:paths)=0.0;
Random1(1:paths)=0;
YYMean(1:TtM)=0;
for tt=1:TtM
t1M=(tt-1)*dtM;
t2M=tt*dtM;
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;

YYPrev(1:paths)=YY(1:paths);

fYYdt(1:paths)=fYYdt(1:paths)+.5*dtM*YY(1:paths).^omega1;

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);

%fYYdt(1:paths)=fYYdt(1:paths)+ t2M.*YY(1:paths).^omega1-t1M.*YYPrev(1:paths).^omega1- ...
%5(omega1.*YYPrev(1:paths).^(omega1-1).*(mu1*YYPrev(1:paths).^beta1+ mu2*YYPrev(1:paths).^beta2) + ...
%    .5*omega1.*(omega1-1).*YYPrev(1:paths).^(omega1-2).*sigma0^2.*YYPrev(1:paths).^(2*gamma))*(t2M^2-t1M^2)/2.0 - ...
%    omega1.*YYPrev(1:paths).^(omega1-1).*(sigma0.*YYPrev(1:paths).^gamma).*sqrt((t2M^3.0-t1M^3.0)/3).*HermiteP1(2,1:paths);

fYYdt(1:paths)=fYYdt(1:paths)+.5*dtM*YY(1:paths).^omega1;

%Uncomment for fourth order monte carlo

%
%   YY(1:paths)=YY(1:paths) + ...
%       (YCoeff0(1,1,2,1).*YY(1:paths).^Fp(1,1,2,1)+YCoeff0(1,2,1,1).*YY(1:paths).^Fp(1,2,1,1)+ ...
%       YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1))*dtM + ...
%       (YCoeff0(1,1,3,1).*YY(1:paths).^Fp(1,1,3,1)+YCoeff0(1,2,2,1).*YY(1:paths).^Fp(1,2,2,1)+ ...
%       YCoeff0(2,1,2,1).*YY(1:paths).^Fp(2,1,2,1)+YCoeff0(1,3,1,1).*YY(1:paths).^Fp(1,3,1,1)+ ...
%       YCoeff0(2,2,1,1).*YY(1:paths).^Fp(2,2,1,1)+YCoeff0(3,1,1,1).*YY(1:paths).^Fp(3,1,1,1))*dtM^2 + ...
%       (YCoeff0(1,1,4,1).*YY(1:paths).^Fp(1,1,4,1)+YCoeff0(1,2,3,1).*YY(1:paths).^Fp(1,2,3,1)+ ...
%       YCoeff0(2,1,3,1).*YY(1:paths).^Fp(2,1,3,1)+YCoeff0(1,3,2,1).*YY(1:paths).^Fp(1,3,2,1)+ ...
%       YCoeff0(2,2,2,1).*YY(1:paths).^Fp(2,2,2,1)+YCoeff0(3,1,2,1).*YY(1:paths).^Fp(3,1,2,1)+ ...
%       YCoeff0(1,4,1,1).*YY(1:paths).^Fp(1,4,1,1)+YCoeff0(2,3,1,1).*YY(1:paths).^Fp(2,3,1,1)+ ...
%       YCoeff0(3,2,1,1).*YY(1:paths).^Fp(3,2,1,1)+YCoeff0(4,1,1,1).*YY(1:paths).^Fp(4,1,1,1))*dtM^3 + ...
%        (YCoeff0(1,1,5,1).*YY(1:paths).^Fp(1,1,5,1)+YCoeff0(1,2,4,1).*YY(1:paths).^Fp(1,2,4,1)+ ...
%     YCoeff0(2,1,4,1).*YY(1:paths).^Fp(2,1,4,1)+YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1)+ ...
%      YCoeff0(2,2,3,1).*YY(1:paths).^Fp(2,2,3,1)+ ...
%      YCoeff0(3,1,3,1).*YY(1:paths).^Fp(3,1,3,1)+YCoeff0(1,4,2,1).*YY(1:paths).^Fp(1,4,2,1)+ ...
%      YCoeff0(2,3,2,1).*YY(1:paths).^Fp(2,3,2,1)+YCoeff0(3,2,2,1).*YY(1:paths).^Fp(3,2,2,1)+ ...
%      YCoeff0(4,1,2,1).*YY(1:paths).^Fp(4,1,2,1)+YCoeff0(1,5,1,1).*YY(1:paths).^Fp(1,5,1,1)+ ...
%      YCoeff0(2,4,1,1).*YY(1:paths).^Fp(2,4,1,1)+ ...
%      YCoeff0(3,3,1,1).*YY(1:paths).^Fp(3,3,1,1)+YCoeff0(4,2,1,1).*YY(1:paths).^Fp(4,2,1,1)+ ...
%      YCoeff0(5,1,1,1).*YY(1:paths).^Fp(5,1,1,1))*dtM^4+ ...
%       ((YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
%     (YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
%     YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5+ ...
%     (YCoeff0(1,1,3,2).*YY(1:paths).^Fp(1,1,3,2)+YCoeff0(1,2,2,2).*YY(1:paths).^Fp(1,2,2,2)+ ...
%     YCoeff0(2,1,2,2).*YY(1:paths).^Fp(2,1,2,2)+YCoeff0(1,3,1,2).*YY(1:paths).^Fp(1,3,1,2)+ ...
%     YCoeff0(2,2,1,2).*YY(1:paths).^Fp(2,2,1,2)+YCoeff0(3,1,1,2).*YY(1:paths).^Fp(3,1,1,2)).*dtM^2.5+ ...
%     (YCoeff0(1,1,4,2).*YY(1:paths).^Fp(1,1,4,2)+YCoeff0(1,2,3,2).*YY(1:paths).^Fp(1,2,3,2)+ ...
%     YCoeff0(2,1,3,2).*YY(1:paths).^Fp(2,1,3,2)+YCoeff0(1,3,2,2).*YY(1:paths).^Fp(1,3,2,2)+ ...
%     YCoeff0(2,2,2,2).*YY(1:paths).^Fp(2,2,2,2)+ YCoeff0(3,1,2,2).*YY(1:paths).^Fp(3,1,2,2)+ ...
%     YCoeff0(1,4,1,2).*YY(1:paths).^Fp(1,4,1,2)+YCoeff0(2,3,1,2).*YY(1:paths).^Fp(2,3,1,2)+ ...
%     YCoeff0(3,2,1,2).*YY(1:paths).^Fp(3,2,1,2)+YCoeff0(4,1,1,2).*YY(1:paths).^Fp(4,1,1,2)).*dtM^3.5) .*HermiteP1(2,1:paths) + ...
%     ((YCoeff0(1,1,1,3).*YY(1:paths).^Fp(1,1,1,3) *dtM) + ...
%     (YCoeff0(1,1,2,3).*YY(1:paths).^Fp(1,1,2,3)+YCoeff0(1,2,1,3).*YY(1:paths).^Fp(1,2,1,3)+ ...
%     YCoeff0(2,1,1,3).*YY(1:paths).^Fp(2,1,1,3)).*dtM^2+ ...
%     (YCoeff0(1,1,3,3).*YY(1:paths).^Fp(1,1,3,3)+YCoeff0(1,2,2,3).*YY(1:paths).^Fp(1,2,2,3)+ ...
%     YCoeff0(2,1,2,3).*YY(1:paths).^Fp(2,1,2,3) + YCoeff0(1,3,1,3).*YY(1:paths).^Fp(1,3,1,3)+ ...
%     YCoeff0(2,2,1,3).*YY(1:paths).^Fp(2,2,1,3)+YCoeff0(3,1,1,3).*YY(1:paths).^Fp(3,1,1,3)).*dtM^3).*HermiteP1(3,1:paths) + ...
%     ((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )+ ...
%     (YCoeff0(1,1,2,4).*YY(1:paths).^Fp(1,1,2,4)+YCoeff0(1,2,1,4).*YY(1:paths).^Fp(1,2,1,4)+ ...
%     YCoeff0(2,1,1,4).*YY(1:paths).^Fp(2,1,1,4))*dtM^2.5).*HermiteP1(4,1:paths) + ...
%     (YCoeff0(1,1,1,5).*YY(1:paths).^Fp(1,1,1,5)*dtM^2.0).*HermiteP1(5,1:paths);
%

%    YYMean(tt)=YYMean(tt)+sum(YY(1:paths))/paths;

end
YY(YY<0)=0;
disp('Original process average from monte carlo');
MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
disp('Original process average from our simulation');
ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates

disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average

fYYdtm=sum(fYYdt(:))/paths   %path integral average from monte carlo
fydtm= sum(fydt(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %path integral from transition

fYYdtvar=sum((fYYdt(:)-fYYdtm).^2)/paths %path integral variance from monte carlo
fydtvar= sum((fydt(wnStart+1:Nn-1)-fydtm).^2.*ZProb(wnStart+1:Nn-1)) %path integral variance from analytical work.

BinSize=.0075*1/2*2;%Please change this bin size accordingly. When data range is too small decrease it.
%When density range is too large increase it. Please change it accordingly
%for a good graph. This is important.
MaxCutOff=40;
[fYYdtDensity,IndexOutfYYdt,IndexMaxfYYdt] = MakeDensityFromSimulation_Infiniti(fYYdt,paths,BinSize,MaxCutOff);
plot(fydt(wnStart+1:Nn-1),pfydt(wnStart+1:Nn-1),'r',IndexOutfYYdt(1:IndexMaxfYYdt),fYYdtDensity(1:IndexMaxfYYdt),'g');
title(sprintf('Path Integral Density, x0 = %.2f,theta=%.2f,kappa=%.2f,gamma=%.2f,sigma=%.2f,T=%.2f,omega1=%.2f', x0,theta,kappa,gamma,sigma0,T,omega1));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

legend({'Path Integral Analytic Density','Path Integral Monte Carlo Density'},'Location','northeast')

str=input('red line is arithmetic dt integral density from transition probability framework , green is monte carlo.');

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_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',yyTr(wnStart+1:Nn-1),pyyTr(wnStart+1:Nn-1),'b');
%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));

%title('text',sprintf('x0 = %f', x0),sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T))
%xlabel(a)
%title('Graphical Comparison of Ito-Hermite Method and Monte Carlo Method Density')
legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')

str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');

%plot((wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r');
end
.
.
Here is a new function that you will need to run this program.
.
function [GridStarts,GridEnds] = CalculateGridStartsAndEnds(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;

x0=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;
[y0] = InterpolateOrderN6(N,x0,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
GridEnds(wnStart)=y0;

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;

x0=Z(Nn-1)+.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;
[y0] = InterpolateOrderN6(N,x0,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);

GridEnds(Nn-1)=y0;

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

end
.
Here is another function that I have also previously posted but older function will create problems and you have to use this new version. You need this function.
.
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

%y0=(x0-x3)/(x2-x3)*y2+(x0-x2)/(x3-x2)*y3;

end
.
.
Here is the output of the new program.
ItoHermiteMean =
0.120300287074824
true Mean only applicable to standard SV mean reverting type models otherwise disregard
TrueMean =
0.120300292485492
Original process average from monte carlo
MCMean =
0.120094622869520
Original process average from our simulation
ItoHermiteMean =
0.120300287074824
true Mean only applicble to standard SV mean reverting type models otherwise disregard
TrueMean =
0.120300292485492
fYYdtm =
0.164804351613393
fydtm =
0.164293039543594
fYYdtvar =
0.002161942305456
fydtvar =
0.002168980290649
IndexMax =
92

Here is a graph I got from this program.

This graph has the same parameters as in the previous post but I am using kappa=2 instead of kappa=1 in previous post.

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

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

Friends, I am taking a small break from stochastic integrals on transition probabilities grid for about two weeks. I have some other ideas about solution of stochastic differential equations using deterministic methods. We have already just done that using solution of fokker-planck equation. But I have a strong belief that we can solve SDEs on the same Z-grid using alternative deterministic methods using similar mathematics with hermite polynomials just like we did in monte carlo simulations.
I have already given solutions to forward starting iterated stochastic integrals here in a previous post
https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=975#p860189
Here is the workmap for the new method.
1. We will use similar Z-grid like we used in solution of Fokker-planck equation. Only on this grid, we will keep track of cumulative Z that will have an increasing variance with each step. We will advance on similar Z-specific lines using this method just as we did for fokker-planck equation.
2. Again we will have to keep track of cumulative Z(that has an increasing along time but uniform variance across the cross-section) and then I also believe we will be able to find transition probabilities for the SDE using just the increment required to connect Z(with increasing variance) at two points in time. I hope that we will be able to find transition between any two points in time very easily.
3. We will have to use the forward starting iterated stochastic integrals that I have shown in the above cited link to a previous post. We will just use simple hermite polynomials based calculus with forward starting iterated stochastic integrals.
4. There might be a way so that we could use brownian motion with increasing variance instead of simple Z with increasing variance as this could greatly simplify and help in finding transition probabilities everywhere.

I want to work on this because it will drastically simplify the solution to SDEs as compared to what we have to currently do for the solution to fokker-planck equation. And it will give us huge insights that we would be able to use when we solve for stochastic integrals on such grids. I even think that we might not even need transition probabilities for stochastic integrals and we might be able to solve them in a straightforward manner using some similar method that we used for the original SDE with a bit of Z-calculus..

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

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

Another benefit of the new method would be that we will be working in original coordinates everywhere unlike we did for fokker-planck equation yet. And we will be taking far larger time steps as compared to very small steps for fokker-planck equation..

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

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

In my current work, I have to deal with a lot of iterated integrals of the form $\int_{t_1}^{t_2} t dz(t) \int_{t_1}^{t} ds$, $\int_{t_1}^{t_2} t dz(t) \int_{t_1}^{t} dz(s)$ and $\int_{t_1}^{t_2} t dt \int_{t_1}^{t} dz(s)$ and even thrice or more iterated integrals like $\int_{t_1}^{t_2} t dz(t) \int_{t_1}^{t} dz(s) \int_{t_1}^{s} dv$. These integrals would not have been difficult if the lower limit of all the integrals involved is zero but once lower limit is different like we have $t_1$ in above integral, they become difficult. However using a trick, we can very easily solve all these integrals. Since we will be dealing with a lot of these integrals in further path integrals work, I decided to share how to solve these integrals with friends in this post. Using this technique, you would be able to solve a large number of difficult stochastic integrals.
Basically all we have to do is to break down the iterated integral into larger number of simpler integrals where the lower limit of all integrals is zero. Let us take this thrice iterated integral  $\int_{t_1}^{t_2} t dz(t) \int_{t_1}^{t} dz(s) \int_{t_1}^{s} dv$
We start breaking it down into more integrals with lower limit zero as
$\int_{t_1}^{t_2} t dz(t) \int_{t_1}^{t} dz(s) \int_{t_1}^{s} dv$   Integral(1)
$=\Big[\int_{0}^{t_2} t dz(t) \int_{t_1}^{t} dz(s) \int_{t_1}^{s} dv$   Integral(2)
$- \int_{0}^{t_1} t dz(t) \int_{t_1}^{t} dz(s) \int_{t_1}^{s} dv \Big]$    Integral(3)
So we have changed the limits of outer integral to zero by breaking it down into two integrals. Now taking the first of the two integrals (Integral2), we have
$\int_{0}^{t_2} t dz(t) \int_{t_1}^{t} dz(s) \int_{t_1}^{s} dv$
$=\Big[\int_{0}^{t_2} t dz(t) \int_{0}^{t} dz(s) \int_{t_1}^{s} dv$
$-\int_{0}^{t_2} t dz(t) \int_{0}^{t_1} dz(s) \int_{t_1}^{s} dv\Big]$
which will be further broken down into four integrals when we change the limits of inner integrals from lower limit zero as
$\Big[\int_{0}^{t_2} t dz(t) \int_{0}^{t} dz(s) \int_{t_1}^{s} dv$
$-\int_{0}^{t_2} t dz(t) \int_{0}^{t_1} dz(s) \int_{t_1}^{s} dv \Big]$
$=\Big[\int_{0}^{t_2} t dz(t) \int_{0}^{t} dz(s) \int_{0}^{s} dv$
$-\int_{0}^{t_2} t dz(t) \int_{0}^{t} dz(s) \int_{0}^{t_1} dv$
$-\int_{0}^{t_2} t dz(t) \int_{0}^{t_1} dz(s) \int_{0}^{s} dv$
$+\int_{0}^{t_2} t dz(t) \int_{0}^{t_1} dz(s) \int_{0}^{t_1} dv\Big]$

so what we have that Integral(2) can be written as
$\int_{0}^{t_2} t dz(t) \int_{t_1}^{t} dz(s) \int_{t_1}^{s} dv$
$=\Big[\int_{0}^{t_2} t dz(t) \int_{0}^{t} dz(s) \int_{0}^{s} dv$
$-\int_{0}^{t_2} t dz(t) \int_{0}^{t} dz(s) \int_{0}^{t_1} dv)$
$-\int_{0}^{t_2} t dz(t) \int_{0}^{t_1} dz(s) \int_{0}^{s} dv$
$+\int_{0}^{t_2} t dz(t) \int_{0}^{t_1} dz(s) \int_{0}^{t_1} dv\Big]$

Similarly Integral(3) can be written as
$\int_{0}^{t_1} t dz(t) \int_{t_1}^{t} dz(s) \int_{t_1}^{s} dv$
$=\Big[\int_{0}^{t_1} t dz(t) \int_{0}^{t} dz(s) \int_{0}^{s} dv$
$-\int_{0}^{t_1} t dz(t) \int_{0}^{t} dz(s) \int_{0}^{t_1} dv$
$-\int_{0}^{t_1} t dz(t) \int_{0}^{t_1} dz(s) \int_{0}^{s} dv$
$+\int_{0}^{t_1} t dz(t) \int_{0}^{t_1} dz(s) \int_{0}^{t_1} dv\Big]$

So our main integral(1) becomes which is sum(or difference) of integral (2) and integral(3)
$\int_{t_1}^{t_2} t dz(t) \int_{t_1}^{t} dz(s) \int_{t_1}^{s} dv$
$=\Big[\int_{0}^{t_2} t dz(t) \int_{0}^{t} dz(s) \int_{0}^{s} dv$
$-\int_{0}^{t_2} t dz(t) \int_{0}^{t} dz(s) \int_{0}^{t_1} dv$
$-\int_{0}^{t_2} t dz(t) \int_{0}^{t_1} dz(s) \int_{0}^{s} dv$
$+\int_{0}^{t_2} t dz(t) \int_{0}^{t_1} dz(s) \int_{0}^{t_1} dv$
$-\int_{0}^{t_1} t dz(t) \int_{0}^{t} dz(s) \int_{0}^{s} dv$
$+\int_{0}^{t_1} t dz(t) \int_{0}^{t} dz(s) \int_{0}^{t_1} dv$
$+\int_{0}^{t_1} t dz(t) \int_{0}^{t_1} dz(s) \int_{0}^{s} dv$
$-\int_{0}^{t_1} t dz(t) \int_{0}^{t_1} dz(s) \int_{0}^{t_1} dv\Big]$

Fortunately we can very easily calculate variances of integrals starting from lower limit zero and it is not difficult at all. I will calculate all these variances here one by one.
$Var[\int_{0}^{t_2} t dz(t) \int_{0}^{t} dz(s) \int_{0}^{s} dv]= {t_2}^6/18$
$-Var[\int_{0}^{t_2} t dz(t) \int_{0}^{t} dz(s) \int_{0}^{t_1} dv]= -{t_2}^4 {t_1}^2/4$
$-Var[\int_{0}^{t_2} t dz(t) \int_{0}^{t_1} dz(s) \int_{0}^{s} dv]= -{t_2}^3 {t_1}^3/9$
$+Var[\int_{0}^{t_2} t dz(t) \int_{0}^{t_1} dz(s) \int_{0}^{t_1} dv]={t_2}^3 {t_1}^3/3$
$-Var[\int_{0}^{t_1} t dz(t) \int_{0}^{t} dz(s) \int_{0}^{s} dv]= -{t_1}^6/18$
$+Var[\int_{0}^{t_1} t dz(t) \int_{0}^{t} dz(s) \int_{0}^{t_1} dv]={t_1}^6/4$
$+Var[\int_{0}^{t_1} t dz(t) \int_{0}^{t_1} dz(s) \int_{0}^{s} dv]={t_1}^6/9$
$-Var[\int_{0}^{t_1} t dz(t) \int_{0}^{t_1} dz(s) \int_{0}^{t_1} dv]=- {t_1}^6/3$
So summing up the variances and re-arranging, we get the total variance as
$({t_2}^6-{t_1}^6)/18 -({t_2}^4 {t_1}^2-{t_1}^6)/4-({t_2}^3 {t_1}^3-{t_1}^6)/9+({t_2}^3 {t_1}^3-{t_1}^6)/3$
or we have
$Var[\int_{t_1}^{t_2} t dz(t) \int_{t_1}^{t} dz(s) \int_{t_1}^{s} dv]$
$=({t_2}^6-{t_1}^6)/18 -({t_2}^4 {t_1}^2-{t_1}^6)/4-({t_2}^3 {t_1}^3-{t_1}^6)/9+({t_2}^3 {t_1}^3-{t_1}^6)/3$

and we can represent our integral as
$\int_{t_1}^{t_2} t dz(t) \int_{t_1}^{t} dz(s) \int_{t_1}^{s} dv$
$=\sqrt{(({t_2}^6-{t_1}^6)/18 -({t_2}^4 {t_1}^2-{t_1}^6)/4-({t_2}^3 {t_1}^3-{t_1}^6)/9+({t_2}^3 {t_1}^3-{t_1}^6)/3)} \frac{H_2(N)}{\sqrt{2!}}$
where $H_2(N)=N^2-1$ and N is a standard Gaussian.

Please pardon any inadvertent mistakes. I think I have done the calculations right.
.
.
Very sorry friends for this error. I just calculated the variances on diagonal right in the above copied post but did not include co-variances between terms. formula for variance of sum of two terms goes by $Var(a+b)=Var(a) + Var(b) +2 COV(a,b)$ and $Var(a-b)=Var(a) + Var(b) -2 COV(a,b)$ and this is the formula that has to be applied to all the terms by adding all individual variances of each term and then accounting for co-variances between each of the pair of two terms. I will fix this error when I come back to transition probabilities based stochastic integrals.
Of course, I knew the variance formulas but for some reason I was totally blank about them when I wrote the above copied post.
My work on the new project is going well and I hope to come back with results in another two to three days.