Serving the Quantitative Finance Community

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

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

Let us suppose we want to solve a general SDE given as
$dX(t)=\mu(X) dt + \sigma(X) dz(t)$          Equation (1)
The fokker planck equation of this SDE is given as
$\frac{\partial p(X,t)}{\partial t} = -\frac{\partial [\mu(X) p(X,t)]}{\partial X} + .5 \frac{\partial^2 [{\sigma(X)}^2 p(X,t)]}{\partial X^2}$ Equation(2)

First a very brief primer on derivatives of normal and derivatives of the SDE variable.
I denote the standard normal density as p(Z)
We want to convert derivatives of the SDE density of X(t) into derivatives of standard normal density. Here are the main equations
$p(Z)=\frac{1}{\sqrt{2 \pi}} exp[-.5 Z^2]$ Equation(3)
derivatives of standard normal density are given in terms of hermite polynomials.
$\frac{\partial p(Z)}{\partial Z}=-Z p(Z)$   Equation(4)
$\frac{\partial^2 p(Z)}{\partial Z^2}=(Z^2-1) p(Z)$  Equation (5)
now we want to convert the density of X(t) and its derivatives in terms of density of Z in which we also use above equations (4-5). In the density of X, each value of X(t) is associated with a particular value of Z through constant CDF points on corresponding densities. So we calculate the density and derivatives of X(t) in terms of Z.
$p(X)=p(Z) |\frac{\partial Z}{\partial X}|$  Equation(6)This follows from the standard change of variables in a density.
$\frac{\partial p(X)}{\partial X}=\frac{\partial p(Z)}{\partial Z} {(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2} =Z p(Z){(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2}$  Equation(7)
$\frac{\partial^2 p(X)}{\partial X^2}=\frac{\partial^2 p(Z)}{\partial Z^2} {(\frac{\partial Z}{\partial X})}^3+3 \frac{\partial p(Z)}{\partial Z} \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2}+p(Z) \frac{\partial^3 Z}{\partial X^3}$
$=(Z^2-1) p(Z) {(\frac{\partial Z}{\partial X})}^3 -3 Z p(Z) \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} +p(Z) \frac{\partial^3 Z}{\partial X^3}$ Equation(8)

$\frac{\partial p(X)}{\partial t}=\frac{\partial [p(Z) \frac{\partial Z}{\partial X}]}{\partial t}= \frac{\partial p(Z)}{\partial Z} {(\frac{\partial Z}{\partial X})}^2 \frac{\partial X}{\partial t}+ p(Z) \frac{\partial^2 Z}{\partial X^2} \frac{\partial X}{\partial t}$
$=Zp(Z) {(\frac{\partial Z}{\partial X})}^2 \frac{\partial X}{\partial t}+ p(Z) \frac{\partial^2 Z}{\partial X^2} \frac{\partial X}{\partial t}$   Equation(9)

Now I write the fokker planck equation after applying the derivatives as
$\frac{\partial p(X,t)}{\partial t} = -\mu(X) \frac{\partial p(X,t)}{\partial X} - \frac{\partial \mu(X)}{\partial X} p(X,t)$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) p(X,t)$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} \frac{\partial p(X,t)}{\partial X}+ .5 {\sigma(X)}^2 \frac{\partial^2 p(X,t)}{\partial X^2}$  Equation(10)

Now we write the density of X(t) in terms of density of standard normal Z by substituting from above equations (6-9) in Fokker-Planck Equation (10)

$Zp(Z) {(\frac{\partial Z}{\partial X})}^2 \frac{\partial X}{\partial t}+ p(Z) \frac{\partial^2 Z}{\partial X^2} \frac{\partial X}{\partial t}$
$= -\mu(X) (Z p(Z){(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2}) - \frac{\partial \mu(X)}{\partial X} p(Z) |\frac{\partial Z}{\partial X}|$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) p(Z) |\frac{\partial Z}{\partial X}|$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} (Z p(Z){(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2})$
$+ .5 {\sigma(X)}^2 ((Z^2-1) p(Z) {(\frac{\partial Z}{\partial X})}^3 -3 Z p(Z) \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} +p(Z) \frac{\partial^3 Z}{\partial X^3})$  Equation(11)

We see that $p(Z)$ which is a static density is common throughout the equations. We eliminate it and write the new equation as

$Z{(\frac{\partial Z}{\partial X})}^2 \frac{\partial X}{\partial t}+ \frac{\partial^2 Z}{\partial X^2} \frac{\partial X}{\partial t}$
$=- \mu(X) (Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2}) - \frac{\partial \mu(X)}{\partial X} |\frac{\partial Z}{\partial X}|$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) |\frac{\partial Z}{\partial X}|$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} (Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2})$
$+ .5 {\sigma(X)}^2 ((Z^2-1) {(\frac{\partial Z}{\partial X})}^3 -3 Z \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} + \frac{\partial^3 Z}{\partial X^3})$ Equation(12)

We multiply both sides of equations with dt and rearrange the above equation(12) to write

$dX(Z{(\frac{\partial Z}{\partial X})}^2 + \frac{\partial^2 Z}{\partial X^2} )$
$= -\mu(X) (Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2}) dt$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} (Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2}) dt$
$+ .5 {\sigma(X)}^2 ((Z^2-1) {(\frac{\partial Z}{\partial X})}^3 -3 Z \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} + \frac{\partial^3 Z}{\partial X^3}) dt$
$- \frac{\partial \mu(X)}{\partial X} |\frac{\partial Z}{\partial X}| dt$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) |\frac{\partial Z}{\partial X}| dt$  Equation(13)

We multiply both sides of the above equation 13 with ${(\frac{\partial X}{\partial Z})}^3$ to get

$dX(Z (\frac{\partial X}{\partial Z}) + {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^2 Z}{\partial X^2} )$
$= -\mu(X) (Z (\frac{\partial X}{\partial Z})+ {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^2 Z}{\partial X^2}) dt$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} (Z (\frac{\partial X}{\partial Z})+ {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^2 Z}{\partial X^2}) dt$
$- .5 {\sigma(X)}^2 dt$
$+ .5 {\sigma(X)}^2 (Z^2 -3 Z {(\frac{\partial X}{\partial Z})}^2 \frac{\partial^2 Z}{\partial X^2} +{(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3}) dt$
$- \frac{\partial \mu(X)}{\partial X} {(\frac{\partial X}{\partial Z})}^2 dt$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) {(\frac{\partial X}{\partial Z})}^2 dt$  Equation(14)

We use the relationship  ${(\frac{\partial X}{\partial Z})}^3 \frac{\partial^2 Z}{\partial X^2} =- \frac{\partial^2 X}{\partial Z^2}$ in above equation (14) to get
the following equation

$dX(Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2} )$
$=- \mu(X) (Z (\frac{\partial X}{\partial Z})- \frac{\partial^2 X}{\partial Z^2}) dt$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} (Z (\frac{\partial X}{\partial Z})- \frac{\partial^2 X}{\partial Z^2}) dt$
$- .5 {\sigma(X)}^2 dt$
$+ .5 {\sigma(X)}^2 (Z^2 +3 Z {(\frac{\partial Z}{\partial X})} \frac{\partial^2 X}{\partial Z^2} +{(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3}) dt$
$- \frac{\partial \mu(X)}{\partial X} {(\frac{\partial X}{\partial Z})}^2 dt$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) {(\frac{\partial X}{\partial Z})}^2 dt$  Equation(15)

The first four terms of the equation can be solved to give a differential of following squared form

$dX(Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2} )$
$= -\mu(X) (Z (\frac{\partial X}{\partial Z})- \frac{\partial^2 X}{\partial Z^2}) dt$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} (Z (\frac{\partial X}{\partial Z})- \frac{\partial^2 X}{\partial Z^2}) dt$
$- .5 {\sigma(X)}^2 dt$
$=.5 d\Big[\int_{t_0}^{t_1} dw - \int_{t_0}^{t_1} \mu(X) ds - \int_{t_0}^{t_1} 2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} ds$
$+ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]^2$  Equation(16)

I want to explain that the above term (equation 16)
$.5 d\Big[\int_{t_0}^{t_1} dX - \int_{t_0}^{t_1} \mu(X) ds - \int_{t_0}^{t_1} 2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} ds$
$+ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]^2$
equals
$+.5 \Big[( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]^2$ at time $t_0$ since the first three integral terms vanish at start time.
and equation 16 equals at time $t_1$
$.5 \Big[X(t_1)-X(t_0) - \int_{t_0}^{t_1} \mu(X) ds - \int_{t_0}^{t_1} 2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} ds$
$+ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]^2$

So we can write the above equation 15 after substituting equation 16 in it as
$.5 d\Big[\int_{t_0}^{t_1} dX - \int_{t_0}^{t_1} \mu(X) ds - \int_{t_0}^{t_1} 2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} ds$
$+ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]^2$  Equation(16)
$=+ .5 {\sigma(X)}^2 (Z^2 +3 Z {(\frac{\partial Z}{\partial X})} \frac{\partial^2 X}{\partial Z^2} +{(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3}) dt$
$- \frac{\partial \mu(X)}{\partial X} {(\frac{\partial X}{\partial Z})}^2 dt$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) {(\frac{\partial X}{\partial Z})}^2 dt$  Equation(17)

We can write the above equation after noting that some terms are getting constant multiplier coefficients not immediately known. So we re-write equation (17) as

$.5 \Big[X(t_1)-X(t_0) - \int_{t_0}^{t_1} \mu(X) ds - \int_{t_0}^{t_1} 2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} ds$
$+ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]^2$
$= +{\Big[ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]}^2$
$+ .5 {\sigma(X)}^2 (Z^2 +3C_2 Z {(\frac{\partial Z}{\partial X})} \frac{\partial^2 X}{\partial Z^2} +C_4 {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3}) dt$
$- C_3\frac{\partial \mu(X)}{\partial X} {(\frac{\partial X}{\partial Z})}^2 dt$
$+C_5({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) {(\frac{\partial X}{\partial Z})}^2 dt$  Equation(18)

After taking square root on both sides and re-arranging, we get the evolution equation for X(t) as

$X(t_1)=X(t_0) + \int_{t_0}^{t_1} \mu(X) ds + \int_{t_0}^{t_1} 2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} ds$
$- Z (\frac{\partial X}{\partial Z}) + \frac{\partial^2 X}{\partial Z^2}$
$+\sqrt{\Big[ +{\big[ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\big]}^2 + {\sigma(X)}^2 \big[ Z^2 +3C_2 Z {(\frac{\partial Z}{\partial X})} \frac{\partial^2 X}{\partial Z^2} +C_4 {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3} \big] dt - C_3\frac{\partial \mu(X)}{\partial X} {(\frac{\partial X}{\partial Z})}^2 dt+C_5({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) {(\frac{\partial X}{\partial Z})}^2 dt \Big]}$
Equation(19)

Please note that all the terms on RHS appear under square-root. The last matlab program I posted has two terms outside the square-root and are being linearly added and that is an error. I will post a new program with this solution in a day and it works far better closer to zero.
Here I mention that Lamperti/Bessel form SDE is given as
$dX(t)=\mu(X) dt + \sigma dz(t)$  Equation(20)
and it has a solution given as

$X(t_1)=X(t_0) + \int_{t_0}^{t_1} \mu(X) ds - Z (\frac{\partial X}{\partial Z}) + \frac{\partial^2 X}{\partial Z^2}$
$+\sqrt{\Big[ +{\big[ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\big]}^2+ {\sigma(X)}^2 \big[Z^2 +3C_2 Z {(\frac{\partial Z}{\partial X})} \frac{\partial^2 X}{\partial Z^2} +C_4 {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3}\big] dt- C_3\frac{\partial \mu(X)}{\partial X} {(\frac{\partial X}{\partial Z})}^2) dt \Big]}$ Equation (21)

Please pardon any error with signs. I quickly wrote the post but I will carefully check it again and fix any errors in a day or two.
Last edited by Amin on August 9th, 2020, 7:18 am, edited 4 times in total.

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

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

I have made some new graphs with the new evolution equation for Fokker-Planck as I wrote in the previous post. I am posting the graphs. I have also given the time step size on graphs so that friends could reproduce them on their computers. My chosen constants are still not perfectly optimal. For many difficult equations we have to take a very small step size for solution of FPE of many difficult SDEs. I have also removed unnecessary computations from the main sub-function and that helps decrease the run-time to one third of previous time required with the same step size.
In the graphs below I have simulated fokker-planck equation of density of following SDE with various values of parameters.

$dx(t)=\kappa (\theta-x(t)) + \sigma x(t)^{\gamma} dz(t)$
with $x(0)=x_0$
You have to interpret the parameters given in title of each graph with reference to above SDE. I have solved the fokker-Planck equation related to above SDE in Lamperti/Bessel form with the last formula given at the end of previous post. On the graph title, M means SDE density mean by new method and TM means true SDE mean from analytical formula.

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

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

Here I am posting my new program for calculation of densities of SDEs from Fokker-Planck equation. I am also posting a sub-function that I have altered from previous version. You should be able to reproduce graphs posted in previous post using this program.
function [TrueMean,ItoHermiteMean,y_w,fy_w,wnStart] = FPERevisited044Wilmott16_2ndD08()

%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.
%In this version, the mean correction has been disabled. If you are
%simulating standard mean-reverting SV type SDEs, please enable the mean
%correction by uncommenting the appropriate line in the body of code below.

dt=.125/16/2/2/2;   % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=128*8/4*2*2*2;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
T=Tt*dt;
OrderA=4;  %
OrderM=4;  %
%dtM=.125/8;%Monte carlo time interval size dtM.
%TtM=8*8;%Monte carlo number of simulation intervals.
dtM=dt*4*4;
TtM=Tt/4/4;

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=.160;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.65;%50;   % volatility power.
kappa=1.0;%.950;   %mean reversion parameter.
theta=.04;%mean reversion target
sigma0=.550;%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
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);

for nn=2:Nn-1
ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1)-normcdf(.5*Z(nn)+.5*Z(nn-1),0,1);
%Above calculate probability mass in each probability subdivision.
end

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

sigma11(1:OrderA+1)=0;
mu11(1:OrderA+1)=0;
mu22(1:OrderA+1)=0;
sigma22(1:OrderA+1)=0;
% index 1 correponds to zero level since matlab indexing starts at one.
sigma11(1)=1;
mu11(1)=1;
mu22(1)=1;
sigma22(1)=1;

for k=1:(OrderA+1)
if sigma0~=0
sigma11(k)=sigma0^(k-1);
end
if mu1 ~= 0
mu11(k)=mu1^(k-1);
end
if mu2 ~= 0
mu22(k)=mu2^(k-1);
end
if sigma0~=0
sigma22(k)=sigma0^(2*(k-1));
end
end
%Ft(1:TtM+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0; %General time powers on hermite polynomials
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers on coefficients of hermite polynomials.
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers for bessel transformed coordinates.

%YCoeff0 and YCoeff are coefficents for original coordinates monte carlo.
%YqCoeff0 and YqCoeff are bessel/lamperti version monte carlo.

YCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
YqCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
%Pre-compute the time and power exponent values in small multi-dimensional arrays
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma); %expand y^alpha where alpha=1;
YqCoeff = ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);%expand y^alpha1 where alpha1=(1-gamma)
YqCoeff=YqCoeff/(1-gamma); %Transformed coordinates coefficients have to be
%further divided by (1-gamma)

for k = 0 : (OrderA)
for m = 0:k
l4 = k - m + 1;
for n = 0 : m
l3 = m - n + 1;
for j = 0:n
l2 = n - j + 1;
l1 = j + 1;
%Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
- (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
- (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));

YCoeff0(l1,l2,l3,l4) =YCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
YqCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
end
end
end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wnStart=1;%
d2wdZ2(1:Nn)=0;
d2wdZ2A(1:Nn)=0;

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

for tt=1:Tt

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

[wMu0dt,dwMu0dtdw,c1] = CalculateDriftAndVolA4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);

%wPrev=w;
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.

%yyMedian=.5*yy(NnMidh)+.5*yy(NnMidl);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[dwdZ,d2wdZ2A,d3wdZ3A] = First3Derivatives4thOrderEqSpaced(wnStart,Nn,dNn,w,Z);

d2wdZ2(wnStart:Nn)=0.0;
d3wdZ3(wnStart:Nn)=0.0;
if (tt>32)

d2wdZ2(wnStart:Nn)=d2wdZ2A(wnStart:Nn);
d2wdZ2(1:10)=d2wdZ2(1:10).*exp(-12*kappa^2*(1-gamma).*dNn*(10-(1:10)));
%Change 1:7/25/2020: I have very heavily weighted the second derivative to
%the left of density peak with an exponential and it allows far larger
%range of parameters with stable execution without sacrificing much
%accuracy. Of course, it is just ad hoc.

d3wdZ3(wnStart:Nn)=d3wdZ3A(wnStart:Nn);
end

%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)
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=1.1/pi/8;
%C33=.9/pi/16;%
%C44=1.1/pi/8;

%  C22=1.1/pi/8*1.5;
%  C33=.5*.5*.707/pi/16;%.9/pi/16;%
%  C44=1.1/pi/8*1.5;

C22=.28125/pi;
C33=.01104/pi;
C44=.28125/pi;

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

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

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)+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. There is instability in the
% SDE representation equation roughly between -.5 SD to +.5 SD in normal terms.
%Instability correction Algorithm takes two ends that are here given as
%lower end=w(NnMidl-2) and upper end=w(NnMidh+3).

%Change 8/3/2020
%I have switched to matlab spline function interpolation by taking five
%points on both sides of instability. This works far better and you can
%now use nuemrical 2nd derivative without any modification for a large
%range of parameters. It truly improves the fit to true density and usually
%works great.

Zi=[Z(NnMidl-6),Z(NnMidl-5),Z(NnMidl-4),Z(NnMidl-3),Z(NnMidl-2),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7)];
wi=[w(NnMidl-6),w(NnMidl-5),w(NnMidl-4),w(NnMidl-3),w(NnMidl-2),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidh+7)];
Zq=[Z(NnMidl-1),Z(NnMidl),Z(NnMidh),Z(NnMidh+1),Z(NnMidh+2)];
wq=spline(Zi,wi,Zq);
w(NnMidl-1:NnMidh+2)=wq(1: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 last line in the block.

end

yyMedian=.5*yy(NnMidh)+.5*yy(NnMidl);
y_w(1:Nn)=0;
y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
Dfy_w(wnStart:Nn)=0;
Dfw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
Dfy_w(nn) = (y_w(nn + 1) - y_w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
Dfw(nn) = (w(nn + 1) - w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
fy_w(1:Nn)=0;
fw(1:Nn)=0;
for nn = wnStart:Nn-1
fy_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
fw(nn) = (normpdf(Z(nn),0, 1))/abs(Dfw(nn));
end

toc

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

rng(29079137, 'twister')
paths=100000;
YY(1:paths)=x0;  %Original process monte carlo.
Random1(1:paths)=0;
YYMean(1:TtM)=0;
for tt=1:TtM
Random1=randn(size(Random1));
HermiteP1(1,1:paths)=1;
HermiteP1(2,1:paths)=Random1(1:paths);
HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;

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

%Uncomment for fourth order monte carlo

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

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

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

MaxCutOff=30;
NoOfBins=round(300*gamma^2*8*sigma0/sqrt(yyMedian)/(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),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');

title(sprintf('x0 = %.4f,theta=%.2f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.5f,TM=%.5f', 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.');

end
Here is the new sub-function that you will also need.
function [wMu0dt,dwMu0dtdw,c1] = CalculateDriftAndVolA4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt)

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

Fp2=Fp1/(1-gamma);

wMu0dt(wnStart:Nn)=(YqCoeff0(1,1,2,1).*yy(wnStart:Nn).^Fp1(1,1,2,1)+ ...
YqCoeff0(1,2,1,1).*yy(wnStart:Nn).^Fp1(1,2,1,1)+ ...
YqCoeff0(2,1,1,1).*yy(wnStart:Nn).^Fp1(2,1,1,1))*dt + ...
(YqCoeff0(1,1,3,1).*yy(wnStart:Nn).^Fp1(1,1,3,1)+ ...
YqCoeff0(1,2,2,1).*yy(wnStart:Nn).^Fp1(1,2,2,1)+ ...
YqCoeff0(2,1,2,1).*yy(wnStart:Nn).^Fp1(2,1,2,1)+ ...
YqCoeff0(1,3,1,1).*yy(wnStart:Nn).^Fp1(1,3,1,1)+ ...
YqCoeff0(2,2,1,1).*yy(wnStart:Nn).^Fp1(2,2,1,1)+ ...
YqCoeff0(3,1,1,1).*yy(wnStart:Nn).^Fp1(3,1,1,1))*dt^2;% + ...
%     (YqCoeff0(1,1,4,1).*yy(wnStart:Nn).^Fp1(1,1,4,1)+ ...
%     YqCoeff0(1,2,3,1).*yy(wnStart:Nn).^Fp1(1,2,3,1)+ ...
%     YqCoeff0(2,1,3,1).*yy(wnStart:Nn).^Fp1(2,1,3,1)+ ...
%     YqCoeff0(1,3,2,1).*yy(wnStart:Nn).^Fp1(1,3,2,1)+ ...
%     YqCoeff0(2,2,2,1).*yy(wnStart:Nn).^Fp1(2,2,2,1)+ ...
%     YqCoeff0(3,1,2,1).*yy(wnStart:Nn).^Fp1(3,1,2,1)+ ...
%     YqCoeff0(1,4,1,1).*yy(wnStart:Nn).^Fp1(1,4,1,1)+ ...
%     YqCoeff0(2,3,1,1).*yy(wnStart:Nn).^Fp1(2,3,1,1)+ ...
%     YqCoeff0(3,2,1,1).*yy(wnStart:Nn).^Fp1(3,2,1,1)+ ...
%     YqCoeff0(4,1,1,1).*yy(wnStart:Nn).^Fp1(4,1,1,1))*dt^3;%+ ...
%     (YqCoeff0(1,1,5,1).*yy(wnStart:Nn).^Fp1(1,1,5,1)+ ...
%      YqCoeff0(1,2,4,1).*yy(wnStart:Nn).^Fp1(1,2,4,1)+ ...
%      YqCoeff0(2,1,4,1).*yy(wnStart:Nn).^Fp1(2,1,4,1)+ ...
%      YqCoeff0(2,1,1,1).*yy(wnStart:Nn).^Fp1(2,1,1,1)+ ...
%       YqCoeff0(2,2,3,1).*yy(wnStart:Nn).^Fp1(2,2,3,1)+ ...
%       YqCoeff0(3,1,3,1).*yy(wnStart:Nn).^Fp1(3,1,3,1)+ ...
%       YqCoeff0(1,4,2,1).*yy(wnStart:Nn).^Fp1(1,4,2,1)+ ...
%       YqCoeff0(2,3,2,1).*yy(wnStart:Nn).^Fp1(2,3,2,1)+ ...
%       YqCoeff0(3,2,2,1).*yy(wnStart:Nn).^Fp1(3,2,2,1)+ ...
%       YqCoeff0(4,1,2,1).*yy(wnStart:Nn).^Fp1(4,1,2,1)+ ...
%       YqCoeff0(1,5,1,1).*yy(wnStart:Nn).^Fp1(1,5,1,1)+ ...
%       YqCoeff0(2,4,1,1).*yy(wnStart:Nn).^Fp1(2,4,1,1)+ ...
%       YqCoeff0(3,3,1,1).*yy(wnStart:Nn).^Fp1(3,3,1,1)+  ...
%       YqCoeff0(4,2,1,1).*yy(wnStart:Nn).^Fp1(4,2,1,1)+ ...
%       YqCoeff0(5,1,1,1).*yy(wnStart:Nn).^Fp1(5,1,1,1))*dt^4;

dwMu0dtdw(wnStart:Nn)=(YqCoeff0(1,1,2,1).*Fp1(1,1,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,1,2,1))+ ...
YqCoeff0(1,2,1,1).*Fp1(1,2,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,2,1,1))+ ...
YqCoeff0(2,1,1,1).*Fp1(2,1,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,1,1,1)))*dt;% + ...
%     (YqCoeff0(1,1,3,1).*Fp1(1,1,3,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,1,3,1))+ ...
%     YqCoeff0(1,2,2,1).*Fp1(1,2,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,2,2,1))+ ...
%     YqCoeff0(2,1,2,1).*Fp1(2,1,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,1,2,1))+ ...
%     YqCoeff0(1,3,1,1).*Fp1(1,3,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,3,1,1))+ ...
%     YqCoeff0(2,2,1,1).*Fp1(2,2,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,2,1,1))+ ...
%     YqCoeff0(3,1,1,1).*Fp1(3,1,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(3,1,1,1)))*dt^2;% + ...
%     (YqCoeff0(1,1,4,1).*Fp1(1,1,4,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,1,4,1))+ ...
%     YqCoeff0(1,2,3,1).*Fp1(1,2,3,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,2,3,1))+ ...
%     YqCoeff0(2,1,3,1).*Fp1(2,1,3,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,1,3,1))+ ...
%     YqCoeff0(1,3,2,1).*Fp1(1,3,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,3,2,1))+ ...
%     YqCoeff0(2,2,2,1).*Fp1(2,2,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,2,2,1))+ ...
%     YqCoeff0(3,1,2,1).*Fp1(3,1,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(3,1,2,1))+ ...
%     YqCoeff0(1,4,1,1).*Fp1(1,4,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,4,1,1))+ ...
%     YqCoeff0(2,3,1,1).*Fp1(2,3,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,3,1,1))+ ...
%     YqCoeff0(3,2,1,1).*Fp1(3,2,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(3,2,1,1))+ ...
%     YqCoeff0(4,1,1,1).*Fp1(4,1,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(4,1,1,1)))*dt^3;

c1(wnStart:Nn)=YqCoeff0(1,1,1,2).*yy(wnStart:Nn).^Fp1(1,1,1,2).*sqrt(dt)+ ...
(YqCoeff0(1,1,2,2).*yy(wnStart:Nn).^Fp1(1,1,2,2)+YqCoeff0(1,2,1,2).*yy(wnStart:Nn).^Fp1(1,2,1,2)+ ...
YqCoeff0(2,1,1,2).*yy(wnStart:Nn).^Fp1(2,1,1,2)).*dt^1.5;%+ ...
%(YqCoeff0(1,1,3,2).*yy(wnStart:Nn).^Fp1(1,1,3,2)+YqCoeff0(1,2,2,2).*yy(wnStart:Nn).^Fp1(1,2,2,2)+ ...
%YqCoeff0(2,1,2,2).*yy(wnStart:Nn).^Fp1(2,1,2,2)+YqCoeff0(1,3,1,2).*yy(wnStart:Nn).^Fp1(1,3,1,2)+ ...
%YqCoeff0(2,2,1,2).*yy(wnStart:Nn).^Fp1(2,2,1,2)+YqCoeff0(3,1,1,2).*yy(wnStart:Nn).^Fp1(3,1,1,2)).*dt^2.5;

end

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

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

Since step size in the matlab program copied to previous post is too small, I want to suggest something. Here is the relevant block of code in previous program

d2wdZ2(wnStart:Nn)=0.0;
d3wdZ3(wnStart:Nn)=0.0;
if (tt>32)

d2wdZ2(wnStart:Nn)=d2wdZ2A(wnStart:Nn);
d2wdZ2(1:10)=d2wdZ2(1:10).*exp(-12*kappa^2*(1-gamma).*dNn*(10-(1:10)));
%Change 1:7/25/2020: I have very heavily weighted the second derivative to
%the left of density peak with an exponential and it allows far larger
%range of parameters with stable execution without sacrificing much
%accuracy. Of course, it is just ad hoc.
d3wdZ3(wnStart:Nn)=d3wdZ3A(wnStart:Nn);
end

if in the above code, you replace  "if (tt>32) " with  "if (tt>64) " , you can take a larger step size in many cases. May be WENO or ENO type non-oscillatory derivatives calculation could help a lot.

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

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

My this post concerns my two previous posts. Here is post # 973 made on Mon Aug 03, 2020 5:01 pm
.
I was trying to post on some linkedin groups asking friends to download the program in previous posts but my posts do not appear on some groups even after I see the notification that my post was successful. It seems that mind control and allied crooks have given money to their horses who want to present my research as their own. Why are they stopping my research from reaching other people?
.
And here is my post # 974 written on Mon Aug 03, 2020 6:29 pm
.
Lot of friends in several countries can recall how crooks in US defence and their allied apparatus were telling my linkedin connections when I did my monte carlo simulation research and approached my connections about my research that I was a very bad guy and "horses" of US defence crooks would present the same (my) research to them. Well, these crooks never grow up and they want to do the same when I have presented solution to fokker-planck equation.
I will try to post latex equations for the solution of Fokker-planck equation later today. But in the meantime, I want to request friends to protest against my persecution. Here I want to share some threads with friends especially my application to human rights watch for help.
My Letter To HRW About My Human Rights Abuse
Muslims, mind control, Muslim-Jewish relationship and American defense
My Open Letter to Prime Minister of Pakistan, Imran Khan, to Intervene In Order to Protect the Interests of Pakistan
.
Friends know that I have been posting my research and code for a long time even while my research would not be perfectly developed and many things in my research would be unknown. I understand that once I post something on internet, it is in public domain and friends are free to do their advanced research about it. I want to say that I encourage friends to freely do their own research even when it borrows something from my research. Only thing I expect that they would give me credit for what I have done.
But what I do not understand why would mind control and related US Army apparatus try to stop my research from reaching other people. I expect that people understand that this is very wrong. Obviously because they do not want my independent research to reach other people before somebody else reaches out. I hope that people understand this was not fair.
I apologize to friends for my remarks and want to tell them that I was already on very bad and evil mind control neurotransmitters that had been given to me through tained food otherwise I would really have remained balanced. I might have remained balanced despite mind control drugs but I became very upset when I knew that mind control agents were blocking my LinkedIn posts to groups and messages. What I am trying to say might be very hard for many friends who have not gone through mind control themselves. But I again apologize to friends and want to tell them that I understand that once my research is in public domain, they are free to do their own research about it but please ask mind control forces in US Army to end their attempts at stopping my research from reaching other people. Since past one week, mind control has dramatically increased and I lost my sanity multiple times and I am relatively sane now only since I have been able to avoid drugged food to some extent. But it is very difficult to avoid drugged food for larger periods of time since my family is a party with US Army and they detain me with psychiatric facilities for several months calling me schizophrenic once I try to avoid their drugged food. Please protest against this racial/religious injustice and torture on me. I am not a practicing Muslim but unfortunately I cannot change my birth history.
In the meantime, I want to request friends to protest against my persecution. Here I again want to share some threads with friends especially my application to human rights watch for help so friends have more information about my persecution and the persecution of thousands of other totally people like me whose only sin was to have some special neurotransmitter in some part of their brain and US Army decided to retard them.
My Letter to HRW About My Human Rights Abuse
Muslims, mind control, Muslim-Jewish relationship and American defense
My Open Letter to Prime Minister of Pakistan, Imran Khan, to Intervene In Order to Protect the Interests of Pakistan

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

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

I am re-writing the post 979 above after making some additions.

Since step size in the matlab program copied to previous post is too small, I want to suggest something. Here is the relevant block of code in previous program

d2wdZ2(wnStart:Nn)=0.0;
d3wdZ3(wnStart:Nn)=0.0;
if (tt>32)

d2wdZ2(wnStart:Nn)=d2wdZ2A(wnStart:Nn);
d2wdZ2(1:10)=d2wdZ2(1:10).*exp(-12*kappa^2*(1-gamma).*dNn*(10-(1:10)));
%Change 1:7/25/2020: I have very heavily weighted the second derivative to
%the left of density peak with an exponential and it allows far larger
%range of parameters with stable execution without sacrificing much
%accuracy. Of course, it is just ad hoc.
d3wdZ3(wnStart:Nn)=d3wdZ3A(wnStart:Nn);
end

if in the above code,
1) you replace  "if (tt>32) " with  "if (tt>64) " ,
and
2) you replace
"d2wdZ2(1:10)=d2wdZ2(1:10).*exp(-12*kappa^2*(1-gamma).*dNn*(10-(1:10)));"
with
"d2wdZ2(1:15)=d2wdZ2(1:15).*exp(-12*kappa^2*(1-gamma).*dNn*(15-(1:15)));"
you can take a larger step size in many cases. May be WENO or ENO type non-oscillatory derivatives calculation could help a lot.

The reason behind first replacement is that 2nd derivative is very small at the start and therefore becomes oscillatory at start and we just add it in calculationsat a later stage when it has stabilized and has become significant.
And the meanings behind 2nd replacement is that we make 2nd derivative zero for three normal standard deviations (-2 SD to -5 SD) instead of (-3 SD to -5 SD). This is again because derivative does not contribute much to accuracy in this area but becomes oscillatory (-2 SD to -5 SD) close to zero.

if you follow this post, you can take up to eight times longer step size in many cases and four time step-size in a lot of other cases. In very extreme scenarios you may only be able to take just two times longer time step size.
I believe that step size restriction is only due to oscillations and instability caused by second derivative. If you can control these oscillations, you could easily take a time step size that is 16 or 32 times larger with good accuracy.

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

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

I am a member of about 50 groups on linkedin. I had posted on a large number of groups but my posts did not appear on most linkedin groups (appeared on < 10 groups.) I thought they were stopping me from posting on groups because of their devices on my computer. So I went out to a netcafe and posted on remaining groups from the netcafe. But as it turns out, my posts still did not appear on any of those groups on which I posted from netcafe. Mind control agents might have prepared the netcafe computers already because of help from crooks in Pakistan Army so it seemed to me that I had posted on groups but my posts were still not on most internet groups. Why do the crooks in American defense have to censor my research before reaching other people? Because they do not want people to see my genuine research before somebody else somewhere presents the research to people because once people know who originally kept on working on the research for several years and took it to almost conclusion, they would have enough regard for my efforts . Please protest against the censoring of my research by US defense crooks and stopping me from reaching out to friends. Those who fashionably keep scolding Chinese, Russians and Turkish, please try to see the dark below the lamp.
I went ahead right now and tried to post on some groups after writing the above lines. I posted on three groups: Computational Physics, Mathematical Finance and Biomedical Engineering. My post actually appeared on all three groups but once I surfed out from the above groups to somewhere else and then came back after a minute to groups to see whether my post was there to realize that there was no post on all three groups. My post was already gone from all three groups.

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

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

I am posting a new version of the FPE solution program. In earlier version, I was interpolating five points in the instability region with matlab specific function "spline". In this new version I am interpolating only three points with Lagrange polynomials interpolation which is just as accurate and faster. I have removed any matlab specific function from this program. I have also removed all unnecessary things and the program is quite fast now. If SDE density does not go into zero, you can easily use the program in an industrial setting. Sometimes, you might have to be careful about step size. Here is the new program.
function [TrueMean,ItoHermiteMean,y_w,fy_w,wnStart] = FPERevisited044Wilmott16_2ndD40()

%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.
%In this version, the mean correction has been disabled. If you are
%simulating standard mean-reverting SV type SDEs, please enable the mean
%correction by uncommenting the appropriate line in the body of code below.

dt=.125/16/2/2;   % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=128*8/4*2*2;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
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*4*2;
TtM=Tt/4/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=.1;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.85;%50;   % volatility power.
kappa=1.0;%.950;   %mean reversion parameter.
theta=.1;%mean reversion target
sigma0=.950;%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
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

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

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.

end

y_w(1:Nn)=0;
y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
Dfy_w(wnStart:Nn)=0;
Dfw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
Dfy_w(nn) = (y_w(nn + 1) - y_w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
Dfw(nn) = (w(nn + 1) - w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
fy_w(1:Nn)=0;
fw(1:Nn)=0;
for nn = wnStart:Nn-1
fy_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
fw(nn) = (normpdf(Z(nn),0, 1))/abs(Dfw(nn));
end

toc

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=100000;
YY(1:paths)=x0;  %Original process monte carlo.
Random1(1:paths)=0;
YYMean(1:TtM)=0;
for tt=1:TtM
Random1=randn(size(Random1));
HermiteP1(1,1:paths)=1;
HermiteP1(2,1:paths)=Random1(1:paths);
HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;

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

%Uncomment for fourth order monte carlo

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

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

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

MaxCutOff=30;
NoOfBins=round(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),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');

title(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 are supporting functions. For Ito-hermite and monte carlo supporting functions, you will have to look at my previous posts.
1.
function [wMu0dt,dwMu0dtdw,c1] = CalculateDriftAndVolA4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt)

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

Fp2=Fp1/(1-gamma);

wMu0dt(wnStart:Nn)=(YqCoeff0(1,1,2,1).*yy(wnStart:Nn).^Fp1(1,1,2,1)+ ...
YqCoeff0(1,2,1,1).*yy(wnStart:Nn).^Fp1(1,2,1,1)+ ...
YqCoeff0(2,1,1,1).*yy(wnStart:Nn).^Fp1(2,1,1,1))*dt + ...
(YqCoeff0(1,1,3,1).*yy(wnStart:Nn).^Fp1(1,1,3,1)+ ...
YqCoeff0(1,2,2,1).*yy(wnStart:Nn).^Fp1(1,2,2,1)+ ...
YqCoeff0(2,1,2,1).*yy(wnStart:Nn).^Fp1(2,1,2,1)+ ...
YqCoeff0(1,3,1,1).*yy(wnStart:Nn).^Fp1(1,3,1,1)+ ...
YqCoeff0(2,2,1,1).*yy(wnStart:Nn).^Fp1(2,2,1,1)+ ...
YqCoeff0(3,1,1,1).*yy(wnStart:Nn).^Fp1(3,1,1,1))*dt^2;% + ...
%     (YqCoeff0(1,1,4,1).*yy(wnStart:Nn).^Fp1(1,1,4,1)+ ...
%     YqCoeff0(1,2,3,1).*yy(wnStart:Nn).^Fp1(1,2,3,1)+ ...
%     YqCoeff0(2,1,3,1).*yy(wnStart:Nn).^Fp1(2,1,3,1)+ ...
%     YqCoeff0(1,3,2,1).*yy(wnStart:Nn).^Fp1(1,3,2,1)+ ...
%     YqCoeff0(2,2,2,1).*yy(wnStart:Nn).^Fp1(2,2,2,1)+ ...
%     YqCoeff0(3,1,2,1).*yy(wnStart:Nn).^Fp1(3,1,2,1)+ ...
%     YqCoeff0(1,4,1,1).*yy(wnStart:Nn).^Fp1(1,4,1,1)+ ...
%     YqCoeff0(2,3,1,1).*yy(wnStart:Nn).^Fp1(2,3,1,1)+ ...
%     YqCoeff0(3,2,1,1).*yy(wnStart:Nn).^Fp1(3,2,1,1)+ ...
%     YqCoeff0(4,1,1,1).*yy(wnStart:Nn).^Fp1(4,1,1,1))*dt^3;%+ ...
%     (YqCoeff0(1,1,5,1).*yy(wnStart:Nn).^Fp1(1,1,5,1)+ ...
%      YqCoeff0(1,2,4,1).*yy(wnStart:Nn).^Fp1(1,2,4,1)+ ...
%      YqCoeff0(2,1,4,1).*yy(wnStart:Nn).^Fp1(2,1,4,1)+ ...
%      YqCoeff0(2,1,1,1).*yy(wnStart:Nn).^Fp1(2,1,1,1)+ ...
%       YqCoeff0(2,2,3,1).*yy(wnStart:Nn).^Fp1(2,2,3,1)+ ...
%       YqCoeff0(3,1,3,1).*yy(wnStart:Nn).^Fp1(3,1,3,1)+ ...
%       YqCoeff0(1,4,2,1).*yy(wnStart:Nn).^Fp1(1,4,2,1)+ ...
%       YqCoeff0(2,3,2,1).*yy(wnStart:Nn).^Fp1(2,3,2,1)+ ...
%       YqCoeff0(3,2,2,1).*yy(wnStart:Nn).^Fp1(3,2,2,1)+ ...
%       YqCoeff0(4,1,2,1).*yy(wnStart:Nn).^Fp1(4,1,2,1)+ ...
%       YqCoeff0(1,5,1,1).*yy(wnStart:Nn).^Fp1(1,5,1,1)+ ...
%       YqCoeff0(2,4,1,1).*yy(wnStart:Nn).^Fp1(2,4,1,1)+ ...
%       YqCoeff0(3,3,1,1).*yy(wnStart:Nn).^Fp1(3,3,1,1)+  ...
%       YqCoeff0(4,2,1,1).*yy(wnStart:Nn).^Fp1(4,2,1,1)+ ...
%       YqCoeff0(5,1,1,1).*yy(wnStart:Nn).^Fp1(5,1,1,1))*dt^4;

dwMu0dtdw(wnStart:Nn)=(YqCoeff0(1,1,2,1).*Fp1(1,1,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,1,2,1))+ ...
YqCoeff0(1,2,1,1).*Fp1(1,2,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,2,1,1))+ ...
YqCoeff0(2,1,1,1).*Fp1(2,1,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,1,1,1)))*dt;% + ...
%    (YqCoeff0(1,1,3,1).*Fp1(1,1,3,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,1,3,1))+ ...
%    YqCoeff0(1,2,2,1).*Fp1(1,2,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,2,2,1))+ ...
%    YqCoeff0(2,1,2,1).*Fp1(2,1,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,1,2,1))+ ...
%    YqCoeff0(1,3,1,1).*Fp1(1,3,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,3,1,1))+ ...
%    YqCoeff0(2,2,1,1).*Fp1(2,2,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,2,1,1))+ ...
%    YqCoeff0(3,1,1,1).*Fp1(3,1,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(3,1,1,1)))*dt^2;% + ...
%     (YqCoeff0(1,1,4,1).*Fp1(1,1,4,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,1,4,1))+ ...
%     YqCoeff0(1,2,3,1).*Fp1(1,2,3,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,2,3,1))+ ...
%     YqCoeff0(2,1,3,1).*Fp1(2,1,3,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,1,3,1))+ ...
%     YqCoeff0(1,3,2,1).*Fp1(1,3,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,3,2,1))+ ...
%     YqCoeff0(2,2,2,1).*Fp1(2,2,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,2,2,1))+ ...
%     YqCoeff0(3,1,2,1).*Fp1(3,1,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(3,1,2,1))+ ...
%     YqCoeff0(1,4,1,1).*Fp1(1,4,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,4,1,1))+ ...
%     YqCoeff0(2,3,1,1).*Fp1(2,3,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,3,1,1))+ ...
%     YqCoeff0(3,2,1,1).*Fp1(3,2,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(3,2,1,1))+ ...
%     YqCoeff0(4,1,1,1).*Fp1(4,1,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(4,1,1,1)))*dt^3;

c1(wnStart:Nn)=YqCoeff0(1,1,1,2).*yy(wnStart:Nn).^Fp1(1,1,1,2).*sqrt(dt)+ ...
(YqCoeff0(1,1,2,2).*yy(wnStart:Nn).^Fp1(1,1,2,2)+YqCoeff0(1,2,1,2).*yy(wnStart:Nn).^Fp1(1,2,1,2)+ ...
YqCoeff0(2,1,1,2).*yy(wnStart:Nn).^Fp1(2,1,1,2)).*dt^1.5;%+ ...

end
.
2.
function [dwdZ,d2wdZ2,d3wdZ3] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z)

[wS] = InterpolateOrderN6(6,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5));
[wE] = InterpolateOrderN6(6,Z(Nn)+dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));

[wS2] = InterpolateOrderN6(6,Z(wnStart)-2*dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5));
[wE2] = InterpolateOrderN6(6,Z(Nn)+2*dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));

d2wdZ2(wnStart)=(wS-2*w(wnStart)+w(wnStart+1))/(dNn.^2);
d2wdZ2(wnStart+1:Nn-1)=(w(wnStart:Nn-2)-2*w(wnStart+1:Nn-1)+1*w(wnStart+2:Nn))/(dNn.^2);
d2wdZ2(Nn)=(wE-2*w(Nn)+w(Nn-1))/(dNn.^2);

dwdZ(wnStart)=(-1*wS+1*w(wnStart+1))/(2*dNn);
dwdZ(wnStart+1:Nn-1)=(-1*w(wnStart:Nn-2)+1*w(wnStart+2:Nn))/(2*dNn);
dwdZ(Nn)=(wE-w(Nn-1))/(2*dNn);

d3wdZ3(wnStart+2:Nn-2)=(-1*w(wnStart:Nn-4)+2*w(wnStart+1:Nn-3)+0*w(wnStart+2:Nn-2)-2*w(wnStart+3:Nn-1)+w(wnStart+4:Nn))/(2*dNn^3);
d3wdZ3(wnStart+1)=(-1*wS+2*w(wnStart)+0*w(wnStart+1)-2*w(wnStart+2)+w(wnStart+3))/(2*dNn^3);
d3wdZ3(wnStart)=(-1*wS2+2*wS+0*w(wnStart)-2*w(wnStart+1)+w(wnStart+2))/(2*dNn^3);
d3wdZ3(Nn-1)=(-1*w(Nn-3)+2*w(Nn-2)+0*w(Nn-1)-2*w(Nn)+wE)/(2*dNn^3);
d3wdZ3(Nn)=(-1*w(Nn-2)+2*w(Nn-1)+0*w(Nn)-2*wE+wE2)/(2*dNn^3);

end
.
3.
function [y0] = InterpolateOrderN8(N,x0,x1,x2,x3,x4,x5,x6,x7,x8,y1,y2,y3,y4,y5,y6,y7,y8)

if(N==8)
y0=(x0-x2)*(x0-x3)*(x0-x4)*(x0-x5)*(x0-x6)*(x0-x7)*(x0-x8)/((x1-x2).*(x1-x3).*(x1-x4).*(x1-x5).*(x1-x6).*(x1-x7).*(x1-x8))*y1+ ...
(x0-x1)*(x0-x3)*(x0-x4)*(x0-x5)*(x0-x6)*(x0-x7)*(x0-x8)/((x2-x1).*(x2-x3).*(x2-x4).*(x2-x5).*(x2-x6).*(x2-x7).*(x2-x8))*y2+ ...
(x0-x1)*(x0-x2)*(x0-x4)*(x0-x5)*(x0-x6)*(x0-x7)*(x0-x8)/((x3-x1).*(x3-x2).*(x3-x4).*(x3-x5).*(x3-x6).*(x3-x7).*(x3-x8))*y3+ ...
(x0-x1)*(x0-x2)*(x0-x3)*(x0-x5)*(x0-x6)*(x0-x7)*(x0-x8)/((x4-x1).*(x4-x2).*(x4-x3).*(x4-x5).*(x4-x6).*(x4-x7).*(x4-x8))*y4+ ...
(x0-x1)*(x0-x2)*(x0-x3)*(x0-x4)*(x0-x6)*(x0-x7)*(x0-x8)/((x5-x1).*(x5-x2).*(x5-x3).*(x5-x4).*(x5-x6).*(x5-x7).*(x5-x8))*y5+ ...
(x0-x1)*(x0-x2)*(x0-x3)*(x0-x4)*(x0-x5)*(x0-x7)*(x0-x8)/((x6-x1).*(x6-x2).*(x6-x3).*(x6-x4).*(x6-x5).*(x6-x7).*(x6-x8))*y6+ ...
(x0-x1)*(x0-x2)*(x0-x3)*(x0-x4)*(x0-x5)*(x0-x6)*(x0-x8)/((x7-x1).*(x7-x2).*(x7-x3).*(x7-x4).*(x7-x5).*(x7-x6).*(x7-x8))*y7+ ...
(x0-x1)*(x0-x2)*(x0-x3)*(x0-x4)*(x0-x5)*(x0-x6)*(x0-x7)/((x8-x1).*(x8-x2).*(x8-x3).*(x8-x4).*(x8-x5).*(x8-x6).*(x8-x7))*y8;

end

if(N==7)
y0=(x0-x2)*(x0-x3)*(x0-x4)*(x0-x5)*(x0-x6)*(x0-x7)/((x1-x2).*(x1-x3).*(x1-x4).*(x1-x5).*(x1-x6).*(x1-x7))*y1+ ...
(x0-x1)*(x0-x3)*(x0-x4)*(x0-x5)*(x0-x6)*(x0-x7)/((x2-x1).*(x2-x3).*(x2-x4).*(x2-x5).*(x2-x6).*(x2-x7))*y2+ ...
(x0-x1)*(x0-x2)*(x0-x4)*(x0-x5)*(x0-x6)*(x0-x7)/((x3-x1).*(x3-x2).*(x3-x4).*(x3-x5).*(x3-x6).*(x3-x7))*y3+ ...
(x0-x1)*(x0-x2)*(x0-x3)*(x0-x5)*(x0-x6)*(x0-x7)/((x4-x1).*(x4-x2).*(x4-x3).*(x4-x5).*(x4-x6).*(x4-x7))*y4+ ...
(x0-x1)*(x0-x2)*(x0-x3)*(x0-x4)*(x0-x6)*(x0-x7)/((x5-x1).*(x5-x2).*(x5-x3).*(x5-x4).*(x5-x6).*(x5-x7))*y5+ ...
(x0-x1)*(x0-x2)*(x0-x3)*(x0-x4)*(x0-x5)*(x0-x7)/((x6-x1).*(x6-x2).*(x6-x3).*(x6-x4).*(x6-x5).*(x6-x7))*y6+ ...
(x0-x1)*(x0-x2)*(x0-x3)*(x0-x4)*(x0-x5)*(x0-x6)/((x7-x1).*(x7-x2).*(x7-x3).*(x7-x4).*(x7-x5).*(x7-x6))*y7;
end
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
.
4.
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 program when you will run it with hardcoded parameters.
MCMean =
0.100203838223922
Original process average from our simulation
ItoHermiteMean =
0.099629749702574
true Mean only applicble to standard SV mean reverting type models otherwise disregard
TrueMean =
0.100000000000000
IndexMax =
1302

Here is the output graph rescaled.

Again sometimes for difficult parameters SDEs, you may have to decrease step size a bit. This new program is a lot faster than previous programs.

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

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

Please protest against use of gases by mind control agencies on innocent civilians who only want to do good research. Read here:
https://forum.wilmott.com/viewtopic.php?f=15&t=94796&p=859944#p859944

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

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

Though my work on Fokker-Planck equation gives a reasonably nice working algorithm, it is not a completely conclusive solution but I will still take a break from the project to move on to two other interesting projects.
First project is to complete my work on path integrals on transition probabilities grid. I have some very good ideas and I hope that I would be able to calculate arithmetic path integrals, dz-integrals and dt-integrals for general SDEs on transition probabilities grid. I plan to calculate transition probabilities on the same expanding/contracting grid that I have used in the solution to Fokker-planck equations of SDEs using Ito-hermite method. Once we have solved FPE of SV type and other SDEs, we will be able to use this project to caclulate the integrals of variance in order to plug the variance integral values in lognormal/CEV/displaced diffusion/normal formulas to calculate the distribution of asset price and also do more interesting analytics.

The second project is to write a general program for monte carlo simulation of SDEs with explicit time dependence. In an earlier post I mentioned how Ito-Taylor method is equally applicable to SDEs that have an explicit time dependence in the drift or volatility. Read here for details:
https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=885#p855709
Since I think it would be nice and valuable, I have decided to write a very general program with explicitly time dependent coefficients. Here I have given the proposed general SDE with explicit time dependence
https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=885#p855857

I had previously made the following progress on the explicit time dependent SDEs project. Here I write the progress and my proposed strategy to solve iterated integrals when time explicit time dependent functions are present in the SDE. As I had mentioned, I had chosen the path of representation of time functions in the form of their Taylor series. Until now I have been dealing with taylor series of exponential function. I have written programs for calculation of variance and other stocahstic integrals. After calculation of variance, I re-convert the variance series back into "volatility series" by taking the square-root of the power series. Square root of a power series usually does not exist if the first few terms of the power series are zero which is mostly the case for the variance series but I just take some required powers of t common out of the series so that leading term is non-zero positive and then we can take the square-root of the variance series to convert it into a volatility series and use it again for higher order calculations of stochastic integrals and this works very well for the exponential functions in the drift and volatility.
Please continue to follow the discussion on the thread, hopefully there is going to be more interesting stuff in coming days.

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

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

Friends, I have had good success in calculating the arithmetic path integrals of the stochastic differential equations on the transition probabilities grid. I have calculated transition probabilities between any two points on the expanding/contracting/drifting grid and initial results to calculate the integrals are very good. It would not have been easily possible if we did not have a grid that keeps probability mass constant within different subdivisions. I hope to be posting a program tomorrow. It is a very fast program and for transition probabilities, there is no need for a miniscule time step as we can calculate transition probabilities between grids that are eight time steps apart or sixteen time steps apart and step size for transition probabilities grid is a multiple of time step we used to calculate the density of SDE using fokker-planck equation solution. I hope to post a working program in another day.

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

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

I have a working program that works very well for easy densities but works not so excellent for densities where mean changes very fast and where kappa is very high. I also still have to work out some integrals as some integrals I used are only approximate. But I will still post the program with explanatory comments tomorrow so friends can play with it on their own. It is a quite reasonable program but still has its shortcomings that I hope to overcome in next few days but I will be posting the preliminary version tomorrow,

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

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

I have posted the following graphs for path integrals of the mean-reverting SDEs. The governing SDE is

$dX(t)=\kappa (\theta- X(t)) dt+ \sigma X(t)^{\gamma} dZ(t)$

And the definition of path integral is   $PI(t+1)=PI(t)+X(t)^{\omega_1} dt$
This is different from dt-integrals as
For monte carlo simulation, I used the simple algorithm

for  t = 1:T

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

YY(1:paths)=YY(1:paths)+ .....Simulation of original SDE.

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

end

Though my algorithm is still not final, major error is coming off a stochastic integral which I could not calculate. I replaced it with another substitute solution which was only reasonably bad but becomes very bad at times. I am very sure that I would solve this integral in another day or two and should not be very difficult. I had promised to post the program so I decided to go ahead otherwise I would have liked to wait for two or three more days. I have mentioned and pointed out this bad integral in comments in sub-function named CalculateFdtDriftAndVol04. This bad integral creates problems with large kappa and very fast moving densities so please be patient. I am perfectly sure and promise friends that we would be solving all different path integrals, dt-integrals and dz-integrals perfectly well in a few weeks.

I will post the matlab program with comments in next post in a few minutes.

Here are some of the graphs of path integral densities comparison of analytic path integral density and monte carlo path integral density as we defined at start. Parameters are at title of the graph. Please notice omega1 in parameters which is the power which defines the function used for path integrals.

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

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

I am posting the program used to make graphs in the above post. Please note again as I mentioned in previous post that there is a bad stochastic integral that affects performance. There may be one or two other minor issues but I will fix everything in next two or three days and come with a new version. Here are the programs.
function [] = FPERevisitedTransProb04()

%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.
%In this version, the mean correction has been disabled. If you are
%simulating standard mean-reverting SV type SDEs, please enable the mean
%correction by uncommenting the appropriate line in the body of code below.

dt=.125/16/2/2;   % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=128*8/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=1.0;   % 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=1.0;%mean reversion target
sigma0=.750;%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]=CalculateFdtDriftAndVol04(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.

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 same
%grid point. It will continue to evolve in the same grid
%point till end.
Efdtnew(mm)=Efdtnew(mm)+Efdtprev(mm)+ds*.5*yy1(mm).^omega1.*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))
%Below calculate the function value and save it in the same
%grid point. It will continue to evolve in the same grid
%point till end.
Efdtnew(mm)=Efdtnew(mm)+Efdtprev(mm)+ds*yy1(mm).^omega1.*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)+ds*(.5*yy2(wnStart:Nn).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm)+ ...
.50/pi^2*dfdtGenerator(mm,wnStart:Nn).*pt(mm,wnStart:Nn).*ZProb(mm);

else
%Below calculate the function value and save it in the same
%grid point. It will continue to evolve in the same grid
%point till end.
Efdtnew(mm)=Efdtnew(mm)+Efdtprev(mm)+ds*yy1(mm).^omega1.*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);
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]=CalculateFdtDriftAndVol04(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);
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=100000;
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*8;%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 first sub-function used. This should ideally be the same as the one originally used in the solution of Fokker-planck and I will replace it with that in next versions.
.
function [fMu0dt,dfz01,df2z02]=CalculateFdtDriftAndVol04(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,w,wnStart,Nn,dt,t1,t2)
%dtsqrt=sqrt(dt);
%dtonehalf=dt^1.5;
%dt2=dt^2/2;
%dt3=dt^3/6.0;
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^2*t2^3/3+t1^5/3-t1^5/5)) + ...
1* (d_dfDD(wnStart:Nn).*VV(wnStart:Nn) + ...
.5*d_d2fQQ(wnStart:Nn).*VV(wnStart:Nn)).*((t2^2.5-t1^2.5)).*((1-1/sqrt(3))*1/sqrt(2)*1/2)/2;
%The very last integral in above expression is totally bad and badly
%affects all calculations especially when kappa is high or when density
%moves around. I will fix it in next version.

%.5*d_d2fQQ(wnStart:Nn).*VV(wnStart:Nn)).*sqrt((1*t2^5/20-1*t1^5/20+t1*t2^4/4-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));

end
.
Here is the second sub-function used.
.
function [wMu0dt,dwMu0dt,dw0]=CalculateDriftAndVol(mu1,mu2,beta1,beta2,gamma,sigma0,w,wnStart,Nn,dt)
dtsqrt=sqrt(dt);
dtonehalf=dt^1.5;
dt2=dt^2/2;
dt3=dt^3/6.0;
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
theta1=mu1;
theta2=mu2;
theta3=.5*sigma0^2*(-gamma);
omega1=beta1-gamma;
omega2=beta2-gamma;
omega3=gamma-1;
GG(wnStart:Nn)=theta1*yy(wnStart:Nn).^omega1+theta2*yy(wnStart:Nn).^omega2+theta3*yy(wnStart:Nn).^omega3;
dGG_dw(wnStart:Nn)=theta1*(1-gamma).^(omega1/(1-gamma))*(omega1/(1-gamma))*w(wnStart:Nn).^(omega1/(1-gamma)-1)+ ...
theta2*(1-gamma).^(omega2/(1-gamma))*(omega2/(1-gamma))*w(wnStart:Nn).^(omega2/(1-gamma)-1)+ ...
theta3*(1-gamma).^(omega3/(1-gamma))*(omega3/(1-gamma))*w(wnStart:Nn).^(omega3/(1-gamma)-1);

dGGDD_dw(wnStart:Nn)= ...
mu1*omega1*theta1.*(omega1+beta1-1).*((1-gamma).*w(wnStart:Nn)).^((omega1+beta1-1)/(1-gamma)-1)+ ...
mu1*omega2*theta2.*(omega2+beta1-1).*((1-gamma).*w(wnStart:Nn)).^((omega2+beta1-1)/(1-gamma)-1)+ ...
mu1*omega3*theta2.*(omega3+beta1-1).*((1-gamma).*w(wnStart:Nn)).^((omega3+beta1-1)/(1-gamma)-1)+ ...
mu2*omega1*theta1.*(omega1+beta2-1).*((1-gamma).*w(wnStart:Nn)).^((omega1+beta2-1)/(1-gamma)-1)+ ...
mu2*omega2*theta2.*(omega2+beta2-1).*((1-gamma).*w(wnStart:Nn)).^((omega2+beta2-1)/(1-gamma)-1)+ ...
mu2*omega3*theta2.*(omega3+beta2-1).*((1-gamma).*w(wnStart:Nn)).^((omega3+beta2-1)/(1-gamma)-1);

dGG2QQ_dw(wnStart:Nn)= ...
theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*((1-gamma).*w(wnStart:Nn)).^((omega1+2*gamma-2)/(1-gamma)-1)+ ...
theta2*sigma0^2.*omega2*(omega2-1).*(omega2+2*gamma-2).*((1-gamma).*w(wnStart:Nn)).^((omega2+2*gamma-2)/(1-gamma)-1)+ ...
theta3*sigma0^2.*omega3*(omega3-1).*(omega3+2*gamma-2).*((1-gamma).*w(wnStart:Nn)).^((omega3+2*gamma-2)/(1-gamma)-1);

dGG(wnStart:Nn)=omega1*theta1*yy(wnStart:Nn).^(omega1-1)+omega2*theta2*yy(wnStart:Nn).^(omega2-1)+ ...
omega3*theta3*yy(wnStart:Nn).^(omega3-1);
DD(wnStart:Nn)=mu1*yy(wnStart:Nn).^beta1+mu2*yy(wnStart:Nn).^beta2;
dGG2(wnStart:Nn)=theta1*omega1*(omega1-1).*yy(wnStart:Nn).^(omega1-2)+ ...
theta2*omega2*(omega2-1).*yy(wnStart:Nn).^(omega2-2)+ ...
theta3*omega3*(omega3-1).*yy(wnStart:Nn).^(omega3-2);
QQ(wnStart:Nn)=sigma0^2*yy(wnStart:Nn).^(2*gamma);
ddGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*yy(wnStart:Nn).^(omega1+beta1-2)+ ...
mu1*omega2*theta2*(omega2+beta1-1).*yy(wnStart:Nn).^(omega2+beta1-2)+ ...
mu1*omega3*theta3*(omega3+beta1-1).*yy(wnStart:Nn).^(omega3+beta1-2)+ ...
mu2*omega1*theta1*(omega1+beta2-1).*yy(wnStart:Nn).^(omega1+beta2-2)+ ...
mu2*omega2*theta2*(omega2+beta2-1).*yy(wnStart:Nn).^(omega2+beta2-2)+ ...
mu2*omega3*theta3*(omega3+beta2-1).*yy(wnStart:Nn).^(omega3+beta2-2);
d2dGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*(omega1+beta1-2).*yy(wnStart:Nn).^(omega1+beta1-3)+ ...
mu1*omega2*theta2*(omega2+beta1-1).*(omega2+beta1-2).*yy(wnStart:Nn).^(omega2+beta1-3)+ ...
mu1*omega3*theta3*(omega3+beta1-1).*(omega3+beta1-2).*yy(wnStart:Nn).^(omega3+beta1-3)+ ...
mu2*omega1*theta1*(omega1+beta2-1).*(omega1+beta2-2).*yy(wnStart:Nn).^(omega1+beta2-3)+ ...
mu2*omega2*theta2*(omega2+beta2-1).*(omega2+beta2-2).*yy(wnStart:Nn).^(omega2+beta2-3)+ ...
mu2*omega3*theta3*(omega3+beta2-1).*(omega3+beta2-2).*yy(wnStart:Nn).^(omega3+beta2-3);
ddGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*yy(wnStart:Nn).^(omega1+2*gamma-3)+ ...
sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*yy(wnStart:Nn).^(omega2+2*gamma-3)+ ...
sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*yy(wnStart:Nn).^(omega3+2*gamma-3);
d2dGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*(omega1+2*gamma-3).*yy(wnStart:Nn).^(omega1+2*gamma-4)+ ...
sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*(omega2+2*gamma-3).*yy(wnStart:Nn).^(omega2+2*gamma-4)+ ...
sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*(omega3+2*gamma-3).*yy(wnStart:Nn).^(omega3+2*gamma-4);

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

wMu0dt01(wnStart:Nn)=GG(wnStart:Nn)*dt+ dGG(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...
.5*dGG2(wnStart:Nn).*QQ(wnStart:Nn).*dt2+ ...
ddGGDD(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
.5*ddGG2QQ(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
.5*d2dGGDD(wnStart:Nn).*QQ(wnStart:Nn).*dt3 + ...
.25*d2dGG2QQ(wnStart:Nn).*QQ(wnStart:Nn).*dt3;
wMu0dt(wnStart:Nn)=GG(wnStart:Nn)*dt;

dwMu0dt(wnStart:Nn)=dGG_dw(wnStart:Nn)*dt+dGGDD_dw(wnStart:Nn)*dt2 +...
.5*dGG2QQ_dw(wnStart:Nn).*dt2;

dw0(1:Nn)=sigma0.*dtsqrt+1*dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*dtonehalf;

end


.
Here is the third sub-function used.
.
function [fMu0dt,dfz0,df2z0]=CalculateFunctionDriftAndVol(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,w,wnStart,Nn,dt)
dtsqrt=sqrt(dt);
dtonehalf=dt^1.5;
dt2=dt^2/2;
dt3=dt^3/6.0;
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;
dDD(wnStart:Nn)=mu1*beta1.*yy(wnStart:Nn).^(beta1-1)+mu2*beta2*yy(wnStart:Nn).^(beta2-1);

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;
dVV(wnStart:Nn)=sigma0.*gamma*yy(wnStart:Nn).^(gamma-1);
d2VV(wnStart:Nn)=sigma0.*gamma*(gamma-1).*yy(wnStart:Nn).^(gamma-2);

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);
d3_dfVV(wnStart:Nn)= ...
sigma0.*omega1.*theta1.*(omega1+gamma-1).*(omega1+gamma-2).*(omega1+gamma-3).*yy(wnStart:Nn).^(omega1+gamma-4);

fMu0dt(wnStart:Nn)=df(wnStart:Nn).*DD(wnStart:Nn).*dt + ...
.5*d2f(wnStart:Nn).*QQ(wnStart:Nn).*dt+ ...
d_dfDD(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...
.5*d_d2fQQ(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...
.5*d2_dfDD(wnStart:Nn).*QQ(wnStart:Nn).*dt2 + ...
.25*d2_d2fQQ(wnStart:Nn).*QQ(wnStart:Nn).*dt2;

dfz0(1:Nn)=df(wnStart:Nn).*VV(wnStart:Nn).*dtsqrt+ ...
d_dfVV(wnStart:Nn).*DD(wnStart:Nn)*1/sqrt(3).*dt^1.5+ ...
.5*d2_dfVV(wnStart:Nn).*QQ(wnStart:Nn)*1/sqrt(3).*dt^1.5+ ...
d_dfDD(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3)).*dt^1.5 + ...
.5*d_d2fQQ(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3)).*dt^1.5;

df2z0(1:Nn)=d_dfVV(wnStart:Nn).*VV(wnStart:Nn).*dt/2+ ...
(d2_dfVV(wnStart:Nn).*VV(wnStart:Nn)+d_dfVV(wnStart:Nn).*dVV(wnStart:Nn)).* ...
DD(wnStart:Nn)./sqrt(3)/2/sqrt(2).*dt^2.0+ ...
.5*(d3_dfVV(wnStart:Nn).*VV(wnStart:Nn)+d2_dfVV(wnStart:Nn).*dVV(wnStart:Nn)+ ...
d2_dfVV(wnStart:Nn).*dVV(wnStart:Nn)+d_dfVV(wnStart:Nn).*d2VV(wnStart:Nn)).* ...
QQ(wnStart:Nn)./sqrt(3)/2/sqrt(2).*dt^2.0+ ...
(d2_dfVV(wnStart:Nn).*DD(wnStart:Nn)+d_dfVV(wnStart:Nn).*dDD(wnStart:Nn)).* ...
VV(wnStart:Nn)*(1-1/sqrt(3))*1/2/sqrt(2)*dt^2.0+ ...
(d2_dfDD(wnStart:Nn).*VV(wnStart:Nn)+d_dfDD(wnStart:Nn).*dVV(wnStart:Nn)).* ...
VV(wnStart:Nn)*1/2*(1-sqrt(2)/2)*dt^2.0;

[size=85][font=Helvetica Neue, Helvetica, Arial, sans-serif]    [/font][/size]
end
Other sub-functions are the same as used for earlier solution of fokker-planck program and I have posted them in previous posts.
Here is the output of the posted program.
MCMean =
0.998162078044431
Original process average from our simulation
ItoHermiteMean =
0.999989249549945
true Mean only applicble to standard SV mean reverting type models otherwise disregard
TrueMean =
1
fYYdtm =
2.000085186273280
fydtm =
2.000010348793002
fYYdtvar =
0.512380044123777
fydtvar =
0.515820772698450
IndexMax =
483
red line is arithmetic dt integral density from transition probability framework , green is monte carlo.
IndexMax =
407
red line is density of SDE from Ito-Hermite method, green is monte carlo.

Amin
Topic Author
Posts: 2721
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.