SERVING THE QUANTITATIVE FINANCE COMMUNITY

 
User avatar
Amin
Topic Author
Posts: 2538
Joined: July 14th, 2002, 3:00 am

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

October 31st, 2020, 3:09 pm

In order to better understand the grid and evolve the functionals of SDEs on it, I have decided to first be able to replicate the distributionof  all stochastic integrals of the form
[$]\int_{t_1}^{t_2} dz(t)[$], [$]\int_{t_1}^{t_2} \int_{t_1}^{t} ds dz(t)[$], [$]\int_{t_1}^{t_2} \int_{t_1}^{t} dz(s) dt[$], [$]\int_{t_1}^{t_2} \int_{t_1}^{t} dz(s) dz(t)[$], and similar other more complicated integrals that are required for evolution of dt-integrals and dz-integrals and match all of these stochastic integrals first with their monte carlo counterparts and try to perfectly reproduce them in our transition probabilities framework. I will have to see if any adjustment is required to properly calculate these integrals. And then we could calculate evolution of these integrals as a function of our Zt(mm,nn) and use the actual value of these integrals as Zt1(t1.mm) advances to Zt2(t2,nn) meaning we will use the precise values (change in particular values of stochastic integrals on point mm on first time grid and point nn on second time grid) of these integrals that corresponds to change in cumulative Z from one point on first time grid to arbitrary second point on the next time grid. I really hope that this exercise will be helpful in understanding the grid better.
In a few days, I will come back with a matlab program comparing these stochastic integrals on transition probabilities grid and those from monte carlo simulation.
 
User avatar
Amin
Topic Author
Posts: 2538
Joined: July 14th, 2002, 3:00 am

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

November 3rd, 2020, 1:28 pm

I was not able to work day before yesterday as I was hit by drugged food. And was able to do some productive work yesterday and today. I went over the program and calculated the integral of generating unit normals driving the mean reverting SDEs both in monte carlo and in transition probabilities framework. First I calculated the conditional expected value of the integral of driving unit normal noise as a function of SDE variable y(t) in original coordinates and then I plotted the integrated normal noise on a linear scale. I also plotted an analytically calculated normal that had variance equal to number of time intervals. I noticed that conditional expected distribution values  of analytically calculated normal and monte carlo simulation integrated normal were almost totally identical. However transition probabilities calculated integrals of driving normal noise lost variance with time which became quite severe as time increased.
Since in plotting conditional expected value of analytically calculated integral of driving noise, I used fokker-planck solution of SDE as x-axis base and it was almost totally identical with monte carlo calculated conditional expected value of integral of driving noise, I think that if the original solution of the SDE by fokker-planck is exact enough, the better way to generate transition probabilities and driving noise would be to use analytic values of integrated driving noise on the underlying normal grid(variance equal to number of time steps) and use that to generate transition probabilities and transition normals between any two points on adjacent grids. We are using a unit scaling for this but time scaling can also be used since transition probabilities between any two times would then have a dt scaling instead of a unit scaling we have now. This can be a very fast and exact method of producing transition probabilities and driving noise for the transition probabilities grid.
In the program, I am posting there is a new function that calculates graph of conditional expected distributed value of one correlated variable as a function of other correlated variable. I wrote this function more than five years ago and it was a sister function to density generating function but I never got to use this function very much. I am sure friends would find it quite useful and handy in their research and I regret not posting it earlier since I think it can be a very helpful function in many ways. It takes one variable on x-axis and returns conditional expected value of other correlated variable on y-axis when both of the variables have been calculated as a pair in each of the path of large number of monte carlo simulations.
Now the graphs
 
There are two sets of graphs, graph-1 and graph-2 on various time levels.

Graph-1 t=8/512
The  graph below plots conditional expected value of integral of unit normals as a function of y(t) which is the value of SDE variable in original coordinates. Blue graph is true analytic normal as calculated from variance, red graph is conditional value of integral calculated from the transition probabilities grid and green value is conditional value of integral of normals as a function of y(t) calculated from monte carlo simulations. Please note that all three graphs are initially very close but red graph slowly becomes different and has lesser variance as time increases. While blue graph of true analytic normal calculated from variance and green graph of conditional integral calculated from monte carlo remain almost exactly the same. Please note that graph of integrated normal will not be a straight line as it is a non-linear function of y(t). Please also note that blue graph of true analytic normal that coincides with monte carlo graph is plotted conditional of values of y(t) calculated in the analytical fokker-planck solution of the SDE.

Image


Once we have established from previous conditional values(as you will see in sets of graph-1 below for larger time levels) graphs-1 that perfect normal blue graph calculated from variance is almost identical with monte carlo calculated integral of normal, we now plot it (graphs-2) on an x-axis value for which normal has to be a linear line so we can compare the true value of blue graph with red graph of our defective transition probabilities calculated integral of normal in a more familiar way. These are second sets of graphs as graph-2.

Graph-2 t=8/512

Image
In the first set of graphs at very initial date all lines are identical since everything is accurate. 
Now we follow with graphs at date  t=128/512.
Graph-1 date t=128/512.
You can see in this graph that analytically calculated normal(blue line) at y-axis with fokker-planck calculated SDE on X-axis is almost identical with conditional value of monte carlo calculated integrated normal(green graph) while variance of transition probabilities calculated normal(red line) has started to lag very slightly.
Image
Linear Graph2 t=128/512
Image
 Graph-1 t=256/512

Here we see that analytically calculated normal conditional values(blue) match with monte carlo calculated conditional values(green) but transition probabilities grid calculated values have started to lag more behind.

Image
Linear version Graph-2 t=256/512

Image
Graph-1 t=384/512
Image
Linear Graph-2 t=384/512

Image

Now I give the output of the program at the end of last time interval t=512/512

ss =
   512
IndexMax =
    51
MCMean =
   1.249844261079288                %Monte carlo mean of the SDE
MCVar =
   0.593135335310294                %Monte carlo variance of the SDE
Original process average from our simulation
ItoHermiteMean =
   1.251551446212164                %Mean of Fokker-planck solution of SDE
ItoHermiteVar =
   0.597372214584988               %Variance of fokker-planck solution of SDE
MuZw =
   0.037044454846106               %Mean of transition probabilities sum of SDE driving unit normals
VZw =
     3.559456867117647e+02      %Variance of transition probabilities sum of SDE driving unit normals
MuIZM =
  -0.020600972280994              %Mean of monte carlo sum of  sum of SDE driving unit normals
VIZM =
     5.127728403892854e+02   %Variance of monte carlo sum of SDE driving unit normals

We see that monte carlo sum of SDE driving unit normals(green line) after 512 intervals is 512.77 while transition probabilities calculated sum(red line) is 355.94. Analytically calculated sum(blue line) is exactly 512. We show conditional values of all three sums in the graph-1 below
Please note that blue line and green line are very close in main body of the graph where monte carlo variance is little. Please also note that blue graph which shows analytic sum on y-axis uses fokker-planck calculated values of y(t)(SDE variable in original coordinates) at x-axis which means that we can use the analytically calculated normal for both transition probabilities and unit variance increments between two adjacent time intervals since blue line when based on fokker-planck SDE variable y(t) at x-axis makes it almost identical with conditional value of monte carlo integrated noise as shown by green line. 
graph-1 t=512/512 


Image
Graph2- t=512/512 linear comparisonof normals
Image
 
User avatar
Amin
Topic Author
Posts: 2538
Joined: July 14th, 2002, 3:00 am

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

November 3rd, 2020, 1:30 pm

Here is the matlab program used for the previous post. It is very slow and is for diagnostic purposes only. It has monte carlo simulations and analytical part embedded in the same time loop so we can compare values on each time step.
function [] = FPERevisitedTransProb08DNew08Wmt02()

%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334 
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=mu1 x(t)^beta1 dt + mu2 x(t)^beta2 dt +sigma x(t)^gamma dz(t)

%I have not directly simulated the SDE but simulated the transformed 
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.

dt=.125/16/2/2;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.
Tt=128*4/4*2*2;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year; 
%Tt=8*60;
T=Tt*dt;
OrderA=4;  %
OrderM=4;  %
%dtM=.125/8;%Monte carlo time interval size dtM.
%TtM=8*8;%Monte carlo number of simulation intervals.
dtM=dt;%*4*2*2;
TtM=Tt;%/4/2/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.2500;   % 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.250;%mean reversion target
sigma0=.850;%Volatility value
omega1=1.50; %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.

%if you alter above parameters and the program blows, disable
%mean-correction since that is only valid for classical mean-reverting
%SDEs.

%you can specify any general mu1 and mu2 and beta1 and beta2.
mu1=+1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
%mu1=0;
%mu2=0;

alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;%This is for expansion of integrals for calculation of drift 
%and volatility coefficients
                
w(1:Nn)=x0^(1-gamma)/(1-gamma);

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

Z
%Znew=sqrt(2)*Z
%sqrt(Znew.^2-Z.^2)
%Znew-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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Current program has 50 grid points. I have divided every grid subdivision
%into three further subdivisions with probabilities P0,Pb,Pa.
 P0(1:Nn)=normcdf(Z(1:Nn)+1/6*dNn)-normcdf(Z(1:Nn)-1/6*dNn);
        Pa(1:Nn-1)=normcdf(Z(1:Nn-1)+1/2*dNn)-normcdf(Z(1:Nn-1)+1/6*dNn);
        Pb(1+1:Nn)=normcdf(Z(1+1:Nn)-1/6*dNn)-normcdf(Z(1+1:Nn)-1/2*dNn);
        Pb(1)=normcdf(Z(1)-1/6*dNn);
        Pa(Nn)=1-normcdf(Z(Nn)+1/6*dNn);
        PSum(1:Nn)=P0(1:Nn)+Pa(1:Nn)+Pb(1:Nn);

ss=0;  %ss is index for transition probability calculation time steps.
Freq=1.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;
   PZprev(1:Nn)=ZProb(1:Nn);%probability mass at starting node.
   PZnew(1:Nn)=0.0;
   Eyyprev(1:Nn)=yy(1:Nn).*ZProb(1:Nn);%Value of SDE variable in original 
   %coordinates at starting nodes.
   Eyynew(1:Nn)=0;
   EZwnew(1:Nn)=0.0;
   EZwprev(1:Nn)=0.0;
   
   
   Efdtprev(1:Nn)=0;%Value of path-integral at starting nodes is zero.
   Efdtnew(1:Nn)=0;
   Egdtprev(1:Nn)=0;%Value of path-integral at starting nodes is zero.
   Egdtnew(1:Nn)=0;
   
   Ef0prev(1:Nn)=0;
   %Ef0prev(Nn0)=0;%Value of noise at initial node is zero.
   Ef0new(1:Nn)=0;
   EIZAprev(1:Nn)=0.0;
   EIZAnew(1:Nn)=0.0;
   
   
    %Monte carl ovariables below
    rng(29079137, 'twister')
    paths=1000000;
    YY(1:paths)=x0;  %Original process monte carlo.
    fYYdt(1:paths)=0.0;
    Random1(1:paths)=0;
    IZM(1:paths)=0.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);
    [wMu0dt,c1,c22] = CalculateDriftAndVolA8Trans(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;%+c22(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
    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*.025;
        C33(wnStart:Nn)=.0125/pi*2*.35; 
        C44(wnStart:Nn)=.0125/pi*25.0*.850;
        
        
        
        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)&&(tt<=128))
         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
     
     if (tt>128)
         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


  %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);
        [wMu0dt2,dw02,c22] = CalculateDriftAndVolA8Trans(W1,wnStart,Nn,YqCoeff0,Fp1,gamma,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,w,wnStart,Nn,dt);
        [fMu0dt,dfz0,df2z0]=CalculateFunctionDriftAndVol(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds);
        
        %below calculates integrals for drift and volatility for
        %yy(s)^omega1 ds. which is used in path integrals calculations
        %on the lattice.
        %[fdtMu0dt,dfdtz01,dfdt2z02]=CalculateFdtDriftAndVol08(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds,t1A,t2A);

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

     %   [fMu0dtdt,dfz0dt,df2z0dt]=CalculateFunctionDriftAndVoldt(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds);

    end

    if(rem(tt,Freq)==0)
    
        %When this above condition is true it means that we are at end of
        %particular transition probability interval and we have end of 
        %interval W2(s)=w(t) while we will use the value W1 that was 
        %previously saved Freq dt time intervals earlier or ds=Freq*dt 
        %time ago.
        
        ss=ss+1; %transition loop variable and is one at the end of first 
        %transition probability interval.
        %Div0=Div0+.5*ds;
        
        
        W2(1:Nn)=w(1:Nn);
        
        ds=dt*Freq;
         
        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.
        
        Yy(ss,1:Nn)=yy2(1:Nn);

        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 

        
        [W2GridStarts,W2GridEnds] = CalculateGridStartsAndEnds(W2,Z,wnStart,Nn,dNn);
        [W1GridStarts,W1GridEnds] = CalculateGridStartsAndEnds(W1,Z,wnStart,Nn,dNn);
         
        %dYy are grid widths in original coordinates. will be handy later.
        dYy(ss,2:Nn-1)=((1-gamma).*W2GridEnds(2:Nn-1)).^(1/(1-gamma))-((1-gamma).*W2GridStarts(2:Nn-1)).^(1/(1-gamma));
        dYy(ss,1)=dYy(ss,2);
        dYy(ss,Nn)=dYy(ss,Nn-1);
        W2GridStarts0(1:Nn)=W2GridStarts(1:Nn);
        W2GridEnds0(1:Nn)=W2GridEnds(1:Nn);
        W2GridStarts0(1)=W2(1)-(W2(2)-W2(1));
        W2GridEnds0(Nn)=W2(Nn)+(W2(Nn)-W2(Nn-1));
        
        dW2(2:Nn-1)=W2GridEnds(2:Nn-1)-W2GridStarts(2:Nn-1);
        dW2(1)=dW2(2);
        dW2(Nn)=dW2(Nn-1);
        
        dW1(2:Nn-1)=W1GridEnds(2:Nn-1)-W1GridStarts(2:Nn-1);
        dW1(1)=dW1(2);
        dW1(Nn)=dW1(Nn-1);
        
        
        
        
        
        %initial grid from where transition probabilities have to be
        %calculated is divided into three further subdivisions. Initial
        %grid subdivision ranges from w(Z_n-.5*dNn) to w(Z_n+.5*dNn). Now
        %it is divided into three further subdivisions with left subdivision 
        %W1GridMidsb centred at w(Z_n-1/3*dNn) spaced between w(Z_n-.5*dNn) to 
        %w(Z_n-1/6*dNn). Second middle subdivision is centred at  w(Z_n) and is
        %spaced between w(Z_n-1/6*dNn) to w(Z_n+1/6*dNn). Third further 
        %subdivision called right subdivision W1GridMidsa is centred at 
        %w(Z_n+1/3*dNn) and ranges from w(Z_n+1/6*dNn) to w(Z_n+.5*dNn). Since our initial full 
        %subdivision can be quite large, We calculate transition
        %probabilities independently from three further subdivision
        %midpoints and then add them according to weight of probability
        %mass of each of the smaler three subdivisions to get one estimate of transition probability. This is why we
        %have to calculated W1GridMidsa(right subdivision)  and 
        %W1GridMidsb(left subdivision) located at w(Z_n+1/3*dNn)
        % and w(Z_n-1/3*dNn) respectively where W1 itself(middle subdivision)
        %is located at w(Z_n).
        
        [W1GridMidsa,W1GridMidsb] = CalculateGridSubdivisions(W1,Z,wnStart,Nn,dNn);
        
        %below we calculate drift and volatility for W1GridMidsa and
        %W1GridMidsb since these would be used in independent calculations
        %of transition probabilities from these points W1GridMidsa and
        %W1GridMidsb
        [wMu0dt2,dw02,c22] = CalculateDriftAndVolA8Trans(W1,wnStart,Nn,YqCoeff0,Fp1,gamma,ds);
        [wMu0dt2a,dw02a,c2a] = CalculateDriftAndVolA8Trans(W1GridMidsa,wnStart,Nn,YqCoeff0,Fp1,gamma,ds);
        [wMu0dt2b,dw02b,c2b] = CalculateDriftAndVolA8Trans(W1GridMidsb,wnStart,Nn,YqCoeff0,Fp1,gamma,ds);

        
        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);
        %[wMu0dt2,c12,c22] = CalculateDriftAndVolA8Trans(W1,wnStart,Nn,YqCoeff0,Fp1,gamma,ds);
    
           %I have used the formula
           %W2=W1+ mu* dt+ sigma *Z + sigma2 *(Z^2-1) and back out Z by use
           %of quadratic formula.
           %I use A lesser known quadratic formula, as used in Muller's method (from wikipedia)
           %in order to obviate round-off errors
           %Here mu*dt=wMu0dt; sigma=dw02; sigma2 =c22; 
            if((W2(wnStart:Nn)-W1(mm)-1*wMu0dt2(mm)+c22(mm))>0)
                %Zt(mm,wnStart:Nn)=(-dw02(mm) + sqrt(dw02(mm).^2+4*c22(mm).*(-W2(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm))))./(2*c22(mm));
                Zt(mm,wnStart:Nn)=(2*(-W2(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm)))./(-dw02(mm) + sqrt(dw02(mm).^2+4*c22(mm).*(-W2(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm))));
            else
                Zt(mm,wnStart:Nn)=(2*(-W2(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm)))./(-dw02(mm) - sqrt(dw02(mm).^2-4*c22(mm).*(-W2(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm))));
            end
            
            
            if((W2GridStarts(wnStart:Nn)-W1(mm)-1*wMu0dt2(mm)+c22(mm))>0)
                %Zt(mm,wnStart:Nn)=(-dw02(mm) + sqrt(dw02(mm).^2+4*c22(mm).*(-W2(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm))))./(2*c22(mm));
                ZtS(mm,wnStart:Nn)=(2*(-W2GridStarts(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm)))./(-dw02(mm) + sqrt(dw02(mm).^2+4*c22(mm).*(-W2GridStarts(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm))));
            else
                ZtS(mm,wnStart:Nn)=(2*(-W2GridStarts(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm)))./(-dw02(mm) - sqrt(dw02(mm).^2-4*c22(mm).*(-W2GridStarts(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm))));
            end
            
            if((W2GridEnds(wnStart:Nn)-W1(mm)-1*wMu0dt2(mm)+c22(mm))>0)
                %Zt(mm,wnStart:Nn)=(-dw02(mm) + sqrt(dw02(mm).^2+4*c22(mm).*(-W2(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm))))./(2*c22(mm));
                ZtE(mm,wnStart:Nn)=(2*(-W2GridEnds(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm)))./(-dw02(mm) + sqrt(dw02(mm).^2+4*c22(mm).*(-W2GridEnds(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm))));
            else
                ZtE(mm,wnStart:Nn)=(2*(-W2GridEnds(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm)))./(-dw02(mm) - sqrt(dw02(mm).^2-4*c22(mm).*(-W2GridEnds(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm))));
            end
            
           
            %Below ZtS indicates the value of Zt at start of W2
            %subdivision and ZtE indicates the value of Zt at the end of W2
            %subdivision
            
            
            if((W2GridStarts(wnStart:Nn)-W1GridMidsa(mm)-1*wMu0dt2a(mm)+c2a(mm))>0)
                %Zt(mm,wnStart:Nn)=(-dw02(mm) + sqrt(dw02(mm).^2+4*c22(mm).*(-W2(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm))))./(2*c22(mm));
                ZtSa(mm,wnStart:Nn)=(2*(-W2GridStarts(wnStart:Nn)+W1GridMidsa(mm)+1*wMu0dt2a(mm)-c2a(mm)))./ ...
                    (-dw02a(mm) + sqrt(dw02a(mm).^2+4*c2a(mm).*(-W2GridStarts(wnStart:Nn)+W1GridMidsa(mm)+1*wMu0dt2a(mm)-c2a(mm))));
            else
                ZtSa(mm,wnStart:Nn)=(2*(-W2GridStarts(wnStart:Nn)+W1GridMidsa(mm)+1*wMu0dt2a(mm)-c2a(mm)))./ ...
                    (-dw02a(mm) - sqrt(dw02a(mm).^2-4*c2a(mm).*(-W2GridStarts(wnStart:Nn)+W1GridMidsa(mm)+1*wMu0dt2a(mm)-c2a(mm))));
            end
            
            if((W2GridEnds(wnStart:Nn)-W1GridMidsa(mm)-1*wMu0dt2a(mm)+c2a(mm))>0)
                %Zt(mm,wnStart:Nn)=(-dw02(mm) + sqrt(dw02(mm).^2+4*c22(mm).*(-W2(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm))))./(2*c22(mm));
                ZtEa(mm,wnStart:Nn)=(2*(-W2GridEnds(wnStart:Nn)+W1GridMidsa(mm)+1*wMu0dt2a(mm)-c2a(mm)))./ ...
                    (-dw02a(mm) + sqrt(dw02a(mm).^2+4*c2a(mm).*(-W2GridEnds(wnStart:Nn)+W1GridMidsa(mm)+1*wMu0dt2a(mm)-c2a(mm))));
            else
                ZtEa(mm,wnStart:Nn)=(2*(-W2GridEnds(wnStart:Nn)+W1GridMidsa(mm)+1*wMu0dt2a(mm)-c2a(mm)))./ ...
                    (-dw02a(mm) - sqrt(dw02a(mm).^2-4*c2a(mm).*(-W2GridEnds(wnStart:Nn)+W1GridMidsa(mm)+1*wMu0dt2a(mm)-c2a(mm))));
            end

            
            if((W2GridStarts(wnStart:Nn)-W1GridMidsb(mm)-1*wMu0dt2b(mm)+c2b(mm))>0)
                %Zt(mm,wnStart:Nn)=(-dw02(mm) + sqrt(dw02(mm).^2+4*c22(mm).*(-W2(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm))))./(2*c22(mm));
                ZtSb(mm,wnStart:Nn)=(2*(-W2GridStarts(wnStart:Nn)+W1GridMidsb(mm)+1*wMu0dt2b(mm)-c2b(mm)))./ ...
                    (-dw02b(mm) + sqrt(dw02b(mm).^2+4*c2b(mm).*(-W2GridStarts(wnStart:Nn)+W1GridMidsb(mm)+1*wMu0dt2b(mm)-c2b(mm))));
            else
                ZtSb(mm,wnStart:Nn)=(2*(-W2GridStarts(wnStart:Nn)+W1GridMidsb(mm)+1*wMu0dt2b(mm)-c2b(mm)))./ ...
                    (-dw02b(mm) - sqrt(dw02b(mm).^2-4*c2b(mm).*(-W2GridStarts(wnStart:Nn)+W1GridMidsb(mm)+1*wMu0dt2b(mm)-c2b(mm))));
            end
            
            if((W2GridEnds(wnStart:Nn)-W1GridMidsb(mm)-1*wMu0dt2b(mm)+c2b(mm))>0)
                %Zt(mm,wnStart:Nn)=(-dw02(mm) + sqrt(dw02(mm).^2+4*c22(mm).*(-W2(wnStart:Nn)+W1(mm)+1*wMu0dt2(mm)-c22(mm))))./(2*c22(mm));
                ZtEb(mm,wnStart:Nn)=(2*(-W2GridEnds(wnStart:Nn)+W1GridMidsb(mm)+1*wMu0dt2b(mm)-c2b(mm)))./ ...
                    (-dw02b(mm) + sqrt(dw02b(mm).^2+4*c2b(mm).*(-W2GridEnds(wnStart:Nn)+W1GridMidsb(mm)+1*wMu0dt2b(mm)-c2b(mm))));
            else
                ZtEb(mm,wnStart:Nn)=(2*(-W2GridEnds(wnStart:Nn)+W1GridMidsb(mm)+1*wMu0dt2b(mm)-c2b(mm)))./ ...
                    (-dw02b(mm) - sqrt(dw02b(mm).^2-4*c2b(mm).*(-W2GridEnds(wnStart:Nn)+W1GridMidsb(mm)+1*wMu0dt2b(mm)-c2b(mm))));
            end
            ZtS(mm,1)=-Inf;
            ZtSa(mm,1)=-Inf;
            ZtSb(mm,1)=-Inf;
            
            ZtS(mm,1)=ZtE(mm,1)-4;
            ZtSa(mm,1)=ZtEa(mm,1)-4;
            ZtSb(mm,1)=ZtEb(mm,1)-4;

            
            ZtE(mm,Nn)=ZtS(mm,Nn)+4;
            ZtEa(mm,Nn)=ZtSa(mm,Nn)+4;
            ZtEb(mm,Nn)=ZtSb(mm,Nn)+4;
            
            %Since ZtS indicates the value of Zt at start of W2
            %subdivision and ZtE indicates the value of Zt at the end of W2
            %subdivision, the probability mass inside can be found by
            %normcdf
            pt0(mm,wnStart:Nn)=normcdf(ZtE(mm,wnStart:Nn))-normcdf(ZtS(mm,wnStart:Nn));
            pta(mm,wnStart:Nn)=normcdf(ZtEa(mm,wnStart:Nn))-normcdf(ZtSa(mm,wnStart:Nn));
            ptb(mm,wnStart:Nn)=normcdf(ZtEb(mm,wnStart:Nn))-normcdf(ZtSb(mm,wnStart:Nn));
            %Zt(mm,wnStart:Nn)=(W2(wnStart:Nn)-W1(mm)-1*wMu0dt22(mm))/c12(mm);    
            
            %Now calculate one transition probability  between two 
            %grid points W1(mm)  at time t1A and W2(wnStart:Nn) at time t2A
            %by adding the transition probabilities from all three
            %subdivision weighted by their probability mass in each of the
            %smaller three subdivisions.
            pt(mm,wnStart:Nn)=P0(mm)./PSum(mm).*pt0(mm,wnStart:Nn)+ ...
                Pa(mm)./PSum(mm).*pta(mm,wnStart:Nn)+ ...
                Pb(mm)./PSum(mm).*ptb(mm,wnStart:Nn);
            
            
            pt(isnan(pt(mm,:))==1)=0;
            pt(isinf(pt(mm,:))==1)=0;
            
            
            Pnew(wnStart:Nn)=Pnew(wnStart:Nn)+Pprev(mm).*pt(mm,wnStart:Nn);

        end
        
              pRatio(wnStart:Nn)=Pprev(wnStart:Nn)./Pnew(wnStart:Nn);
              
     %   pRatio
     %   ss
     %   str=input('Look at pRatio');
        Pnew(wnStart:Nn)=0.0;
        for mm=1:Nn
            pt(mm,wnStart:Nn)=pRatio(wnStart:Nn).*pt(mm,wnStart:Nn);
            Pnew(wnStart:Nn)=Pnew(wnStart:Nn)+Pprev(mm).*pt(mm,wnStart:Nn);
        %    sum(pt(mm,:))
        %    mm
          %  str=input('look at pt-sum')
        end
        
        
        
        %Pnew(wnStart:Nn)=0.0;
        for mm=1:Nn
            %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)=pt(mm,wnStart:Nn).*(yyMu0dt(mm)+dyyz0(mm).*Zt(mm,wnStart:Nn)+dyy2z0(mm).*(Zt(mm,wnStart:Nn).^2-1));
              dyyGenerator0(mm,wnStart:Nn)=(yyMu0dt(mm)+dyyz0(mm).*Zt(mm,wnStart:Nn)+dyy2z0(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.
            df0Generator(mm,1:Nn)=dfz0(mm).*Zt(mm,1:Nn)+df2z0(mm).*(Zt(mm,1:Nn).^2-1);

            
            %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).*Pprev(mm);
            EZwnew(wnStart:Nn)=EZwnew(wnStart:Nn)+EZwprev(mm).*pt(mm,wnStart:Nn)+ ...
               Zt(mm,wnStart:Nn).*pt(mm,wnStart:Nn).*Pprev(mm);
         
        end
     %    EZwprev(wnStart:Nn)./EZwnew(wnStart:Nn)
     %   plot((wnStart:Nn),EZwnew(wnStart:Nn)./ZProb(wnStart:Nn),'r')
     %   str=input('Look at plot');
  
        
        
        for mm=1:Nn
            if(ss>1)
                Egdtnew(wnStart:Nn)=Egdtnew(wnStart:Nn)+  Egdtprev(mm).*pt(mm,wnStart:Nn);%...
                Egdtnew(wnStart:Nn)=Egdtnew(wnStart:Nn)+ ...
                1*(.5*ds*yy1(mm).^omega1+0*ds*.5*yy2(wnStart:Nn).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm);
            
               Efdtnew(mm)=Efdtnew(mm)+Efdtprev(mm);
             
                Efdtnew(wnStart:Nn)=Efdtnew(wnStart:Nn)+ ...
               1*(ds*.5*yy1(mm).^omega1+ds*.5*yy2(wnStart:Nn).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm);
           
                     
               Ef0new(1:Nn)=Ef0new(1:Nn)+1*Ef0prev(mm).*pt(mm,1:Nn);
                
            end
            if(ss==1)
            Egdtnew(wnStart:Nn)=Egdtnew(wnStart:Nn)+ ...
                (ds*.5*yy1(mm).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm);
             Efdtnew(wnStart:Nn)=Efdtnew(wnStart:Nn)+ ...
                (ds*.5*yy1(mm).^omega1+1*ds*.5*yy2(wnStart:Nn).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm);
            end
  

        end
        %%%%Line 636 starts. Here I have shown how to calculate the
        %%%%multiplier that would transport the functional data specified
        %%%%on one time grid faithfully to next time grid.
        gdtnew(wnStart:Nn)=Egdtnew(wnStart:Nn)./ZProb(wnStart:Nn);
        gdtprev(wnStart:Nn)=Egdtprev(wnStart:Nn)./ZProb(wnStart:Nn);
        
        Cnn(wnStart:Nn)=(gdtprev(wnStart:Nn)+.5*ds*yy1(wnStart:Nn).^omega1)./gdtnew(wnStart:Nn);
        %.5*ds*yy1(wnStart:Nn).^omega1 is added to gdtprev since this value
        %is added at the start of the time interval.
        
      %  Cnn   % Cnn is the multiplier that is used to transport the 
        %function value specified on one grid faithfully to the other grid
        %for function gdtprev.
      %  str=input('Look at Cnn');
        f0new(wnStart:Nn)=Ef0new(wnStart:Nn)./ZProb(wnStart:Nn);
        f0prev(wnStart:Nn)=Ef0prev(wnStart:Nn)./ZProb(wnStart:Nn);
        
        Cf0nn(wnStart:Nn)=(f0prev(wnStart:Nn))./f0new(wnStart:Nn);
        %Cf0nn is the multiplier that is used to transport the 
        %function value specified on one grid faithfully to the other grid
        %for function f0prev.
        %Line 657 Starts
        
    %    plot((wnStart:Nn),Cnn(wnStart:Nn),'g',(wnStart:Nn),Cf0nn(wnStart:Nn),'r')
    %    str=input('Look at correction constants');
        Egdtnew(wnStart:Nn)=0.0;
        Ef0new(wnStart:Nn)=0.0;
                for mm=1:Nn
                    
            if(ss>1)
                Egdtnew(wnStart:Nn)=Egdtnew(wnStart:Nn)+  Egdtprev(mm).*pt(mm,wnStart:Nn).*Cnn(wnStart:Nn).^0;%...
                Egdtnew(wnStart:Nn)=Egdtnew(wnStart:Nn)+  Egdtprev(mm).*pt(mm,wnStart:Nn).*(Cnn(mm)-1);%...
          %      Egdtnew(wnStart:Nn)=Egdtnew(wnStart:Nn)+ ...
          %          1*(ds*yy1(mm).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm).*Cnn(wnStart:Nn);
                Egdtnew(wnStart:Nn)=Egdtnew(wnStart:Nn)+ ...
                    .5*(ds*yy1(mm).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm).*Cnn(wnStart:Nn);
                Egdtnew(wnStart:Nn)=Egdtnew(wnStart:Nn)+ ...
                    .5*(ds*yy2(wnStart:Nn).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm);
                Egdtnew(wnStart:Nn)=Egdtnew(wnStart:Nn)- ...
                (EZwprev(mm)/EZwnew(mm)).*Ef0prev(mm).*pt(mm,1:Nn).*Cf0nn(wnStart:Nn);
           %      Egdtnew(wnStart:Nn)=Egdtnew(wnStart:Nn)+ ...
           %         (-.06-(ss*ds)*16.666*sum(fMu0dt(wnStart:Nn).*ZProb(wnStart:Nn))).*Ef0prev(mm).*pt(mm,1:Nn).*Cf0nn(wnStart:Nn);
                %Cf0nn is multiplier associated with Ef0prev
           %     Egdtnew(wnStart:Nn)=Egdtnew(wnStart:Nn)+ ...
           %         (.012*ss*ds).*df0Generator(mm,1:Nn).*pt(mm,1:Nn).*ZProb(mm);
                 
 
       
               Ef0new(1:Nn)=Ef0new(1:Nn)+1*Ef0prev(mm).*pt(mm,1:Nn).*Cf0nn(wnStart:Nn);
               Ef0new(1:Nn)=Ef0new(1:Nn)+ds*df0Generator(mm,1:Nn).*pt(mm,1:Nn).*ZProb(mm);
            end
            if(ss==1)
            Egdtnew(wnStart:Nn)=Egdtnew(wnStart:Nn)+ ...
                (ds*.5*yy1(mm).^omega1+ds*.5*yy2(wnStart:Nn).^omega1).*pt(mm,wnStart:Nn).*ZProb(mm);
           Ef0new(1:Nn)=Ef0new(1:Nn)+ds*df0Generator(mm,1:Nn).*pt(mm,1:Nn).*ZProb(mm);
            end
                end

 
%        plot((wnStart:Nn),Efdtnew(wnStart:Nn),'r',(wnStart:Nn),Egdtnew(wnStart:Nn),'g',(wnStart:Nn),Ef0prev(wnStart:Nn),'b')
% str=input('Look at graphs');
ss 
    %plot((wnStart:Nn),Efdtnew(wnStart:Nn)./Pnew(wnStart:Nn),'r',(wnStart:Nn),Egdtnew(wnStart:Nn)./Pnew(wnStart:Nn),'g');%,(wnStart:Nn),Ef0prev(wnStart:Nn)./Pnew(wnStart:Nn),'b')
 %str=input('Look at graphs');


 
       Pprev(1:Nn)=Pnew(1:Nn);
       Eyyprev(1:Nn)=Eyynew(1:Nn);    
       EZwprev(1:Nn)=EZwnew(1:Nn);    
       
       
       Efdtprev(1:Nn)=Efdtnew(1:Nn);
       Egdtprev(1:Nn)=Egdtnew(1:Nn);
       
       Ef0prev(1:Nn)=Ef0new(1:Nn);
      
       Ef0new(1:Nn)=0.0;
       
       Pnew(1:Nn)=0;
    %   P1new(1:Nn)=0;
       PZnew(1:Nn)=0;
       Eyynew(1:Nn)=0.0;
       
       EZwnew(1:Nn)=0.0;
       Efdtnew(1:Nn)=0.0;
       Egdtnew(1:Nn)=0.0;

       %Below calculate the start and end time t1A and t2A for the next
       %transition probability time intervals.
       t1A=(ss)*ds;
       t2A=(ss+1)*ds;
        %Below set the starting value W1 for next time interval equal to end
        %value of the current time interval W2 and calculate various
        %integrals. All these integrals only need the starting value W1 
        W1(1:Nn)=W2(1:Nn);

       
      %  [fdtMu0dt,dfdtz01,dfdt2z02]=CalculateFdtDriftAndVol08(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds,t1A,t2A);
     %   [wMu0dt2,dwMu0dt2,dw02]=CalculateDriftAndVol(mu1,mu2,beta1,beta2,gamma,sigma0,W1,1,Nn,ds);
        [yyMu0dt,dyyz0,dyy2z0]=CalculateFunctionDriftAndVol(1.0,1.0,mu1,mu2,beta1,beta2,gamma,sigma0,W1,wnStart,Nn,ds);
        %[fMu0dtdt,dfz0dt,df2z0dt]=CalculateFunctionDriftAndVoldt(theta1,omega1,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,dt);
    end
    
    
    
    
    
    
    %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;
    
    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)+.5*dtM*YY(1:paths).^omega1;
     
     IZM(1:paths)=IZM(1:paths)+Random1(1:paths);
     
 if(rem(tt,8)==0)
     
     
    MaxCutOff=30;
    %NoOfBins=round(300*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too
    NoOfBins=50;
    %[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
%IZM
     
     
    [YDensity,IZMConditional,IndexOutY,IndexMaxY] = MakeDensityFromSimulationConditional_Infiniti(YY,IZM,paths,NoOfBins,MaxCutOff);
     

    MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
    MCVar=sum((YY(:)-MCMean).^2)/paths
    disp('Original process average from our simulation');
    ItoHermiteMean=sum(yy2(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
    ItoHermiteVar=sum((yy2(wnStart+1:Nn-1)-ItoHermiteMean).^2.*ZProb(wnStart+1:Nn-1)) 

    Zw(1:Nn)=EZwprev(1:Nn)./ZProb(1:Nn);
    MuZw=sum(Zw(1:Nn).*ZProb(1:Nn))
    VZw=sum((Zw(1:Nn)-MuZw).^2.*ZProb(1:Nn))
    
    MuIZM=sum(IZM(:))./paths
    VIZM=sum((IZM(:)-MuIZM).^2)/paths
    
    
    ZPerfect(1:Nn)=sqrt(ss)*Z(1:Nn);
    
    plot(yy2(1:Nn),ZPerfect(1:Nn),'b',yy2(1:Nn),Zw(1:Nn),'r',IndexOutY(1:IndexMaxY),IZMConditional(1:IndexMaxY),'g');
    %The above graph plots conditional expected value of integral of 
    %unit normals as a function of y(t) which is the value of SDE variable 
    %in original coordinates. Blue graph is true analytic normal as calculated 
    %from variance, red graph is conditional value of integral calculated 
    %from the transition probabilities grid and green value is conditional 
    %value of integral of normals as a function of y(t) calculated from
    %monte carlo simulations. Please note that all three graphs are
    %initially very close but red graph slowly becomes different and has
    %lesser variance as time increases. While blue graph of true analytic normal
    %calculated from variance and green graph of conditional integral 
    %calculated from monte carlo remain almost exactly the same. Please
    %note that graph of integrated normal will not be a straight line as it is
    %a non-linear function of y(t). Please also note that blue graph of
    %true analytic normal that coincides with monte carlo graph is plotted 
    %conditional of values of y(t) calculated in the analytical 
    %fokker-planck solution of the SDE.
    title('Comparison of Conditional Values of integrated Values of SDE driving unit normals as a function of y(t) SDE variable in original coordinates');
    legend({'True Analytic Conditional Value','Transition Grid Integrated Conditional Value','Monte Carlo Simulation Conditional Integrated Value'},'Location','northwest')
 
    str=input('Look at graph Comparison of Z-integral--first');
    
     plot((1:Nn),ZPerfect(1:Nn),'b',(1:Nn),Zw(1:Nn),'r');
     %Once we have established from previous conditional values graphs 
     %that perfect normal blue graph calculated 
     %from variance is almost identical with monte carlo calculated 
     %integral of normal, we now plot it on an x-axis value for which 
     %normal has to be a linear line so we can compare the true value of 
     %blue graph with red graph of our defective transition probabilities
     %calculated integral of normal in a more familiar way.
    title('Comparison of true analytic value of integrated SDE driving unit normals VS transition probabilities grid integrated value on a linear scale');
    legend({'True Analytic Value','Transition Grid Integrated Value'},'Location','northwest')
 
    str=input('Look at graph Comparison of Z-integral--second');
 end
end

plot((1:Nn),ZProb(1:Nn),'g',(1:Nn),PZprev(1:Nn),'b',(1:Nn),Pprev(1:Nn),'r');
%plot((1:Nn),ZProb(1:Nn),'g',(1:Nn),Pprev(1:Nn),'r');
str=input('Look at distributions');

Pnew=Pprev;
Eyynew=Eyyprev;
Efdtnew=Efdtprev;
Egdtnew=Egdtprev;

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)  ;
yyTr2(wnStart:Nn)=Eyynew(wnStart:Nn)./ZProb(wnStart:Nn)  ;

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

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


%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 D's (the names of variables starting with D) are 
%change of probability derivatives.
%Dfy_w(wmStart:wmEnd)=0;




fydt(wnStart:Nn)=Efdtnew(wnStart:Nn)./ZProb(wnStart:Nn);

gydt(wnStart:Nn)=Egdtnew(wnStart:Nn)./ZProb(wnStart:Nn);
%str=input('We have reached point 1');
%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;
Dgydt(1:Nn)=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));
    Dgydt(mm) = (gydt(mm + 1) - gydt(mm - 1))/(Z(mm + 1) - Z(mm - 1));

end
%str=input('We have reached point 2');
%below pfy and pfydt are respective density amplitudes.
pfy(1:Nn)=0;
pfydt(1:Nn)=0;
pgydt(1:Nn)=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));
    pgydt(mm) = (normpdf(Z(mm),0, 1))/abs(Dgydt(mm));

end



DyyTr(wnStart:Nn)=0.0;
for mm=wnStart+1:Nn-1
    DyyTr(mm) = (yyTr(mm + 1) - yyTr(mm - 1))/(Z(mm + 1) - Z(mm - 1));
end

%below pfy and pfydt are respective density amplitudes.
pyyTr(1:Nn)=0.0;
for mm = wnStart+1:Nn-1
    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



% 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;
%     
%      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)+.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.
MCVar=sum((YY(:)-MCMean).^2)/paths
disp('Original process average from our simulation');
ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
ItoHermiteVar=sum((y_w(wnStart+1:Nn-1)-ItoHermiteMean).^2.*ZProb(wnStart+1:Nn-1)) 
disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average

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
gydtm= sum(gydt(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
gydtvar= sum((gydt(wnStart+1:Nn-1)-gydtm).^2.*ZProb(wnStart+1:Nn-1)) %path integral variance from analytical work.
fydtvar= sum((fydt(wnStart+1:Nn-1)-fydtm).^2.*ZProb(wnStart+1:Nn-1)) %path integral variance from analytical work.

 BinSize=.0075*1/2*2;%Please change this bin size accordingly. When data range is too small decrease it.
 %When density range is too large increase it. Please change it accordingly
 %for a good graph. This is important.
 MaxCutOff=40;
 [fYYdtDensity,IndexOutfYYdt,IndexMaxfYYdt] = MakeDensityFromSimulation_Infiniti(fYYdt,paths,BinSize,MaxCutOff);
 plot(gydt(wnStart+1:Nn-1),pgydt(wnStart+1:Nn-1),'b',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('blue 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');

title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
 
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
end
.
Here is a new function for calculation of conditional distributed values that we will need. It is a very helpful function and I hope friends will like this very small handy function.
.
function [XDensity,YConditional,IndexOut,IndexMax] = MakeDensityFromSimulationConditional_Infiniti(X,Y,Paths,NoOfBins,MaxCutOff)

% 


Xmin=X(1);
Xmax=X(1);
for p=1:Paths
if(X(p)>MaxCutOff)
X(p)=MaxCutOff;
end
%if(X(p)<0)
%X(p)=0;
%end
if(Xmin>real(X(p)))
Xmin=real(X(p));
end
if(Xmax<real(X(p)))
Xmax=real(X(p));
end
end

BinSize=(Xmax-Xmin)/NoOfBins;


IndexMax=floor((Xmax-Xmin)/BinSize+.5)+1
XDensity(1:IndexMax)=0.0;
YConditional(1:IndexMax)=0.0;

for p=1:Paths
index=real(floor(real(X(p)-Xmin)/BinSize+.5)+1);


if(real(index)<1)
index=1;
end
if(real(index)>IndexMax)
index=IndexMax;
end

XDensity(index)=XDensity(index)+1.0/Paths/BinSize;

YConditional(index)=YConditional(index)+Y(p);
%You can calculate many other conditional variables here below
%WConditional(index)=WConditional(index)+W (p)/paths;
%YConditional(index)=YConditional(index)+Y (p)/paths;

end

for nn=1:IndexMax
    if(XDensity(nn)>0)
        YConditional(nn)=YConditional(nn)./XDensity(nn)/Paths./BinSize;
    end
end

IndexOut(1:IndexMax)=Xmin+(0:(IndexMax-1))*BinSize;

end
	

.

When you will run this program, you will get the same output at ss=512 (without comments that I added to explain the meanings of variables) as I have shown in numbers in the previous post at t=512.


.
 
User avatar
Amin
Topic Author
Posts: 2538
Joined: July 14th, 2002, 3:00 am

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

November 6th, 2020, 3:16 am

Friends, I have not been able to work for past few days due to election anxiety. I was a bit nervous yesterday as the lead of my favorite Joe Biden was declining in Arizona. This election means too much to me since it holds a real hope for end of my mind control persecution. My persecution already has lasted more than 23 years and I was extremely worried at the idea that it will continue for another four years or for the rest of my life as some people want. I made this post yesterday in off topic forum:  https://forum.wilmott.com/viewtopic.php?f=15&p=861804#p861804 
I am supposed to see a psychiatrist today and will tell friends in off-topic how it goes. I will also get an antipsychotic injection probably later today and that will likely slow the speed of my new research.
My research has already reached an interesting stage and I hope to continue start working possibly later today and post my programs on the weekend. I really hope that many friends will be interested in my new research and programs as I have some very interesting ideas.
 
User avatar
Amin
Topic Author
Posts: 2538
Joined: July 14th, 2002, 3:00 am

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

November 6th, 2020, 3:21 pm

Please read my post in off-topic
viewtopic.php?f=15&t=94796&p=861849#p861849
 
User avatar
Amin
Topic Author
Posts: 2538
Joined: July 14th, 2002, 3:00 am

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

November 9th, 2020, 2:21 pm

Please read my request for support here and see if you can make a difference:
https://forum.wilmott.com/viewtopic.php?f=15&t=94796&start=1170#p861938
 
User avatar
Amin
Topic Author
Posts: 2538
Joined: July 14th, 2002, 3:00 am

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

November 17th, 2020, 3:09 am

I have made some progress in my research. I am trying to recreate the density of SDE variable in original coordinates using unit normals that are calculated as a difference of two normals with increasing variance but unit difference between them on the fokker-planck grid. The calculation of unit normals on each step is totally independent of dynamics of  the original SDE. I am getting a roughly close match with the original density when unit normals calculated this way are used to generate the SDE variable through the SDE in original coordinates but there are slight discrepancies as there are some inaccuracies in calculation of integrals especially close to the boundary. Now I am trying to make a fokker-planck grid of 100 subdivisions for more accuracy and see if I can reproduce the distribution of SDE variable using the analytically calculated unit normals on the fokker-planck SDE grid. I will come back with results in a few days.
 
User avatar
Amin
Topic Author
Posts: 2538
Joined: July 14th, 2002, 3:00 am

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

November 18th, 2020, 3:06 am

Friends, I want to say that my theory that transition normal between any two points can be solved as a difference of two normal with a unit difference in appropriate variances is probably not true. When I solve for the normal like that and advance the SDE of the SDE variable y(t), I do not get an exact match with true density but I am quite close though. That leads me to think that there is a contribution from second derivative in dynamics that has to be considered when using this normal. Friends would recall from the solution of fokker-planck that derivative with respect to normal(multiplied by normal) takes an extra second order derivative that has to be considered for appropriate solution. I strongly believe that this derivative is also in play when we are considering transition probabilities and their solution and just considering only difference in transition normals would not solve for the equation. Another food for thought and idea that I would also be interested in working is that solution of fokker-planck we considered was along diagonal lines but I believe that a similar form of the equation might possibly hold for transition probabilities when we are starting from one Z(m) on t1 grid and transitioning to a different Z(n) on t2 grid. Again PDE should be valid for all Z's on different time grids and I am not interested in very small peripheral terms(with constants C1, C2 and C3 in my program) that are considered when PDE is solved. Once the PDE has been solved I believe that some relationship should hold in the principal diffusion term that is defined by derivative with respect to normal(multiplied by normal) and an extra second order derivative and added in a squared fashion.

For example when we consider the solution to SDE
[$]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)

We are primarily concerned now in our transition probabilities framework in how to find out a transition probabilities version of the following part of the equation or its modification for two different values of Z's as I believe that effect of other minor terms would have be included when we solved the whole PDE by our original fokker-planck algorithm.
[$]  -  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 Z^2 dt \Big]} [$]

I strongly believe that even in our transition probabilities framework, we have to consider this [$] \frac{\partial^2 X}{\partial Z^2}  [$] term and its effect on the transition. What could be particularly interesting that on the same Z grid, this term would have a slightly different effect for two different diffusions even thought their driving normals would be the same.

I hope to continue to work and try to take this to conclusion in coming days and weeks but friends are also encouraged to experiment. Stochastics, probability theory and finance have a lot of very intelligent people so hopefully we would have a good understanding of the whole thing very soon. 
 
User avatar
Amin
Topic Author
Posts: 2538
Joined: July 14th, 2002, 3:00 am

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

November 21st, 2020, 3:51 pm

Dear friends, I have decided to change my strategy for solution to densities of stochastic integrals a little bit. I have decided to investigate the fokker-planck solution algorithm route simultaneously. Here is the proposed strategy for right now. Obviously things will evolve as the research progresses.
Here are the points of the strategy.
1. Solve the fokker-planck in bessel coordinates first.
2. On the solution grid given by point one, solve for the PDE in original coordinates and functions of the SDE variable in original coordinates. I hope we will learn many interesting things as we solve for the SDE/PDE in original coordinates on already solved bessel grid. Obviously solution algorithm could be briefly written as

[$]f(X(t+1))=\mu(f(X(t)))dt +2 \frac{\partial \nu(f(X(t))) dt}{\partial f}+\sqrt{{\Big [\frac{\partial f}{\partial Z}Z - \frac{\partial^2 f}{\partial Z ^2} \Big]}^2+ {\nu(X)}^2 Z^2 dt }[$]
[$]-\Big [\frac{\partial f}{\partial Z}Z - \frac{\partial^2 f}{\partial Z ^2} \Big]+ Less Significant Terms[$]

where I have omitted less significant terms. here [$]\nu(X)[$] is effective squared volatility. This may be reasonably valid since we have accurately solved our PDE in bessel coordinates and we are solving for original coordinates based on bessel solved grid.

3. From the lessons learnt from the solution of PDE in original coordinates on bessel solved fokker-planck grid, we could try to solve for other functions/integrals (that have explicit dependence on t )of the SDE based on the proposed formula
[$]f(X(t+1),t+1)=\mu(f(X(t),t))dt +2 \frac{\partial \nu(f(X(t),t)) dt}{\partial f}+\sqrt{{\Big [\frac{\partial f}{\partial Z}Z - \frac{\partial^2 f}{\partial Z ^2} \Big]}^2+ {\nu(X,t)}^2 Z^2 dt }[$]
[$]-\Big [\frac{\partial f(X,t)}{\partial Z}Z - \frac{\partial^2 f}{\partial Z ^2} \Big]+ Less Significant Terms[$]

Instead of trying to solve dt-integrals at the start, I have decided to solve the dependent SDEs of the king

[$]dY(t)=t \mu(X) dt + t \sigma(X) dZ(t)[$] 
where X has already been solved using bessel fokker-planck grid. Once we have found solution to above SDE with t, I am sure it will be easier to find a way to solve dt-integrals for their densities leveraging the solution of above SDE with t.

So if we are successful, we will have one bessel solution grid and the dependent function and integrals solution grids superimposed on the same Z-grid.
Once we have solved it, I think it will be easier to convert the solution into transition probabilities grid and try to understand transition probabilities better.
 
User avatar
Amin
Topic Author
Posts: 2538
Joined: July 14th, 2002, 3:00 am

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

November 26th, 2020, 9:40 am

Dear friends, I think I know why the variance of summation of standard normals calculated on SDE transition probabilities grid is lesser than analytic variance. I also have some good sense of why the variance of payoffs/Data specified on various time points on grid continues to slightly decrease with each time step as we move ahead on the transition probabilities grid in time. It would be impossible to properly understand the transition probabilities grid without this new information. I want to do some further experiments and then come back with technical post and the programs in a few days.
Please pressurize the mind control agencies to not actively drug my food since mind control agents  sharply increased activity after I made my last post post on off-topic forum. If things go well and I could properly pursue my research, I will come back in a few days with more results. 
ABOUT WILMOTT

PW by JB

Wilmott.com has been "Serving the Quantitative Finance Community" since 2001. Continued...


Twitter LinkedIn Instagram

JOBS BOARD

JOBS BOARD

Looking for a quant job, risk, algo trading,...? Browse jobs here...


GZIP: On