Serving the Quantitative Finance Community

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

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

April 20th, 2021, 12:57 pm

Friends, I have been working with densities of SDEs that go into zero but I have had limited success when a lot of density goes into zero since both drift and volatility become very unstable especially at volatility power gamma close to .5. However when a small portion of density goes to zero, the whole method works quite well. I have simply turned volatility to first order constant value [$]\sigma \sqrt{dt}[$] close to zero and drift to zero when parameters become unstable but when a lot of density goes into zero, this creates problems with right extreme boundary of the density. I will be posting my programs with graphs either later today or early tomorrow.
However, I am writing this post to let friends know about my future plans. After completing this part of work on SDEs and some of its extensions, I plan to continue for another 3-9 months working on parameter estimation and Bayesian filtering with SDEs and densities in general. Later I want to work on "Statistical information calculus" with densities and Artificial intelligence when information (various variables and related parameters) is given/or inferred in terms of densities that are changing all the time and then artificial intelligence is used to make best sense(risk estimation, projection, prediction and control of desired variables) of all the information we have. I hope that it will continue to become more exciting  and more applied as time passes. I hope to continue for another 5-6 years working with densities and artificial intelligence and will continue posting all my programs here on this forum. I am thinking of some kind of business model so I can make some decent income to live comfortably while I could post all my new research on internet as I have been doing for past several years and I will avoid any activity that will stop me from sharing  my new research and programs with friends here. I will disclose my future research plans in even more detail in coming weeks and would ask friends to wish me well in this endeavor to explore new things together with friends.
 
User avatar
Amin
Topic Author
Posts: 2782
Joined: July 14th, 2002, 3:00 am

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

April 21st, 2021, 7:32 am

Dear friends, I am posting some graphs that I have taken with my program. I will be posting the program in next post.

Image

Image

Image

Image

Image

Image

Image

Image

Image

Image

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

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

April 21st, 2021, 9:04 am

Friends, here is the program. It is a bit messy but I will post a new neat program with explanatory comments in a few hours.
.
function [] = SDETransProb08WmtGrid200b76DriftO0()

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

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

%In this program, I have simulated the density for first few very small time
%intervals using a simple method. After first few intervals, I start
%simulating with probability mass transfer method.  Please experiment with step size and
%you will find the method robust enough. This is just an early version and
%I will be coming with better versions soon.

%dt=.125/16/2/2/2;   % Regular Simulation time interval to advance the density before large step by CDF method.
dt=1/48*1/2*1/32;             %decrease dt for accuracy.

Tt1=32;
Tt2=32;

T=1;
dt2=(T-dt*Tt1)/Tt2;
Tt=Tt1+Tt2;
T=dt*Tt1+dt2*Tt2;
OrderA=4;  %
%dtM=dt/8;%.125/2/4;%Monte carlo time interval size dtM.
%TtM=Tt*8;%T*2*8*4;%Monte carlo number of simulation intervals.
dtM=T/32;%.125/2/4;%Monte carlo time interval size dtM.
TtM=32;


dNn=.05;   % Normal density subdivisions width. would change with number of subdivisions
Nn=200;  % No of normal density subdivisions
NnMidl=100;%One half density Subdivision left from mid of normal density(low)
NnMidh=101;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;

Z(1:Nn)=(((1:Nn)-20.5)*dNn-NnMid)

%There are two indices. One goes to Nn=200. This is used for calculation of
%SDE simulation at next time level t+1. This covers -5 to +5 SDs of
%underlying Z.
%Second index Nn1 below goes from 1 to 240. This covers -6 to +6 SD of
%underlying Z. In probability mass transfer methos, we simulate the SDE
%from -5 SD to +5 SD and then extrapolate on both sides so that we have a
%new grid from -6 SD to +6 SD(1:Nn1). This extrapolated grid(which is assigned time t)
%is used to calculate the time t+1 grid from -5 SD to +5 SD(1:Nn).
%Grid spacing is the same on both grids.
Nn1=240;
%dNn1=.05;
Nn1Midl=120;
Nn1Midh=121;
Nn1Mid=6.025;

Z1(1:Nn1)=((1:Nn1)*dNn-Nn1Mid)
Z1Prob(1)=normcdf(-5.95)-normcdf(-6.0);
Z1Prob(Nn1)=normcdf(6.0)-normcdf(5.95);
Z1Prob(2:Nn1-1)=normcdf(.5*Z1(2:Nn1-1)+.5*Z1(3:Nn1),0,1)-normcdf(.5*Z1(2:Nn1-1)+.5*Z1(1:Nn1-2),0,1);

ZProb(1:Nn)=Z1Prob(21:Nn1-20);
Z(1:Nn)=Z1(21:Nn1-20);

Z1Prob
sum(Z1Prob(1:Nn1))
sum(ZProb(1:Nn))
str=input('Look at Zs');

x0=.01;   % starting value of SDE
beta1=0.0;
beta2=1.00;   % Second drift term power.
gamma=.75;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=1.0;%.950;   %mean reversion parameter.
theta=.02;%mean reversion target
sigma0=.90;%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
Z1
str=input('Look at Zs');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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;


mu33(1)=1;
mu44(1)=1;
mu33(2:OrderA+1)=0;
mu44(2:OrderA+1)=0;

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);
                YqzCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu33(l1).*mu44(l2).*sigma22(l3).*sigma11(l4);
            end
        end
    end
end

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


wnStart=1;%
wnStart1=1;
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
%The above yy(wnStart:Nn)=x0;

tic


for tt=1:Tt 
  
    
   if(tt==1)
   
        w_00=x0^(1-gamma)/(1-gamma);
        %[wMu0dt,c1,c2,c3] = CalculateDriftAndVolA404B4Drift02(w_00,1,1,YqCoeff0,Fp1,gamma,dt);
        [wMu0dt,c1] = CalculateDriftAndVolA404(w_00,1,1,YqCoeff0,Fp1,gamma,dt);
        w(1:Nn)=x0^(1-gamma)/(1-gamma)+wMu0dt(1)+c1(1).*Z(1:Nn);%+c2(1).*(Z(1:Nn).^2-1)+c3(1).*(Z(1:Nn).^3-3*Z(1:Nn));
        
   
    end
     if((tt>1) && (tt<=Tt1))
         %The first Tt1 steps are simulated using a simple and slightly
         %altered version of the method I have used for simple CEV noises
         %and brownian motions. I used it only for starting few time steps
         %since density remains pretty much very close to linear here. If
         %that is not the case, you would have to alter this.
        [wMu0dt,c1] = CalculateDriftAndVolA404(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
         %dw(wnStart:Nn)=c1(wnStart:Nn);%.*Z(wnStart:Nn) ;% ...
        dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn) ;% ...
        dw2(wnStart:Nn)=dw(wnStart:Nn).^2;
        [wMid] = InterpolateOrderN8(8,0,Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidl),Z(NnMidh),Z(NnMidh+1),Z(NnMidh+2),Z(NnMidh+3),w(NnMidl-3),w(NnMidl-2),w(NnMidl-1),w(NnMidl),w(NnMidh),w(NnMidh+1),w(NnMidh+2),w(NnMidh+3));
        B(1:Nn)=w(1:Nn)-wMid;
        %[dBdZ,d2BdZ2,d3BdZ3] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,B,Z);
        B(wnStart:Nn)=sign(B(wnStart:Nn)+dw(wnStart:Nn)).* ...
             sqrt(abs(sign(B(wnStart:Nn)).*(B(wnStart:Nn)).^2+ ...
             sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)));
         
         w(wnStart:Nn)=wMid+wMu0dt(wnStart:Nn)+B(wnStart:Nn);
         
         
         
     end
    %Now starts the calculation for SDE simulation using probability mass
    %transfer.
    if(tt>Tt1)
        %Nn=170;
        %Nn1=190;
         dt=dt2;
       
                
        [wMu0dt,c1] = CalculateDriftAndVolA404(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt); 
        %[wMu0dt,wMu0O0dt,dwMu0dtdw,c1,dc1dw,Sigma2w] = CalculateDriftAndVolA404B4DriftO0(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
        [dwdZ_0,d2wdZ2_0,d3wdZ3_0] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z);
        [dwMu0dtdZ_0,d2wMu0dtdZ2_0,d3wMu0dtdZ3_0] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,wMu0dt,Z);
        %BB1(Nn+20+1:Nn+20+20)=BB1(Nn+20)+(1:20).*dwdZ_0(Nn)*dNn;      
        
        if(wnStart>=1)
        wnStart1=wnStart+20;
        end
        BB1(20+wnStart:Nn+20)=w(wnStart:Nn);
        if(wnStart1>=21)
            %BB1(1:20)=BB1(21)-(20:-1:1).*dwdZ_0(1)*dNn;
            BB1(wnStart1-20:wnStart1-1)=BB1(wnStart1)-(20:-1:1).*dwdZ_0(1)*dNn;
            %BB1(1:wnStart1-1)=BB1(wnStart1+10)-(wnStart1+10-1:-1:11).*dwdZ_0(wnStart+10)*dNn;
        
            wnStart1=wnStart1-20;
            LoopFlag=0;
            for mm=1:20
                if((BB1(mm)>=0)&&(LoopFlag==0))
                    wnStart1=mm;
                    LoopFlag=1;
                end
            end
           % for mm=Nnmid1:-1:wnStart1
           %     if(BB1(mm)>BB1(mm+1))
           %         wnStart1=mm+1;
           %         BB1(1:mm1)=0.0;
           %     end
           % end
            for mm=wnStart1:Nn1Midl
                if(BB1(mm)>=BB1(mm+1))
                    wnStart1=mm+1;
                    BB1(1:mm)=0.0;
                end
            end
        end
        BB1(Nn-10+20+1:Nn+20+20)=BB1(Nn-10+20)+(1:30).*dwdZ_0(Nn-10)*dNn;      

        
        %wnStart=1;
            ZZ1=Z1;
            NN1=Nn1;
            ZZProb(1:NN1)=Z1Prob(1:Nn1);
            if(wnStart1>1)
                ZZProb(wnStart1)=Z1Prob(wnStart1)+sum(Z1Prob(1:wnStart1-1));
                ZZProb(1:wnStart1-1)=0.0;
            end
            
     
            %plot(ZZ1(1:NN1),BB1(1:NN1),'b',ZZ1(1:NN1),BB1a(1:NN1),'y',ZZ1(1:NN1),BB1b(1:NN1),'k',Z(1:Nn),B1(1:Nn),'r')
            %str=input('Look at B1 and BB1, original and extrapolated version and boundaries');
            
            %Below calculate first four derivatives. Please replace with your favorite
            %robust function. 
      %      [dBdZ,d2BdZ2,d3BdZ3] = First3Derivatives2ndOrderEqSpacedA(1,NN1,dNn,BB1,Z1);
      %      [dwdZAA,d2wdZ2AA,d3wdZ3AA,d4BdZ4] = First4Derivatives2ndOrderEqSpaced(1,NN1,dNn,BB1,Z1);
        
            %plot((1:NN1),dBdZ(1:NN1),'r')
            %str=input('Look at 1st derivative');
            %plot((1:NN1),d2BdZ2(1:NN1),'r')
            %str=input('Look at 2nd derivative');
            %plot((1:NN1),d3BdZ3(1:NN1),'r')
            %str=input('Look at 3rd derivative');
        
            
            ww1(wnStart1:NN1)=BB1(wnStart1:NN1);%+wMid;
            ww1(ww1<=0)=0.000000000001;
            %I will improve the above improvisation in a day.
            %SigmawA below is volatility associated with each grid cell at time
            %t Bm(t) or BB1.
         %   [wMu0dtA,dwMu0dtdwA,d2wMu0dtdw2A,d3wMu0dtdw3A,d4wMu0dtdw4A,d5wMu0dtdw5A,d6wMu0dtdw6A,d7wMu0dtdw7A,d8wMu0dtdw8A,SigmawA,dSigmawdwA,d2Sigmawdw2A,Sigma2wA,Sigma3wA] = CalculateDriftAndVolA404B4Drift(ww1,wnStart,NN1,YqCoeff0,Fp1,gamma,dt);
 
            [Muw,MuwO0,dMuwdw,Sigmaw,dSigmawdw,Sigma2w] = CalculateDriftAndVolA404B4DriftO0(ww1,wnStart1,NN1,YqCoeff0,Fp1,gamma,dt);

    %        plot(ZZ1(wnStart1:NN1),ww1(wnStart1:NN1),'b',Z(wnStart1:Nn),w(wnStart1:Nn),'r')
    %        tt
    %        str=input('Look at comparison of interpolated and original w-0');

             %plot(ZZ1(1:NN1),Sigmaw(1:NN1),'b',Z(1:Nn),c1(1:Nn),'r')
    %         plot(ZZ1(1:NN1),Sigmaw(1:NN1),'b',ZZ1(1:NN1),dSigmawdw(1:NN1),'r')
    %        str=input('Look at comparison of original and interpolated vol coefficient-1');
            %plot(ZZ1(1:NN1),cc1(1:NN1),'b',Z(1:Nn),c1(1:Nn),'r')
            %str=input('Look at comparison of original and interpolated vol coefficient-1A');
            %plot(ZZ1(1:NN1),Muw(1:NN1),'b',Z(1:Nn),wMu0dt(1:Nn),'r')
     %       plot(ZZ1(1:NN1),Muw(1:NN1),'b',ZZ1(1:NN1),dMuwdw(1:NN1),'r')
     %       str=input('Look at drifts origianl and interpolated-2');
            %plot(ZZ1(1:NN1),ZZProb(1:NN1),'b',Z(1:Nn),ZProb(1:Nn),'r')
            %str=input('Look at Zprobss -3');
            %plot((1:Nn),Z(1:Nn),'r',(1:NN1)-20,ZZ1(1:NN1),'b')
            %str=input('Look at Zs -4');
       
            %Below are calculations of CDF and its derivatives to calculate
            %the mid point of the grid at time t+1. This mid point should
            %lie exactly on median which means CDF there would be exactly
            %equal to .5 or CDF=.5.
            
            c1Guess(wnStart:Nn)=sigma0*sqrt(dt);
            %Sigmaw_0(wnStart:Nn)=sqrt(dwdZ_0(wnStart:Nn).^2+ c1(wnStart:Nn).^2);
            Sigmaw_0(wnStart:Nn)=sqrt(dwdZ_0(wnStart:Nn).^2+ c1Guess(wnStart:Nn).^2);
            
           
%            plot(Z(wnStart:Nn),dwdZ_0(wnStart:Nn),'r');
%            str=input('Look at dwdZ_0 before addition to Bn');
            
%            plot(Z(wnStart:Nn),wMu0dt(wnStart:Nn),'r');
%            str=input('Look at dwdZ_0 before addition to Bn');
            
            
            for nn=wnStart:Nn
                if(w(nn)<=0)
                  wMu0dt(nn)=0.0;
                end
            end
            
            if(nn<Nn)
          if(5*abs(wMu0dt(nn)-wMu0dt(nn+1))<abs(dwMu0dtdw(nn)-dwMu0dtdw(nn+1)))
           %   wMu0dt(nn)=0.0;
          end
          %if(5*abs(Sigmaw(mm)-Sigmaw(mm+1))<abs(dSigmawdw(mm)-dSigmawdw(mm+1)))
          %    Sigmaw(mm)=sigma0*sqrt(dt);
          %end
            end
          wMu0dtA=wMu0dt;
            if(gamma>=.75)
                    %wMu0dt(wnStart:wnStart+30)=wMu0dt(wnStart+30)-dwMu0dtdw(wnStart+30)*(w(wnStart+30)-w(w(wnStart:-1:wnStart-30)));
                    wMu0dtA(wnStart:wnStart+9)=wMu0dt(wnStart+10)-(10:-1:1).*dNn.*dwMu0dtdZ_0(wnStart+10);
                    %wMu0dt(wnStart:wnStart+9)=wMu0dt(wnStart+10)-dwMu0dtdw(wnStart+10)*-1*(w(wnStart:wnStart+4)-w(wnStart+5));
            end
                
           plot((1:Nn),wMu0dt(1:Nn),'r',(1:Nn),wMu0dtA(1:Nn),'b');
           str=input('Look at wMu0dt linearly interpolated');
                
            
            
            
            
            
            
            Bn(wnStart:Nn)=w(wnStart:Nn)+wMu0dtA(wnStart:Nn)+(Sigmaw_0(wnStart:Nn)-dwdZ_0(wnStart:Nn)).*Z(wnStart:Nn);
            Bn(Bn<=0)=0.0;%000000000001;
            
            plot((wnStart:Nn),Bn(wnStart:Nn),'r',(wnStart:Nn),w(wnStart:Nn),'g')
            str=input('Look at Bn and w');
            
            
            
            wnStart
%            plot((1:Nn),Bn(1:Nn),'r');
%            str=input('Look at graphs Bn-2--before');
            for nn=NnMidl:-1:wnStart+1
                if(Bn(nn)<Bn(nn-1))
                    %Bn(nn-1)=Bn(nn);
                    Bn(1:nn-1)=Bn(nn);
                    wnStart=nn;
                end
            end
%            plot((1:Nn),Bn(1:Nn),'r');
%            str=input('Look at graphs Bn-2--After');
%            plot((1:Nn),Z(1:Nn),'r');
%            str=input('Look at graphs Z--1');
            
            
            
            
            
            %for nn=NnMid:-1:2
            %    if(Bn(nn)<Bn(nn-1))
            %        Bn(nn-1)=Bn(nn);
            %    end
            %end
            %for nn=1:Nn
            %    if(Bn(nn)<=0)
            %        Bn(nn)=0.0;
            %    end
            %end
            
            
%            plot(Z(wnStart:Nn),Bn(wnStart:Nn),'r');
%            str=input('Look at Bn');
            
            
            %I will improve the above improvisation in a day.
            
            Bm0(wnStart1:NN1)=BB1(wnStart1:NN1);%
      %      Bm1(1:NN1)=BB1a(1:NN1);%
      %      Bm2(1:NN1)=BB1b(1:NN1);%
            Bm0(Bm0<0)=0.0;%00000000001;
            for mm=Nn1Midl:-1:wnStart1+1
                if(Bm0(mm)<Bm0(mm-1))
                    %Bn(nn-1)=Bn(nn);
                    Bm0(1:mm-1)=Bm0(mm);
                    wnStart1=mm;
                end
            end
            
            
%            plot(ZZ1(wnStart1:NN1),Bm0(wnStart1:NN1),'r');
%            str=input('Look at Bm');
            
       %     plot((1:Nn),c1(1:Nn),'r');
       %     str=input('Look at Bn plot---0');
       %     plot((1:Nn),SigmaB(1:Nn),'r',(1:Nn),dBdZ_0(1:Nn),'b');
       %     str=input('Look at Bn plot---1');
       %     plot((1:Nn),Bn(1:Nn));
       %     str=input('Look at Bn plot');
            
           Pn(1:Nn)=0.0;
           pn(1:Nn)=0.0;
           ESecondHeTotal(1:Nn)=0.0;
           %EThirdHeTotal(1:Nn)=0.0;
           
           
           %for mm0=wnStart1:NN1
           %     if(Bm0(mm0)<=0)
           %             Muw(mm)=0.0;
           %             Sigmaw(mm)=sigma0*sqrt(dt);
           %             Sigma2w(mm)=0.0;
           %             
           %     end
           % end
       %    SumThreshld=0.0;
       %     yy(wnStart1:NN1)=((1-gamma)*Bm0(wnStart1:NN1)).^(1/(1-gamma));
       %     for mm=wnStart1:NN1
       %         Threshld(mm)=floor(gamma*sigma0*yy(mm)^(gamma-1)+1);
       %         if(Threshld(mm)>1)
       %             SumThreshld=SumThreshld+Threshld(mm);
       %         end
       %     end
                
            
            ZZProb(1:NN1)=Z1Prob(1:Nn1);
            if(wnStart1>1)
                ZZProb(wnStart1)=Z1Prob(wnStart1)+sum(Z1Prob(1:wnStart1-1));
                ZZProb(1:wnStart1-1)=0.0;
            end
            
            

            for mm=wnStart1:NN1

                
                
                
                %if(Bm0(mm)<=0)
                %        Muw(mm)=0.0;
                %        Sigmaw(mm)=sigma0*sqrt(dt);
                %        Sigma2w(mm)=0.0;
                %        str=input('I am inside the loop you are looking for');
                %end
                %mm
                
               % MuwMedian=Muw(Nn1Midh);
               % if(Muw(mm)>1.5*MuwMedian)
               %     Threshld=mm-wnStart1;
               % end
                
                Threshld=0.0;
                yy(mm)=((1-gamma)*Bm0(mm)).^(1/(1-gamma));
                if(sigma0*(yy(mm)^(gamma-1))>2)
                  %  Threshld=mm-wnStart1;
                    %Threshld=floor(gamma*sigma0*yy(mm)^(gamma-1)+1)-1;
                    ThreshldT=mm-wnStart1;
                    %Threshld=floor(.375.^0*2*gamma.^2*sigma0.^2*yy(mm)^(2*gamma-2)+1);
                    if(ThreshldT>Threshld)
                        Threshld=ThreshldT;
                    end
                end
      %          Threshld
                %Bm0(mm)
                %yy(mm)
                
                
                if(wnStart1<40)
                Threshld=floor(wnStart1/(25+1));
                %Threshld=floor(wnStart1/(15+1));
                %Threshld=floor(wnStart1/(10+1));
                
                end
                if((wnStart1>40)&&(wnStart1<95))
                %Threshld=floor(wnStart1/18+1);
                Threshld=floor(wnStart1/(25+1));
                Threshld=floor(wnStart1/20+1);
 
                end

       %         if(Threshld(mm)>=1)
       %             Muw(mm)=(1-Threshld(mm)/SumThreshld)*Muw(mm)   ;
       %              Sigmaw(mm)=sigma0*sqrt(dt);
       %              Sigma2w(mm)=0.0;
 
       %         end

       
                %frac=(mm-wnStart1)/Threshld;
                bound=wnStart1+Threshld;
                keata=log(.125);%./Threshld;
              
                 if(mm<=wnStart1+Threshld)
                     %mm
                     %bound-mm
                     %    (exp(keata*(bound-mm)))
                         Muw(mm)=0.0;%MuwO0(mm);%Muw(mm).*(exp(keata*(bound-mm)));%Muw(mm)./(1+Threshld);
                         Sigmaw(mm)=sigma0*sqrt(dt);
                         Sigma2w(mm)=0.0;
  %                        str=input('I am inside the loop you are looking for');
%                 elseif(yy(mm)<.005)
               %      Muw(mm)=MuwO0(mm);
               %      Sigmaw(mm)=sigma0*sqrt(dt);
                 end

          %    if(yy(mm)<.0175)
          %      Muw(mm)=0;
          %      Sigmaw(mm)=sigma0*sqrt(dt);
          %      Sigma2w(mm)=0.0;
          %    end
          if(mm<NnMidl)
          if(mm<=wnStart1+6)    
          %if((4*abs(Muw(mm)-Muw(mm+1))<abs(dMuwdw(mm)-dMuwdw(mm+1)))||(4*abs(Sigmaw(mm)-Sigmaw(mm+1))<abs(dSigmawdw(mm)-dSigmawdw(mm+1))))
 %             Muw(mm)=0.0;
      %    end
      %    if(3*abs(Sigmaw(mm)-Sigmaw(mm+1))<abs(dSigmawdw(mm)-dSigmawdw(mm+1)))
 %             Sigmaw(mm)=sigma0*sqrt(dt);
          end
          end
                
           %     if((mm<=wnStart1+Threshld)&&(mm>=wnStart1+.5*Threshld))
           %             Muw(mm)=Muw(mm)./(2+Threshld);
           %             Sigmaw(mm)=sigma0*sqrt(dt);
           %             Sigma2w(mm)=0.0;
 %         %               str=input('I am inside the loop you are looking for');
           %     end
                
             %    if((mm<=wnStart1+Threshld)&&(mm>=wnStart1+.75*Threshld))
             %           Muw(mm)=Muw(mm)./(1+Threshld);
             %           Sigmaw(mm)=sigma0*sqrt(dt);
             %           Sigma2w(mm)=0.0;
 %           %             str=input('I am inside the loop you are looking for');
             %   end
                
                %if(Bm0(mm)<=Bm0(wnStart1)+.01)
                %       nnExp=nn
                %end
                %if(Bm0(mm)<=Bm0(wnStart1)+.01)
                %if((mm)<=(wnStart1)+10)    
                %        Muw(mm)=Muw(mm).*exp(-.35*(10-(mm-wnStart1)));
                %        Sigmaw(mm)=sigma0*sqrt(dt);
                %        Sigma2w(mm)=0.0;
%               %         str=input('I am inside the loop you are looking for');
                %end

                
                nn1=wnStart;
                for nn=wnStart:Nn
                    if( (Bn(nn)-Bm0(mm)-Muw(mm)) <-5*abs(Sigmaw(mm)))
                        nn1=nn;
                    end
                end

                nn2=Nn;
                for nn=Nn:-1:wnStart
                    if( (Bn(nn)-Bm0(mm)-Muw(mm)) >5*abs(Sigmaw(mm)))
                        nn2=nn;
                    end
                end
                Zt(1:Nn)=0.0;
                Pmn(1:Nn)=0.0;
                ESecondHe(1:Nn)=0.0;
%Bt(mm1:mm2)=Bn(mm1:mm2)-Bm0;
                Zt(nn1:nn2)=(Bn(nn1:nn2)-Bm0(mm)-Muw(mm))/Sigmaw(mm);
                Zt(nn1:nn2)=real(Zt(nn1:nn2));
                Pm0n(nn1:nn2)=normcdf(Zt(nn1:nn2),0,1);
                pm0n(nn1:nn2)=normpdf(Zt(nn1:nn2),0,1);
                %ESecondHe(nn1:nn2)=(Sigma2w(mm).*(Zt(nn1:nn2).^2-1)).*pmn(nn1:nn2);
                
                Pmn(nn1:nn2)=Pm0n(nn1:nn2).*ZZProb(mm);
                pmn(nn1:nn2)=pm0n(nn1:nn2).*ZZProb(mm);
                %[Pmn,pmn,ESecondHe] = CalculateCDFArraySignIndexDrift(Bn,Nn,wMu0dtA(mm),dwMu0dtdwA(mm),d2wMu0dtdw2A(mm),d3wMu0dtdw3A(mm),d4wMu0dtdw4A(mm),d5wMu0dtdw5A(mm),d6wMu0dtdw6A(mm),d7wMu0dtdw7A(mm),d8wMu0dtdw8A(mm),SigmawA(mm),dSigmawdwA(mm),d2Sigmawdw2A(mm),Sigma2wA(mm),Bm0(mm),Bm1(mm),Bm2(mm),ZZProb(mm),pm0(mm),dpm0(mm),d2pm0(mm),d3pm0(mm),d4pm0(mm),d5pm0(mm),d6pm0(mm),d7pm0(mm));
                
                ESecondHe(nn1:nn2)=(Sigma2w(mm).*(Zt(nn1:nn2).^2-1)).*pmn(nn1:nn2);
                
                
                
                if(nn1>wnStart)
                    Pmn(1:nn1-1)=0;
                    pmn(1:nn1-1)=0;
                    ESecondHe(1:nn1-1)=0;
                  %  EThirdHe(1:nn1-1)=0;
                end
                if (nn2<Nn)
                    Pmn(nn2+1:Nn)=ZZProb(mm);
                    pmn(nn2+1:Nn)=0;
                    ESecondHe(nn2+1:Nn)=0;
                 %   EThirdHe(nn2+1:Nn)=0;

                end
                
                
                %Pn is total CDF at the point Bn while Pmn is contribution
                %to CDF from mth grid cell at previous time point.
                Pmn(isnan(Pmn)==1)=0.0;
                pmn(isnan(pmn)==1)=0.0;
                ESecondHe(isnan(ESecondHe)==1)=0.0;
                
                Pn(wnStart:Nn)=Pn(wnStart:Nn)+Pmn(wnStart:Nn);
                pn(wnStart:Nn)=pn(wnStart:Nn)+pmn(wnStart:Nn);
                ESecondHeTotal(wnStart:Nn)=ESecondHeTotal(wnStart:Nn)+ESecondHe(wnStart:Nn);
                
                
            end
            
%Calculate Znew which is Znew-grid that corresponds to our Guess Bn where CDF is now available.                
            Znew(wnStart:Nn)=norminv(Pn(wnStart:Nn));
            
            %Bn2(wnStart:Nn)=Bn(wnStart:Nn)+ESecondHeTotal(wnStart:Nn)./pn(wnStart:Nn);%
            %for nn=wnStart+1:Nn
            Bn2(wnStart:Nn)=Bn(wnStart:Nn);%+ESecondHeTotal(wnStart:Nn)./pn(wnStart:Nn);%
            001
            wnStart
            for nn=1:NnMid
                if(Bn2(nn)<0)
                    wnStart=nn+1;
                end
            end
            002
            wnStart
            
            %Bn2(1:Nn)=Bn(1:Nn)+ESecondHeTotal(1:Nn);%./pn(1:Nn);%
     %       wnStart
     %       Nn
%            plot((wnStart:Nn),Znew(wnStart:Nn),'r');
%            str=input('Look at Bn graphs');
%            plot((wnStart:Nn),Bn(wnStart:Nn),'g',(wnStart:Nn),Bn2(wnStart:Nn),'r');
%            str=input('Look at Bn graphs');

%Below we Interpolate values of w at Z-grid (on variable Z)from our 
%calculated values of Znew that correspond to CDF at Bn.      
       %[w] = InterpolateNewGrid(Bn,Znew,wnStart,Nn,dNn,Z)
       [w] = InterpolateNewGrid(Bn2,Znew,wnStart,Nn,dNn,Z);
        if(wnStart>5)
            [w(wnStart-1)] = InterpolateOrderN6(6,Z(wnStart-1),Znew(wnStart),Znew(wnStart+1),Znew(wnStart+2),Znew(wnStart+3),Znew(wnStart+4),Znew(wnStart+5),Bn2(wnStart),Bn2(wnStart+1),Bn2(wnStart+2),Bn2(wnStart+3),Bn2(wnStart+4),Bn2(wnStart+5));
            [w(wnStart-2)] = InterpolateOrderN6(6,Z(wnStart-2),Znew(wnStart),Znew(wnStart+1),Znew(wnStart+2),Znew(wnStart+3),Znew(wnStart+4),Znew(wnStart+5),Bn2(wnStart),Bn2(wnStart+1),Bn2(wnStart+2),Bn2(wnStart+3),Bn2(wnStart+4),Bn2(wnStart+5));
            %[w(wnStart-3)] = InterpolateOrderN6(6,Z(wnStart-3),Znew(wnStart),Znew(wnStart+1),Znew(wnStart+2),Znew(wnStart+3),Znew(wnStart+4),Znew(wnStart+5),Bn2(wnStart),Bn2(wnStart+1),Bn2(wnStart+2),Bn2(wnStart+3),Bn2(wnStart+4),Bn2(wnStart+5));
            %[w(wnStart-4)] = InterpolateOrderN6(6,Z(wnStart-4),Znew(wnStart),Znew(wnStart+1),Znew(wnStart+2),Znew(wnStart+3),Znew(wnStart+4),Znew(wnStart+5),Bn2(wnStart),Bn2(wnStart+1),Bn2(wnStart+2),Bn2(wnStart+3),Bn2(wnStart+4),Bn2(wnStart+5));
            %[w(wnStart-5)] = InterpolateOrderN6(6,Z(wnStart-5),Znew(wnStart),Znew(wnStart+1),Znew(wnStart+2),Znew(wnStart+3),Znew(wnStart+4),Znew(wnStart+5),Bn2(wnStart),Bn2(wnStart+1),Bn2(wnStart+2),Bn2(wnStart+3),Bn2(wnStart+4),Bn2(wnStart+5));
        wnStart=wnStart-2;
        end
       
    %   if(wnStart>5)
    %       [w(wnStart-1)] = InterpolateOrderN6(6,Z(wnStart-1),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));
    %       [w(wnStart-2)] = InterpolateOrderN6(6,Z(wnStart-2),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));
    %       [w(wnStart-3)] = InterpolateOrderN6(6,Z(wnStart-3),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));
    %       [w(wnStart-4)] = InterpolateOrderN6(6,Z(wnStart-4),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));
    %       [w(wnStart-5)] = InterpolateOrderN6(6,Z(wnStart-5),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));
    %   wnStart=wnStart-5;
    %   end
       
       
       
       
       wnStartA=wnStart
       Nn
       1
       
       for nn=1:NnMid
                if(w(nn)<0)
                    wnStart=nn;
                end
            end
            
       wnStartB=wnStart
       Nn
       2
       if(wnStartA~=wnStartB)
           str=input('Different values of wnStart detected');
       end
       
%       plot((wnStart:Nn),Znew(wnStart:Nn),'g',(wnStart:Nn),Z(wnStart:Nn),'r');
%       str=input('Look at old and interpolated Zs');
     
       w(w<0)=0.0;
       
        plot((1:Nn),Bn(1:Nn),'g',(1:Nn),w(1:Nn),'r');
            str=input('Look at w-Bn graphs after iterpolation');
       
       
       
       [dwdZ_1,d2wdZ2_1,d3wdZ3_1] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z);
       
%       plot(Z(1:Nn),dwdZ_1(1:Nn),'r');
%       str=input('Look at dwdZ graph');             
       %plot((1:Nn),Bn(1:Nn),'g',(1:Nn),Bn2(1:Nn),'r');
       %     str=input('Look at Bn graphs');

       
       if(wnStart1<=35)
       w(Nn-20+1:Nn)=w(Nn-20)+(1:20).*dwdZ_1(Nn-20)*dNn;
       end
       
  %     if(wnStart1>35)
       w(Nn-30+1:Nn)=w(Nn-30)+(1:30).*dwdZ_1(Nn-30)*dNn;
  %     end
       % w(Nn-40+1:Nn)=w(Nn-40)+(1:40).*dwdZ_1(Nn-40)*dNn;      
%      plot((1:Nn),w(1:Nn),'r');
%     str=input('Look at w graph');

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

  %      w(isnan(w)==1)=0.0;
        
  %      w(w<0)=0.0;
  %      for nn=wnStart:Nn
  %          if(w(nn)<=0)
  %             % wnStart=nn+1;
  %              w(1:wnStart-1)=0.0;
  %          end
  %      end

%Below is the grsph of newly calculated value of w.       
%       plot((1:Nn),w(1:Nn),'r')
%       str=input('Calculated Value of w');
    end
end
    


%below D's (the names of variables starting with D) are 
%change of probability derivatives.

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

rng(29079137, 'twister')
paths=200000;
YY(1:paths)=x0;  %Original process monte carlo.
Random1(1:paths)=0;
%for tt=1:(Tt-1)+8
for tt=1:TtM    
if(tt>Tt1)
    %dtM=dt2/8;
end
    Random1=randn(size(Random1));
    HermiteP1(1,1:paths)=1;
    HermiteP1(2,1:paths)=Random1(1:paths);
    HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
    HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
    HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;
    

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

end

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

MaxCutOff=30;
NoOfBins=round(300*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too 
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );

%plot(y_w1(wnStart+1:Nn-1),py_w1(wnStart+1:Nn-1),'b',y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
plot(y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');

title(sprintf('x0 = %.4f,theta=%.4f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.3f,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));
 
legend({'New Model Density','Monte Carlo Density'},'Location','northeast')
 
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
end
.
One sub function
.
function [wMu0dt,c1] = CalculateDriftAndVolA404(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;

wMu0dtOdt(wnStart:Nn)=(YqCoeff0(1,1,2,1).*yy(wnStart:Nn).^Fp1(1,1,2,1)*0+ ...
    YqCoeff0(1,2,1,1).*yy(wnStart:Nn).^Fp1(1,2,1,1)+ ...
    YqCoeff0(2,1,1,1).*yy(wnStart:Nn).^Fp1(2,1,1,1))*dt ;

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



end
.
Second sub function.
.
function [wMu0dt,wMu0dtO0,dwMu0dtdw,c1,dc1dw,c2] = CalculateDriftAndVolA404B4DriftO0(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt)

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

Fp2=Fp1/(1-gamma);
wMu0dtO0(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;

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,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 + ...
     (YqCoeff0(1,1,5,1).*Fp1(1,1,5,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,1,5,1))+ ...
      YqCoeff0(1,2,4,1).*Fp1(1,2,4,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,2,4,1))+ ...
      YqCoeff0(2,1,4,1).*Fp1(2,1,4,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,1,4,1))+ ...
       YqCoeff0(2,2,3,1).*Fp1(2,2,3,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,2,3,1))+ ...
       YqCoeff0(3,1,3,1).*Fp1(3,1,3,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(3,1,3,1))+ ...
       YqCoeff0(1,4,2,1).*Fp1(1,4,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,4,2,1))+ ...
       YqCoeff0(2,3,2,1).*Fp1(2,3,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,3,2,1))+ ...
       YqCoeff0(3,2,2,1).*Fp1(3,2,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(3,2,2,1))+ ...
       YqCoeff0(4,1,2,1).*Fp1(4,1,2,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(4,1,2,1))+ ...
       YqCoeff0(1,5,1,1).*Fp1(1,5,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,5,1,1))+ ...
       YqCoeff0(2,4,1,1).*Fp1(2,4,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,4,1,1))+ ...
       YqCoeff0(3,3,1,1).*Fp1(3,3,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(3,3,1,1))+  ...
       YqCoeff0(4,2,1,1).*Fp1(4,2,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(4,2,1,1))+ ...
       YqCoeff0(5,1,1,1).*Fp1(5,1,1,1).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(5,1,1,1)))*dt^4;
 
   
   
   
   
   
   
   

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



 dc1dw(wnStart:Nn)=((YqCoeff0(1,1,1,2).*Fp1(1,1,1,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,1,1,2)).*sqrt(dt))+ ...
     (YqCoeff0(1,1,2,2).*Fp1(1,1,2,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,1,2,2))+ ...
     YqCoeff0(1,2,1,2).*Fp1(1,2,1,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,2,1,2))+ ...
     YqCoeff0(2,1,1,2).*Fp1(2,1,1,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,1,1,2))).*dt^1.5+ ...
     (YqCoeff0(1,1,3,2).*Fp1(1,1,3,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,1,3,2))+ ...
     YqCoeff0(1,2,2,2).*Fp1(1,2,2,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,2,2,2))+ ...
     YqCoeff0(2,1,2,2).*Fp1(2,1,2,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,1,2,2))+ ...
     YqCoeff0(1,3,1,2).*Fp1(1,3,1,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,3,1,2))+ ...
     YqCoeff0(2,2,1,2).*Fp1(2,2,1,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,2,1,2))+ ...
     YqCoeff0(3,1,1,2).*Fp1(3,1,1,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(3,1,1,2))).*dt^2.5+ ...
     (YqCoeff0(1,1,4,2).*Fp1(1,1,4,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,1,4,2))+ ...
     YqCoeff0(1,2,3,2).*Fp1(1,2,3,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,2,3,2))+ ...
     YqCoeff0(2,1,3,2).*Fp1(2,1,3,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,1,3,2))+ ...
     YqCoeff0(1,3,2,2).*Fp1(1,3,2,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,3,2,2))+ ...
     YqCoeff0(2,2,2,2).*Fp1(2,2,2,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,2,2,2))+ ...
     YqCoeff0(3,1,2,2).*Fp1(3,1,2,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(3,1,2,2))+ ...
     YqCoeff0(1,4,1,2).*Fp1(1,4,1,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(1,4,1,2))+ ...
     YqCoeff0(2,3,1,2).*Fp1(2,3,1,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(2,3,1,2))+ ...
     YqCoeff0(3,2,1,2).*Fp1(3,2,1,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(3,2,1,2))+ ...
     YqCoeff0(4,1,1,2).*Fp1(4,1,1,2).*((1-gamma).*w(wnStart:Nn)).^(-1+Fp2(4,1,1,2))).*dt^3.5);



 
 c2(wnStart:Nn)=((YqCoeff0(1,1,1,3).*yy(wnStart:Nn).^Fp1(1,1,1,3) *dt) + ...
     (YqCoeff0(1,1,2,3).*yy(wnStart:Nn).^Fp1(1,1,2,3)+ ...
     YqCoeff0(1,2,1,3).*yy(wnStart:Nn).^Fp1(1,2,1,3)+ ...
     YqCoeff0(2,1,1,3).*yy(wnStart:Nn).^Fp1(2,1,1,3)).*dt^2+ ...
     (YqCoeff0(1,1,3,3).*yy(wnStart:Nn).^Fp1(1,1,3,3)+ ...
     YqCoeff0(1,2,2,3).*yy(wnStart:Nn).^Fp1(1,2,2,3)+ ...
     YqCoeff0(2,1,2,3).*yy(wnStart:Nn).^Fp1(2,1,2,3) + ...
     YqCoeff0(1,3,1,3).*yy(wnStart:Nn).^Fp1(1,3,1,3)+ ...
     YqCoeff0(2,2,1,3).*yy(wnStart:Nn).^Fp1(2,2,1,3)+ ...
     YqCoeff0(3,1,1,3).*yy(wnStart:Nn).^Fp1(3,1,1,3)).*dt^3);
 
 c3(wnStart:Nn)=((YqCoeff0(1,1,1,4).*yy(wnStart:Nn).^Fp1(1,1,1,4)*dt^1.5 )+ ...
     (YqCoeff0(1,1,2,4).*yy(wnStart:Nn).^Fp1(1,1,2,4)+ ...
     YqCoeff0(1,2,1,4).*yy(wnStart:Nn).^Fp1(1,2,1,4)+ ...
     YqCoeff0(2,1,1,4).*yy(wnStart:Nn).^Fp1(2,1,1,4))*dt^2.5);
 

end
.
.
.Here is the output of the program with given hardcoded parameters.
.
ItoHermiteMean =
   0.016292515154324
true Mean only applicable to standard SV mean reverting type models otherwise disregard
TrueMean =
   0.016321205588286


Original process average from monte carlo
MCMean =
  0.016422305493483 - 0.000000075516855i
MCVar =
  0.001154848355913 - 0.000000008907504i
Original process average from our simulation
ItoHermiteMean =
   0.016292515154324
ItoHermiteVar =
   0.001114551566303
true Mean only applicble to standard SV mean reverting type models otherwise disregard
TrueMean =
   0.016321205588286
IndexMax =
        2371
red line is density of SDE from Ito-Hermite method, green is monte carlo.>> 

Here is the output graph.
Image


Please note that I very slightly modified the program after making graphs in previous post. Values there should be closely but not identically retrieved. The above graph is exactly from the same program posted.
I will be posting another neat program with comments in a few hours(same thing but mess deleted and comments added).


Friends, I quickly noticed that first graph in previous post could not be produced by this program and it blew up. 
Please change in the body of the main function
 
         wMu0dtA=wMu0dt;
            if(gamma>=.75)
                    %wMu0dt(wnStart:wnStart+30)=wMu0dt(wnStart+30)-dwMu0dtdw(wnStart+30)*(w(wnStart+30)-w(w(wnStart:-1:wnStart-30)));
                    wMu0dtA(wnStart:wnStart+9)=wMu0dt(wnStart+10)-(10:-1:1).*dNn.*dwMu0dtdZ_0(wnStart+10);
                    %wMu0dt(wnStart:wnStart+9)=wMu0dt(wnStart+10)-dwMu0dtdw(wnStart+10)*-1*(w(wnStart:wnStart+4)-w(wnStart+5));
            end
                
           plot((1:Nn),wMu0dt(1:Nn),'r',(1:Nn),wMu0dtA(1:Nn),'b');
           str=input('Look at wMu0dt linearly interpolated');


to 


 
         wMu0dtA=wMu0dt;
         %   if(gamma>=.75)
         %         %wMu0dt(wnStart:wnStart+30)=wMu0dt(wnStart+30)-dwMu0dtdw(wnStart+30)*(w(wnStart+30)-w(w(wnStart:-1:wnStart-30)));
         %           wMu0dtA(wnStart:wnStart+9)=wMu0dt(wnStart+10)-(10:-1:1).*dNn.*dwMu0dtdZ_0(wnStart+10);
          %          %wMu0dt(wnStart:wnStart+9)=wMu0dt(wnStart+10)-dwMu0dtdw(wnStart+10)*-1*(w(wnStart:wnStart+4)-w(wnStart+5));
           % end
                
           plot((1:Nn),wMu0dt(1:Nn),'r',(1:Nn),wMu0dtA(1:Nn),'b');
           str=input('Look at wMu0dt linearly interpolated');


and the above program will work perfectly. Just comment out the above if loop and program will work perfectly. Sorry about this little problem.
 
User avatar
Amin
Topic Author
Posts: 2782
Joined: July 14th, 2002, 3:00 am

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

April 21st, 2021, 12:41 pm

Friends, here is the neat and commented version of the main program. I have commented only the part that corresponds to main new algorithm while comments at older section of the program are only legacy. I am sure friends will easily understand the program.
Here is the program.
function [] = SDETransProb08WmtGrid200b80DriftO0()

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

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

%In this program, I have simulated the density for first few very small time
%intervals using a simple method. After first few intervals, I start
%simulating with probability mass transfer method.  Please experiment with step size and
%you will find the method robust enough. This is just an early version and
%I will be coming with better versions soon.

%dt=.125/16/2/2/2;   % Regular Simulation time interval to advance the density before large step by CDF method.
dt=1/48*1/2*1/32;             %decrease dt for accuracy.

Tt1=32;
Tt2=32;

T=1;
dt2=(T-dt*Tt1)/Tt2;
Tt=Tt1+Tt2;
T=dt*Tt1+dt2*Tt2;
OrderA=4;  %
%dtM=dt/8;%.125/2/4;%Monte carlo time interval size dtM.
%TtM=Tt*8;%T*2*8*4;%Monte carlo number of simulation intervals.
dtM=T/32;%.125/2/4;%Monte carlo time interval size dtM.
TtM=32;


dNn=.05;   % Normal density subdivisions width. would change with number of subdivisions
Nn=200;  % No of normal density subdivisions
NnMidl=100;%One half density Subdivision left from mid of normal density(low)
NnMidh=101;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;

Z(1:Nn)=(((1:Nn)-20.5)*dNn-NnMid)

%There are two indices. One goes to Nn=200. This is used for calculation of
%SDE simulation at next time level t+1. This covers -5 to +5 SDs of
%underlying Z.
%Second index Nn1 below goes from 1 to 240. This covers -6 to +6 SD of
%underlying Z. In probability mass transfer methos, we simulate the SDE
%from -5 SD to +5 SD and then extrapolate on both sides so that we have a
%new grid from -6 SD to +6 SD(1:Nn1). This extrapolated grid(which is assigned time t)
%is used to calculate the time t+1 grid from -5 SD to +5 SD(1:Nn).
%Grid spacing is the same on both grids.
Nn1=240;
%dNn1=.05;
Nn1Midl=120;
Nn1Midh=121;
Nn1Mid=6.025;

Z1(1:Nn1)=((1:Nn1)*dNn-Nn1Mid)
Z1Prob(1)=normcdf(-5.95)-normcdf(-6.0);
Z1Prob(Nn1)=normcdf(6.0)-normcdf(5.95);
Z1Prob(2:Nn1-1)=normcdf(.5*Z1(2:Nn1-1)+.5*Z1(3:Nn1),0,1)-normcdf(.5*Z1(2:Nn1-1)+.5*Z1(1:Nn1-2),0,1);

ZProb(1:Nn)=Z1Prob(21:Nn1-20);
Z(1:Nn)=Z1(21:Nn1-20);

Z1Prob
sum(Z1Prob(1:Nn1))
sum(ZProb(1:Nn))
str=input('Look at Zs');

x0=.2;   % starting value of SDE
beta1=0.0;
beta2=1.00;   % Second drift term power.
gamma=.75;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=1.0;%.950;   %mean reversion parameter.
theta=.01;%mean reversion target
sigma0=1.0;%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
Z1
str=input('Look at Zs');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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;


mu33(1)=1;
mu44(1)=1;
mu33(2:OrderA+1)=0;
mu44(2:OrderA+1)=0;

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);
                YqzCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu33(l1).*mu44(l2).*sigma22(l3).*sigma11(l4);
            end
        end
    end
end

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


wnStart=1;%
wnStart1=1;
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
%The above yy(wnStart:Nn)=x0;
NN1=Nn1;
tic


for tt=1:Tt 
  
    
   if(tt==1)
   
        w_00=x0^(1-gamma)/(1-gamma);
        %[wMu0dt,c1,c2,c3] = CalculateDriftAndVolA404B4Drift02(w_00,1,1,YqCoeff0,Fp1,gamma,dt);
        [wMu0dt,c1] = CalculateDriftAndVolA404(w_00,1,1,YqCoeff0,Fp1,gamma,dt);
        w(1:Nn)=x0^(1-gamma)/(1-gamma)+wMu0dt(1)+c1(1).*Z(1:Nn);%+c2(1).*(Z(1:Nn).^2-1)+c3(1).*(Z(1:Nn).^3-3*Z(1:Nn));
        
   
    end
     if((tt>1) && (tt<=Tt1))
         %The first Tt1 steps are simulated using a simple and slightly
         %altered version of the method I have used for simple CEV noises
         %and brownian motions. I used it only for starting few time steps
         %since density remains pretty much very close to linear here. If
         %that is not the case, you would have to alter this.
        [wMu0dt,c1] = CalculateDriftAndVolA404(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
         %dw(wnStart:Nn)=c1(wnStart:Nn);%.*Z(wnStart:Nn) ;% ...
        dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn) ;% ...
        dw2(wnStart:Nn)=dw(wnStart:Nn).^2;
        [wMid] = InterpolateOrderN8(8,0,Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidl),Z(NnMidh),Z(NnMidh+1),Z(NnMidh+2),Z(NnMidh+3),w(NnMidl-3),w(NnMidl-2),w(NnMidl-1),w(NnMidl),w(NnMidh),w(NnMidh+1),w(NnMidh+2),w(NnMidh+3));
        B(1:Nn)=w(1:Nn)-wMid;
        %[dBdZ,d2BdZ2,d3BdZ3] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,B,Z);
        B(wnStart:Nn)=sign(B(wnStart:Nn)+dw(wnStart:Nn)).* ...
             sqrt(abs(sign(B(wnStart:Nn)).*(B(wnStart:Nn)).^2+ ...
             sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)));
         
         w(wnStart:Nn)=wMid+wMu0dt(wnStart:Nn)+B(wnStart:Nn);
         
         
         
     end
     %Below there are two grids. The starting grid and target grid. Probability 
     %at starting grid(grid at time: t) is used with SDE to predict the probability mass in 
     %target grid(grid at time: t+1). Target
     %grid is w/Bn/Bn2 while starting grid is BB1/Bm0. Starting grid is 
     %larger while target grid is smaller with the same grid cell spacing.
     %Starting index for starting grid is wnStart1 and starting index for
     %target grid is wnStart. End index for starting grid is Nn1/NN1 while
     %that of target grid is Nn. Index variable for starting grid is mm
     %while that of target grid is nn.
     %At every time point target grid w from previous time is used to
     %create a larger grid by extrapolation and this larger grid
     %will be used as starting grid to calculate the target grid using the
     %SDE analytics.
     %When target grid becomes unstable, the first grid cell that remains
     %good and stable is assigned to wnStart.
     %Similalry first good and stable cell for the starting grid is
     %wnStart1.
     %At start starting grid and target grid are so constructed that
     %starting grid has 20 extra cells padded on both boundaries of target grid.
     %S0 w(nn)=BB1(nn+20) because 20 points have been padded on both sides
     %of w(target grid which is also Bn) to construct the target grid BB1(or Bm0).
     if(tt>Tt1)    %Tt1 is when calculation with new method starts.

         dt=dt2;   %assign an appropriate step size to dt in this part of simulation.
       
                
        [wMu0dt,c1] = CalculateDriftAndVolA404(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt); 
        [dwdZ_0,d2wdZ2_0,d3wdZ3_0] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z);
        
        if(wnStart>=1)
            wnStart1=wnStart+20;%Because of 20 padded cells at start.
        end
        BB1(20+wnStart:Nn+20)=w(wnStart:Nn);  %Target grid from previous 
        %step  is assigned to construct the starting grid. Index is changed
        %by amount of padding=20 on the left boundary od starting grid.
        if(wnStart1>=21)
            BB1(wnStart1-20:wnStart1-1)=BB1(wnStart1)-(20:-1:1).*dwdZ_0(wnStart)*dNn;
 
            %Above 20 cells are padded at starting boundary of starting
            %grid using linear extrapolation from the derivative on cell at
            %the boundary of the target grid from previous point.
            %"dwdZ_0(wnStart)" is derivative on first point on the left
            %boundary of previous grid.
                
            wnStart1=wnStart1-20;
            %Above:Now wnStart1 is appropriately decreased by 20 points after
            %padding 20 points using linear extrapolation.
   
            
            %Below: BB1 should strictly increase as we increase the index.
            %If some value is smaller than value to the left, it means that
            %grid has become corrupt there so we simply set that point as
            %starting point of our calculations and assign it as wnStart1.
            for mm=wnStart1:Nn1Midl
                if(BB1(mm)>=BB1(mm+1))
                    wnStart1=mm+1;
                    BB1(1:mm)=0.0;
                end
            end
        end
        BB1(Nn-10+20+1:Nn+20+20)=BB1(Nn-10+20)+(1:30).*dwdZ_0(Nn-10)*dNn;      

        %Above extrapolated values using derivative at ten points to the 
        %left of right boundary are padded to starting grid. Since
        %sometimes right side also becomes somewhat corrupt due to CDF
        %miscalculations on the right, I add thirty values to the right
        %using a value from ten points deep inside the left boundary. 20
        %values are to be padded to the right and ten extra points are
        %extrapolated in case they have gone bad and a total of thirty
        %points are being extrapolated.
        %wnStart=1;
           
            
            
     
            %plot(ZZ1(1:NN1),BB1(1:NN1),'b',ZZ1(1:NN1),BB1a(1:NN1),'y',ZZ1(1:NN1),BB1b(1:NN1),'k',Z(1:Nn),B1(1:Nn),'r')
            %str=input('Look at B1 and BB1, original and extrapolated version and boundaries');
            
            ww1(wnStart1:NN1)=BB1(wnStart1:NN1);
 
            [Muw,Sigmaw] = CalculateDriftAndVolA404(ww1,wnStart1,NN1,YqCoeff0,Fp1,gamma,dt);

    %        plot(ZZ1(wnStart1:NN1),ww1(wnStart1:NN1),'b',Z(wnStart1:Nn),w(wnStart1:Nn),'r')
            tt
    %        str=input('Look at comparison of interpolated and original w-0');

             %plot(ZZ1(1:NN1),Sigmaw(1:NN1),'b',Z(1:Nn),c1(1:Nn),'r')
    %         plot(ZZ1(1:NN1),Sigmaw(1:NN1),'b',ZZ1(1:NN1),dSigmawdw(1:NN1),'r')
    %        str=input('Look at comparison of original and interpolated vol coefficient-1');
            %plot(ZZ1(1:NN1),cc1(1:NN1),'b',Z(1:Nn),c1(1:Nn),'r')
            %str=input('Look at comparison of original and interpolated vol coefficient-1A');
            %plot(ZZ1(1:NN1),Muw(1:NN1),'b',Z(1:Nn),wMu0dt(1:Nn),'r')
     %       plot(ZZ1(1:NN1),Muw(1:NN1),'b',ZZ1(1:NN1),dMuwdw(1:NN1),'r')
     %       str=input('Look at drifts origianl and interpolated-2');
            %plot(ZZ1(1:NN1),ZZProb(1:NN1),'b',Z(1:Nn),ZProb(1:Nn),'r')
            %str=input('Look at Zprobss -3');
            %plot((1:Nn),Z(1:Nn),'r',(1:NN1)-20,ZZ1(1:NN1),'b')
            %str=input('Look at Zs -4');
       
            
            c1Guess(wnStart:Nn)=sigma0*sqrt(dt);
            Sigmaw_0(wnStart:Nn)=sqrt(dwdZ_0(wnStart:Nn).^2+ c1Guess(wnStart:Nn).^2);
            
           
%            plot(Z(wnStart:Nn),dwdZ_0(wnStart:Nn),'r');
%            str=input('Look at dwdZ_0 before addition to Bn');
            
%            plot(Z(wnStart:Nn),wMu0dt(wnStart:Nn),'r');
%            str=input('Look at dwdZ_0 before addition to Bn');
            

%Below: we project our (guess) target grid using estimate of volatility and
%drift which are being added to target grid at previous point.
            
            Bn(wnStart:Nn)=w(wnStart:Nn)+wMu0dt(wnStart:Nn)+(Sigmaw_0(wnStart:Nn)-dwdZ_0(wnStart:Nn)).*Z(wnStart:Nn);
            Bn(Bn<=0)=0.0;%000000000001;
            
      %      plot((wnStart:Nn),Bn(wnStart:Nn),'r',(wnStart:Nn),w(wnStart:Nn),'g')
      %      str=input('Look at Bn and w');
            
            
            
            wnStart
            
            %Below: After calculating target grid, we recalculate wnStart
            %based on the criteria that as we move to right on the grid,
            %value of grid has to increase and if decreases somewhere, our
            %grid has become corrupt there and we make that point as
            %starting point of caclulations and assign it wnStart.
            
            for nn=NnMidl:-1:wnStart+1
                if(Bn(nn)<Bn(nn-1))
                    %Bn(nn-1)=Bn(nn);
                    Bn(1:nn-1)=Bn(nn);
                    wnStart=nn;
                end
            end
%            plot((1:Nn),Bn(1:Nn),'r');
%            str=input('Look at graphs Bn-2--After');
%            plot((1:Nn),Z(1:Nn),'r');
%            str=input('Look at graphs Z--1');
            
            
            
            
            
%            plot(Z(wnStart:Nn),Bn(wnStart:Nn),'r');
%            str=input('Look at Bn');
            
%Below we have the starting grid that is used as starting point for 
%probability calculations.            
            Bm0(wnStart1:NN1)=BB1(wnStart1:NN1);%
   
            Bm0(Bm0<0)=0.0;%00000000001;

            %Below: Again on starting grid, we recalculate wnStart1
            %based on the criteria that as we move to right on the grid,
            %value of grid has to increase and if decreases somewhere, our
            %grid has become corrupt there and we make that point as
            %starting point of caclulations and assign it wnStart1.
            
            for mm=Nn1Midl:-1:wnStart1+1
                if(Bm0(mm)<Bm0(mm-1))
                    Bm0(1:mm-1)=Bm0(mm);
                    wnStart1=mm;
                end
            end
            
            
%            plot(ZZ1(wnStart1:NN1),Bm0(wnStart1:NN1),'r');
%            str=input('Look at Bm');
            
       %     plot((1:Nn),c1(1:Nn),'r');
       %     str=input('Look at Bn plot---0');
       %     plot((1:Nn),SigmaB(1:Nn),'r',(1:Nn),dBdZ_0(1:Nn),'b');
       %     str=input('Look at Bn plot---1');
       %     plot((1:Nn),Bn(1:Nn));
       %     str=input('Look at Bn plot');
            
           Pn(1:Nn)=0.0;
            
            ZZProb(1:NN1)=Z1Prob(1:Nn1);
            
       %Below we assign all the probability mass to the left of the left 
       %boundary of starting grid to the first cell with index wnStart1 that 
       %is used to make the calculations. This way we make sure that lost
       %probability(that is not on the calculation grid) is accounted for 
       %in the calculations of CDF
            if(wnStart1>1)
                ZZProb(wnStart1)=Z1Prob(wnStart1)+sum(Z1Prob(1:wnStart1-1));
                ZZProb(1:wnStart1-1)=0.0;
            end
            
            

            for mm=wnStart1:NN1

                %Below threshold is calculated using number of points on
                %the padded starting grid that are extra on top of wnstart
                %of the target grid of the previous time. These points were
                %padded using linear extrapolation and now we want to use them 
                %with zero drift and constant first order volatility.
                Threshold=wnStart-wnStart1+20-1;
              
                %Below drift is turned to zero in the linearly extrapolated
                %padded zone and volatility is turned to first order value.
                 if(mm<=wnStart1+Threshold)
                         Muw(mm)=0.0;
                         Sigmaw(mm)=sigma0*sqrt(dt);
                 end

                
                nn1=wnStart;
                for nn=wnStart:Nn
                    if( (Bn(nn)-Bm0(mm)-Muw(mm)) <-5*abs(Sigmaw(mm)))
                        nn1=nn;
                    end
                end

                nn2=Nn;
                for nn=Nn:-1:wnStart
                    if( (Bn(nn)-Bm0(mm)-Muw(mm)) >5*abs(Sigmaw(mm)))
                        nn2=nn;
                    end
                end
                Zt(1:Nn)=0.0;
                Pmn(1:Nn)=0.0;

                Zt(nn1:nn2)=(Bn(nn1:nn2)-Bm0(mm)-Muw(mm))/Sigmaw(mm);
                %Zt(nn1:nn2)=real(Zt(nn1:nn2));
                Pm0n(nn1:nn2)=normcdf(Zt(nn1:nn2),0,1);
                
                Pmn(nn1:nn2)=Pm0n(nn1:nn2).*ZZProb(mm);
                
                
                if(nn1>wnStart)
                    Pmn(1:nn1-1)=0;
                end
                if (nn2<Nn)
                    Pmn(nn2+1:Nn)=ZZProb(mm);
                
                end
                
                
                %Pn is total CDF at the point Bn while Pmn is contribution
                %to CDF from mth grid cell at previous time point.
                Pmn(isnan(Pmn)==1)=0.0;
                
                Pn(wnStart:Nn)=Pn(wnStart:Nn)+Pmn(wnStart:Nn);
                
                
            end
            
%Calculate Znew which is Znew-grid that corresponds to our Guess Bn where CDF is now available.                
            Znew(wnStart:Nn)=norminv(Pn(wnStart:Nn));
            wnStart
            for nn=1:NnMid
                if(Bn(nn)<0)
                    wnStart=nn+1;
                end
            end
       [w] = InterpolateNewGrid(Bn,Znew,wnStart,Nn,dNn,Z);
       
       
       wnStartA=wnStart
      % Nn
      % 1
       
       for nn=1:NnMid
                if(w(nn)<0)
                    wnStart=nn;
                end
            end
            
       wnStartB=wnStart
      % Nn
      % 2
      % if(wnStartA~=wnStartB)
      %     str=input('Different values of wnStart detected');
      % end
       
%       plot((wnStart:Nn),Znew(wnStart:Nn),'g',(wnStart:Nn),Z(wnStart:Nn),'r');
%       str=input('Look at old and interpolated Zs');
     
       w(w<0)=0.0;
       
 %       plot((1:Nn),Bn(1:Nn),'g',(1:Nn),w(1:Nn),'r');
 %           str=input('Look at w-Bn graphs after iterpolation');
       
       
       
       [dwdZ_1,d2wdZ2_1,d3wdZ3_1] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z);
       
%       plot(Z(1:Nn),dwdZ_1(1:Nn),'r');
%       str=input('Look at dwdZ graph');             
       %plot((1:Nn),Bn(1:Nn),'g',(1:Nn),Bn2(1:Nn),'r');
       %     str=input('Look at Bn graphs');

       
       
       %Below: We linearly re-calculate by linear extrapolation 
       %the last 20/30 cells on the grid 
       %next to the right boundary since sometimes rarely due to CDF
       %miscalculations, and when these points are recalculated, the grid 
       %remains stable. This is usually not needed but when the process
       %comes under stress due to miscalculations on the left side, this
       %helps to keep everything in sync.
       
       
       if(wnStart1<=35)
       w(Nn-20+1:Nn)=w(Nn-20)+(1:20).*dwdZ_1(Nn-20)*dNn;
       end
       
  %     if(wnStart1>35)
       w(Nn-30+1:Nn)=w(Nn-30)+(1:30).*dwdZ_1(Nn-30)*dNn;
  %     end
       
       
 
%Below is the grsph of newly calculated value of w.       
%       plot((1:Nn),w(1:Nn),'r')
%       str=input('Calculated Value of w');
    end
end
    


%below D's (the names of variables starting with D) are 
%change of probability derivatives.

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

rng(29079137, 'twister')
paths=200000;
YY(1:paths)=x0;  %Original process monte carlo.
Random1(1:paths)=0;
%for tt=1:(Tt-1)+8
for tt=1:TtM    
if(tt>Tt1)
    %dtM=dt2/8;
end
    Random1=randn(size(Random1));
    HermiteP1(1,1:paths)=1;
    HermiteP1(2,1:paths)=Random1(1:paths);
    HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
    HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
    HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;
    

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

end

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

MaxCutOff=30;
NoOfBins=round(300*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too 
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );

%plot(y_w1(wnStart+1:Nn-1),py_w1(wnStart+1:Nn-1),'b',y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
plot(y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');

title(sprintf('x0 = %.4f,theta=%.4f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.3f,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));
 
legend({'New Model Density','Monte Carlo Density'},'Location','northeast')
 
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
end
.
.
Here is the last output when you run this program:

Elapsed time is 2.476761 seconds.
ItoHermiteMean =
   0.080058929576259
true Mean only applicable to standard SV mean reverting type models otherwise disregard
TrueMean =
   0.079897093822574
Original process average from monte carlo
MCMean =
  0.080205149759109 - 0.000000084956820i
MCVar =
  0.025910900616490 + 0.000000011733038i
Original process average from our simulation
ItoHermiteMean =
   0.080058929576259
ItoHermiteVar =
   0.025206289664637
true Mean only applicble to standard SV mean reverting type models otherwise disregard
TrueMean =
   0.079897093822574
IndexMax =
        1193
red line is density of SDE from Ito-Hermite method, green is monte carlo.

And here is the graph, you will see when you run the program. I have adjusted the Axes scales to visualize the main body of the graph. 

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

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

April 23rd, 2021, 7:31 am

Friends, I want to alter the above program and try to add second hermite polynomial in calculations using quadratic formula(and not earlier approximation). I want to see whether adding second hermite polynomial improves the density on the right boundary where it sometimes becomes unstable since in some of my other experiments I have seen that adding second hermite polynomial adds to stability on the right. I am not sure if it would work but want to explore this and will let friends know how it goes.
 
User avatar
Amin
Topic Author
Posts: 2782
Joined: July 14th, 2002, 3:00 am

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

April 23rd, 2021, 9:56 am

Here I am posting the modified version of the previous program in which I have used quadratic formula with second hermite polynomial for simulation of SDEs. I have removed(commented) the linear reconstruction on the right boundary of the density with linear extrapolation and it seems to work perfectly well(without any extrapolation reconstruction) for a lot of cases I ran with simulations but you can still uncomment it if you need it somewhere. After the final time of the simulation, a lot of points on the left starting boundary of the density have been cut, I reconstruct them with linear extrapolation that I did not do before. You can replace it with your favorite method of extrapolation.
Here is the program.

function [] = SDETransProb08WmtGrid200b80DriftO0()

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

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

%In this program, I have simulated the density for first few very small time
%intervals using a simple method. After first few intervals, I start
%simulating with probability mass transfer method.  Please experiment with step size and
%you will find the method robust enough. This is just an early version and
%I will be coming with better versions soon.

%dt=.125/16/2/2/2;   % Regular Simulation time interval to advance the density before large step by CDF method.
dt=1/48*1/2*1/32;             %decrease dt for accuracy.

Tt1=32;
Tt2=32;

T=1;
dt2=(T-dt*Tt1)/Tt2;
Tt=Tt1+Tt2;
T=dt*Tt1+dt2*Tt2;
OrderA=4;  %
%dtM=dt/8;%.125/2/4;%Monte carlo time interval size dtM.
%TtM=Tt*8;%T*2*8*4;%Monte carlo number of simulation intervals.
dtM=T/32;%.125/2/4;%Monte carlo time interval size dtM.
TtM=32;


dNn=.05;   % Normal density subdivisions width. would change with number of subdivisions
Nn=200;  % No of normal density subdivisions
NnMidl=100;%One half density Subdivision left from mid of normal density(low)
NnMidh=101;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;

Z(1:Nn)=(((1:Nn)-20.5)*dNn-NnMid)

%There are two indices. One goes to Nn=200. This is used for calculation of
%SDE simulation at next time level t+1. This covers -5 to +5 SDs of
%underlying Z.
%Second index Nn1 below goes from 1 to 240. This covers -6 to +6 SD of
%underlying Z. In probability mass transfer methos, we simulate the SDE
%from -5 SD to +5 SD and then extrapolate on both sides so that we have a
%new grid from -6 SD to +6 SD(1:Nn1). This extrapolated grid(which is assigned time t)
%is used to calculate the time t+1 grid from -5 SD to +5 SD(1:Nn).
%Grid spacing is the same on both grids.
Nn1=240;
%dNn1=.05;
Nn1Midl=120;
Nn1Midh=121;
Nn1Mid=6.025;

Z1(1:Nn1)=((1:Nn1)*dNn-Nn1Mid)
Z1Prob(1)=normcdf(-5.95)-normcdf(-6.0);
Z1Prob(Nn1)=normcdf(6.0)-normcdf(5.95);
Z1Prob(2:Nn1-1)=normcdf(.5*Z1(2:Nn1-1)+.5*Z1(3:Nn1),0,1)-normcdf(.5*Z1(2:Nn1-1)+.5*Z1(1:Nn1-2),0,1);

ZProb(1:Nn)=Z1Prob(21:Nn1-20);
Z(1:Nn)=Z1(21:Nn1-20);

Z1Prob
sum(Z1Prob(1:Nn1))
sum(ZProb(1:Nn))
str=input('Look at Zs');

x0=.014;   % starting value of SDE
beta1=0.0;
beta2=1.00;   % Second drift term power.
gamma=.85;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=1.0;%.950;   %mean reversion parameter.
theta=.014;%mean reversion target
sigma0=1.0;%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
Z1
str=input('Look at Zs');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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;


mu33(1)=1;
mu44(1)=1;
mu33(2:OrderA+1)=0;
mu44(2:OrderA+1)=0;

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);
                YqzCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu33(l1).*mu44(l2).*sigma22(l3).*sigma11(l4);
            end
        end
    end
end

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


wnStart=1;%
wnStart1=1;
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
%The above yy(wnStart:Nn)=x0;
NN1=Nn1;
tic


for tt=1:Tt 
  
    
   if(tt==1)
   
        w_00=x0^(1-gamma)/(1-gamma);
        %[wMu0dt,c1,c2,c3] = CalculateDriftAndVolA404B4Drift02(w_00,1,1,YqCoeff0,Fp1,gamma,dt);
        [wMu0dt,c1] = CalculateDriftAndVolA404(w_00,1,1,YqCoeff0,Fp1,gamma,dt);
        w(1:Nn)=x0^(1-gamma)/(1-gamma)+wMu0dt(1)+c1(1).*Z(1:Nn);%+c2(1).*(Z(1:Nn).^2-1)+c3(1).*(Z(1:Nn).^3-3*Z(1:Nn));
        
   
    end
     if((tt>1) && (tt<=Tt1))
         %The first Tt1 steps are simulated using a simple and slightly
         %altered version of the method I have used for simple CEV noises
         %and brownian motions. I used it only for starting few time steps
         %since density remains pretty much very close to linear here. If
         %that is not the case, you would have to alter this.
        [wMu0dt,c1] = CalculateDriftAndVolA404(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
         %dw(wnStart:Nn)=c1(wnStart:Nn);%.*Z(wnStart:Nn) ;% ...
        dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn) ;% ...
        dw2(wnStart:Nn)=dw(wnStart:Nn).^2;
        [wMid] = InterpolateOrderN8(8,0,Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidl),Z(NnMidh),Z(NnMidh+1),Z(NnMidh+2),Z(NnMidh+3),w(NnMidl-3),w(NnMidl-2),w(NnMidl-1),w(NnMidl),w(NnMidh),w(NnMidh+1),w(NnMidh+2),w(NnMidh+3));
        B(1:Nn)=w(1:Nn)-wMid;
        %[dBdZ,d2BdZ2,d3BdZ3] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,B,Z);
        B(wnStart:Nn)=sign(B(wnStart:Nn)+dw(wnStart:Nn)).* ...
             sqrt(abs(sign(B(wnStart:Nn)).*(B(wnStart:Nn)).^2+ ...
             sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)));
         
         w(wnStart:Nn)=wMid+wMu0dt(wnStart:Nn)+B(wnStart:Nn);
         
         
         
     end
     %Below there are two grids. The starting grid and target grid. Probability 
     %at starting grid(grid at time: t) is used with SDE to predict the probability mass in 
     %target grid(grid at time: t+1). Target
     %grid is w/Bn/Bn2 while starting grid is BB1/Bm0. Starting grid is 
     %larger while target grid is smaller with the same grid cell spacing.
     %Starting index for starting grid is wnStart1 and starting index for
     %target grid is wnStart. End index for starting grid is Nn1/NN1 while
     %that of target grid is Nn. Index variable for starting grid is mm
     %while that of target grid is nn.
     %At every time point target grid w from previous time is used to
     %create a larger grid by extrapolation and this larger grid
     %will be used as starting grid to calculate the target grid using the
     %SDE analytics.
     %When target grid becomes unstable, the first grid cell that remains
     %good and stable is assigned to wnStart.
     %Similalry first good and stable cell for the starting grid is
     %wnStart1.
     %At start starting grid and target grid are so constructed that
     %starting grid has 20 extra cells padded on both boundaries of target grid.
     %S0 w(nn)=BB1(nn+20) because 20 points have been padded on both sides
     %of w(target grid which is also Bn) to construct the target grid BB1(or Bm0).
     if(tt>Tt1)    %Tt1 is when calculation with new method starts.

         dt=dt2;   %assign an appropriate step size to dt in this part of simulation.
       
                
        [wMu0dt,c1] = CalculateDriftAndVolA404(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt); 
        [dwdZ_0,d2wdZ2_0,d3wdZ3_0] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z);
        
        if(wnStart>=1)
            wnStart1=wnStart+20;%Because of 20 padded cells at start.
        end
        BB1(20+wnStart:Nn+20)=w(wnStart:Nn);  %Target grid from previous 
        %step  is assigned to construct the starting grid. Index is changed
        %by amount of padding=20 on the left boundary od starting grid.
        if(wnStart1>=21)
            BB1(wnStart1-20:wnStart1-1)=BB1(wnStart1)-(20:-1:1).*dwdZ_0(wnStart)*dNn;
 
            %Above 20 cells are padded at starting boundary of starting
            %grid using linear extrapolation from the derivative on cell at
            %the boundary of the target grid from previous point.
            %"dwdZ_0(wnStart)" is derivative on first point on the left
            %boundary of previous grid.
                
            wnStart1=wnStart1-20;
            %Above:Now wnStart1 is appropriately decreased by 20 points after
            %padding 20 points using linear extrapolation.
   
            
            %Below: BB1 should strictly increase as we increase the index.
            %If some value is smaller than value to the left, it means that
            %grid has become corrupt there so we simply set that point as
            %starting point of our calculations and assign it as wnStart1.
            for mm=wnStart1:Nn1Midl
                if(BB1(mm)>=BB1(mm+1))
                    wnStart1=mm+1;
                    BB1(1:mm)=0.0;
                end
            end
        end
        BB1(Nn-10+20+1:Nn+20+20)=BB1(Nn-10+20)+(1:30).*dwdZ_0(Nn-10)*dNn;      

        %Above extrapolated values using derivative at ten points to the 
        %left of right boundary are padded to starting grid. Since
        %sometimes right side also becomes somewhat corrupt due to CDF
        %miscalculations on the right, I add thirty values to the right
        %using a value from ten points deep inside the left boundary. 20
        %values are to be padded to the right and ten extra points are
        %extrapolated in case they have gone bad and a total of thirty
        %points are being extrapolated.
        %wnStart=1;
           
            
            
     
            %plot(ZZ1(1:NN1),BB1(1:NN1),'b',ZZ1(1:NN1),BB1a(1:NN1),'y',ZZ1(1:NN1),BB1b(1:NN1),'k',Z(1:Nn),B1(1:Nn),'r')
            %str=input('Look at B1 and BB1, original and extrapolated version and boundaries');
            
            %ww1(wnStart1:NN1)=BB1(wnStart1:NN1);
 
            %[Muw,Sigmaw] = CalculateDriftAndVolA404(ww1,wnStart1,NN1,YqCoeff0,Fp1,gamma,dt);
            
            [Muw,Sigmaw,Sigma2w,Sigma3w] = CalculateDriftAndVolA404B4Drift02(BB1,wnStart1,NN1,YqCoeff0,Fp1,gamma,dt);

    %        plot(ZZ1(wnStart1:NN1),ww1(wnStart1:NN1),'b',Z(wnStart1:Nn),w(wnStart1:Nn),'r')
            tt
    %        str=input('Look at comparison of interpolated and original w-0');

             %plot(ZZ1(1:NN1),Sigmaw(1:NN1),'b',Z(1:Nn),c1(1:Nn),'r')
    %         plot(ZZ1(1:NN1),Sigmaw(1:NN1),'b',ZZ1(1:NN1),dSigmawdw(1:NN1),'r')
    %        str=input('Look at comparison of original and interpolated vol coefficient-1');
            %plot(ZZ1(1:NN1),cc1(1:NN1),'b',Z(1:Nn),c1(1:Nn),'r')
            %str=input('Look at comparison of original and interpolated vol coefficient-1A');
            %plot(ZZ1(1:NN1),Muw(1:NN1),'b',Z(1:Nn),wMu0dt(1:Nn),'r')
     %       plot(ZZ1(1:NN1),Muw(1:NN1),'b',ZZ1(1:NN1),dMuwdw(1:NN1),'r')
     %       str=input('Look at drifts origianl and interpolated-2');
            %plot(ZZ1(1:NN1),ZZProb(1:NN1),'b',Z(1:Nn),ZProb(1:Nn),'r')
            %str=input('Look at Zprobss -3');
            %plot((1:Nn),Z(1:Nn),'r',(1:NN1)-20,ZZ1(1:NN1),'b')
            %str=input('Look at Zs -4');
       
            
            c1Guess(wnStart:Nn)=sigma0*sqrt(dt);
            Sigmaw_0(wnStart:Nn)=sqrt(dwdZ_0(wnStart:Nn).^2+ c1Guess(wnStart:Nn).^2);
            
           
%            plot(Z(wnStart:Nn),dwdZ_0(wnStart:Nn),'r');
%            str=input('Look at dwdZ_0 before addition to Bn');
            
%            plot(Z(wnStart:Nn),wMu0dt(wnStart:Nn),'r');
%            str=input('Look at dwdZ_0 before addition to Bn');
            

%Below: we project our (guess) target grid using estimate of volatility and
%drift which are being added to target grid at previous point.
            
            Bn(wnStart:Nn)=w(wnStart:Nn)+wMu0dt(wnStart:Nn)+(Sigmaw_0(wnStart:Nn)-dwdZ_0(wnStart:Nn)).*Z(wnStart:Nn);
            Bn(Bn<=0)=0.0;%000000000001;
            
      %      plot((wnStart:Nn),Bn(wnStart:Nn),'r',(wnStart:Nn),w(wnStart:Nn),'g')
      %      str=input('Look at Bn and w');
            
            
            
            wnStart
            
            %Below: After calculating target grid, we recalculate wnStart
            %based on the criteria that as we move to right on the grid,
            %value of grid has to increase and if decreases somewhere, our
            %grid has become corrupt there and we make that point as
            %starting point of caclulations and assign it wnStart.
            
            for nn=NnMidl:-1:wnStart+1
                if(Bn(nn)<Bn(nn-1))
                    %Bn(nn-1)=Bn(nn);
                    Bn(1:nn-1)=Bn(nn);
                    wnStart=nn;
                end
            end
%            plot((1:Nn),Bn(1:Nn),'r');
%            str=input('Look at graphs Bn-2--After');
%            plot((1:Nn),Z(1:Nn),'r');
%            str=input('Look at graphs Z--1');
            
            
            
            
            
%            plot(Z(wnStart:Nn),Bn(wnStart:Nn),'r');
%            str=input('Look at Bn');
            
%Below we have the starting grid that is used as starting point for 
%probability calculations.            
            Bm0(wnStart1:NN1)=BB1(wnStart1:NN1);%
   
            Bm0(Bm0<0)=0.0;%00000000001;

            %Below: Again on starting grid, we recalculate wnStart1
            %based on the criteria that as we move to right on the grid,
            %value of grid has to increase and if decreases somewhere, our
            %grid has become corrupt there and we make that point as
            %starting point of caclulations and assign it wnStart1.
            
            for mm=Nn1Midl:-1:wnStart1+1
                if(Bm0(mm)<Bm0(mm-1))
                    Bm0(1:mm-1)=Bm0(mm);
                    wnStart1=mm;
                end
            end
            
            
%            plot(ZZ1(wnStart1:NN1),Bm0(wnStart1:NN1),'r');
%            str=input('Look at Bm');
            
       %     plot((1:Nn),c1(1:Nn),'r');
       %     str=input('Look at Bn plot---0');
       %     plot((1:Nn),SigmaB(1:Nn),'r',(1:Nn),dBdZ_0(1:Nn),'b');
       %     str=input('Look at Bn plot---1');
       %     plot((1:Nn),Bn(1:Nn));
       %     str=input('Look at Bn plot');
            
           Pn(1:Nn)=0.0;
            
            ZZProb(1:NN1)=Z1Prob(1:Nn1);
            
       %Below we assign all the probability mass to the left of the left 
       %boundary of starting grid to the first cell with index wnStart1 that 
       %is used to make the calculations. This way we make sure that lost
       %probability(that is not on the calculation grid) is accounted for 
       %in the calculations of CDF
            if(wnStart1>1)
                ZZProb(wnStart1)=Z1Prob(wnStart1)+sum(Z1Prob(1:wnStart1-1));
                ZZProb(1:wnStart1-1)=0.0;
            end
            
            

            for mm=wnStart1:NN1

                %Below threshold is calculated using number of points on
                %the padded starting grid that are extra on top of wnstart
                %of the target grid of the previous time. These points were
                %padded using linear extrapolation and now we want to use them 
                %with zero drift and constant first order volatility.
                Threshold=wnStart-wnStart1+20-1;
              
                %Below drift is turned to zero in the linearly extrapolated
                %padded zone and volatility is turned to first order value.
                 if(mm<=wnStart1+Threshold)
                         Muw(mm)=0.0;
                         Sigmaw(mm)=sigma0*sqrt(dt);
                 end

                
                nn1=wnStart;
                for nn=wnStart:Nn
                    if( (Bn(nn)-Bm0(mm)-Muw(mm)) <-5*abs(Sigmaw(mm)))
                        nn1=nn;
                    end
                end

                nn2=Nn;
                for nn=Nn:-1:wnStart
                    if( (Bn(nn)-Bm0(mm)-Muw(mm)) >5*abs(Sigmaw(mm)))
                        nn2=nn;
                    end
                end
                Zt(1:Nn)=0.0;
                Pmn(1:Nn)=0.0;

                if(mm<=wnStart1+Threshold+1)
                    Zt(nn1:nn2)=(Bn(nn1:nn2)-Bm0(mm)-Muw(mm))/Sigmaw(mm);
                else
                    %Below is quadratic formula related to Muller's method
                    %I beleive. It used to be on wikipedia but is not
                    %there now and somebody has probably removed it from 
                    %there.
                    if((Bn(nn1:nn2)-Bm0(mm)-Muw(mm)+Sigma2w(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(nn1:nn2)=(2*(-Bn(nn1:nn2)+Bm0(mm)+Muw(mm)-Sigma2w(mm)))./(-Sigmaw(mm) + sqrt(Sigmaw(mm).^2+4*Sigma2w(mm).*(-Bn(nn1:nn2)+Bm0(mm)+Muw(mm)-Sigma2w(mm))));
                    else
                        Zt(nn1:nn2)=(2*(-Bn(nn1:nn2)+Bm0(mm)+Muw(mm)-Sigma2w(mm)))./(-Sigmaw(mm) - sqrt(Sigmaw(mm).^2-4*Sigma2w(mm).*(-Bn(nn1:nn2)+Bm0(mm)+Muw(mm)-Sigma2w(mm))));
                    end

                    
                end
                Zt(nn1:nn2)=real(Zt(nn1:nn2));
                Pm0n(nn1:nn2)=normcdf(Zt(nn1:nn2),0,1);
                
                Pmn(nn1:nn2)=Pm0n(nn1:nn2).*ZZProb(mm);
                
                
                if(nn1>wnStart)
                    Pmn(1:nn1-1)=0;
                end
                if (nn2<Nn)
                    Pmn(nn2+1:Nn)=ZZProb(mm);
                
                end
                
                
                %Pn is total CDF at the point Bn while Pmn is contribution
                %to CDF from mth grid cell at previous time point.
                Pmn(isnan(Pmn)==1)=0.0;
                
                Pn(wnStart:Nn)=Pn(wnStart:Nn)+Pmn(wnStart:Nn);
                
                
            end
            
%Calculate Znew which is Znew-grid that corresponds to our Guess Bn where CDF is now available.                
            Znew(wnStart:Nn)=norminv(Pn(wnStart:Nn));
            wnStart
            for nn=1:NnMid
                if(Bn(nn)<0)
                    wnStart=nn+1;
                end
            end
       [w] = InterpolateNewGrid(Bn,Znew,wnStart,Nn,dNn,Z);
       
       
       wnStartA=wnStart
      % Nn
      % 1
       
       for nn=1:NnMid
                if(w(nn)<0)
                    wnStart=nn;
                end
            end
            
       wnStartB=wnStart
      % Nn
      % 2
      % if(wnStartA~=wnStartB)
      %     str=input('Different values of wnStart detected');
      % end
       
%       plot((wnStart:Nn),Znew(wnStart:Nn),'g',(wnStart:Nn),Z(wnStart:Nn),'r');
%       str=input('Look at old and interpolated Zs');
     
       w(w<0)=0.0;
       
 %       plot((1:Nn),Bn(1:Nn),'g',(1:Nn),w(1:Nn),'r');
 %           str=input('Look at w-Bn graphs after iterpolation');
       
       
       
       [dwdZ_1,d2wdZ2_1,d3wdZ3_1] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z);
       
%       plot(Z(1:Nn),dwdZ_1(1:Nn),'r');
%       str=input('Look at dwdZ graph');             
       %plot((1:Nn),Bn(1:Nn),'g',(1:Nn),Bn2(1:Nn),'r');
       %     str=input('Look at Bn graphs');

       
       
       %Below: We linearly re-calculate by linear extrapolation 
       %the last 20/30 cells on the grid 
       %next to the right boundary since sometimes rarely due to CDF
       %miscalculations, and when these points are recalculated, the grid 
       %remains stable. This is usually not needed but when the process
       %comes under stress due to miscalculations on the left side, this
       %helps to keep everything in sync.
       
       
       if(wnStart1<=35)
%       w(Nn-20+1:Nn)=w(Nn-20)+(1:20).*dwdZ_1(Nn-20)*dNn;
       end
       
  %     if(wnStart1>35)
%       w(Nn-30+1:Nn)=w(Nn-30)+(1:30).*dwdZ_1(Nn-30)*dNn;
  %     end
       
       
 
%Below is the grsph of newly calculated value of w.       
%       plot((1:Nn),w(1:Nn),'r')
%       str=input('Calculated Value of w');
    end
end

%Below: After the final point in simulation has reached, 
%Since a lot of points on the left side of the density have been
%cut, we reconstruct them with linear extrapolation. This extrapolation is
%relatively poor, replace with your favorite method to extrapolate the
%first few points of the density.
[dwdZ_2,d2wdZ2_2,d3wdZ3_2] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z);
  %w(1:wnStart-1)=w(wnStart)-(wnStart-1:-1:1).*dwdZ_2(wnStart)*dNn;
  w(1:wnStart+1)=w(wnStart+2)-(wnStart+1:-1:1).*dwdZ_2(wnStart+2)*dNn;%
  wnStart=1;
     for nn=1:NnMid
                if(w(nn)<=0)
                    wnStart=nn;
                end
     end
  



%below D's (the names of variables starting with D) are 
%change of probability derivatives.

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

rng(29079137, 'twister')
paths=200000;
YY(1:paths)=x0;  %Original process monte carlo.
Random1(1:paths)=0;
%for tt=1:(Tt-1)+8
for tt=1:TtM    
if(tt>Tt1)
    %dtM=dt2/8;
end
    Random1=randn(size(Random1));
    HermiteP1(1,1:paths)=1;
    HermiteP1(2,1:paths)=Random1(1:paths);
    HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
    HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
    HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;
    

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

end

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

MaxCutOff=30;
NoOfBins=round(300*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too 
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );

%plot(y_w1(wnStart+1:Nn-1),py_w1(wnStart+1:Nn-1),'b',y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
plot(y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');

title(sprintf('x0 = %.4f,theta=%.4f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.3f,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));
 
legend({'New Model Density','Monte Carlo Density'},'Location','northeast')
 
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
end
.
Here is another function you will need.
.
.
function [wMu0dt,c1,c2,c3] = CalculateDriftAndVolA404B4Drift02(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,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;

   
   

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

 c2(wnStart:Nn)=((YqCoeff0(1,1,1,3).*yy(wnStart:Nn).^Fp1(1,1,1,3) *dt) + ...
     (YqCoeff0(1,1,2,3).*yy(wnStart:Nn).^Fp1(1,1,2,3)+ ...
     YqCoeff0(1,2,1,3).*yy(wnStart:Nn).^Fp1(1,2,1,3)+ ...
     YqCoeff0(2,1,1,3).*yy(wnStart:Nn).^Fp1(2,1,1,3)).*dt^2+ ...
     (YqCoeff0(1,1,3,3).*yy(wnStart:Nn).^Fp1(1,1,3,3)+ ...
     YqCoeff0(1,2,2,3).*yy(wnStart:Nn).^Fp1(1,2,2,3)+ ...
     YqCoeff0(2,1,2,3).*yy(wnStart:Nn).^Fp1(2,1,2,3) + ...
     YqCoeff0(1,3,1,3).*yy(wnStart:Nn).^Fp1(1,3,1,3)+ ...
     YqCoeff0(2,2,1,3).*yy(wnStart:Nn).^Fp1(2,2,1,3)+ ...
     YqCoeff0(3,1,1,3).*yy(wnStart:Nn).^Fp1(3,1,1,3)).*dt^3);
 
 c3(wnStart:Nn)=((YqCoeff0(1,1,1,4).*yy(wnStart:Nn).^Fp1(1,1,1,4)*dt^1.5 )+ ...
     (YqCoeff0(1,1,2,4).*yy(wnStart:Nn).^Fp1(1,1,2,4)+ ...
     YqCoeff0(1,2,1,4).*yy(wnStart:Nn).^Fp1(1,2,1,4)+ ...
     YqCoeff0(2,1,1,4).*yy(wnStart:Nn).^Fp1(2,1,1,4))*dt^2.5);
 

end
.
.
Here is the final output of the program when you will run it.

ItoHermiteMean =
   0.014086993873228
true Mean only applicable to standard SV mean reverting type models otherwise disregard
TrueMean =
   0.014000000000000
Original process average from monte carlo
MCMean =
   0.014070563722125
MCVar =
     6.298015695360379e-04
Original process average from our simulation
ItoHermiteMean =
   0.014086993873228
ItoHermiteVar =
     6.098922224266122e-04
true Mean only applicble to standard SV mean reverting type models otherwise disregard
TrueMean =
   0.014000000000000
IndexMax =
        3656
red line is density of SDE from Ito-Hermite method, green is monte carlo.

and here is the graph you will get.

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

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

April 23rd, 2021, 2:55 pm

Many times when SDE is in a very close range to zero and you want to use very high volatility and you notice that simulation has become unstable, it helps to decrease the time step size for simulation. You can sometimes drastically increase the vol and still be able to simulate if you drastically decrease the step size. For example I could never simulate the two high volatility scenarios below but when I decreased step size from 1/32 to 1/100, I was able to simulate them perfectly. Here are the graphs with parameters for two high volatility scenarios I am mentioning.

Image

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

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

April 24th, 2021, 3:55 am

With this method in previous posts, we can also tackle probably SDEs that have more than one diffusion term but it would require some hack. It is probably not simple to convert such SDEs to Bessel coordinates with Lamperti transform but we do not necessarily need everything in Bessel coordinates. We can still convert one of the two diffusion terms to a constant form considering which one of the two diffusion terms is easier to deal with.
If the target point at time t+1 is Bn and starting point at time t is Bm. we can have the relationship

[$]B_n=B_m + \mu(B_m) + \sigma_{11}(B_m) Z_1 + \sigma_{12}(B_m) Z_2[$]
The above relationship can be tackled since coefficients of Z's become constant at any particular point on starting density grid and can be added in squared fashion and replaced by a single Z.
But that will work well only for low volatilities and we will have to add a second hermite polynomial to above expression when volatilities are even slightly above average and our equation would become
[$]B_n=B_m + \mu(B_m) + \sigma_{11}(B_m) Z_1 + \sigma_{12}(B_m) Z_2+ \sigma_{21}(B_m) ({Z_1}^2-1) + \sigma_{22}(B_m) ({Z_2}^2-1)[$]
but here sum of two squared normals with different variances would be a generalized chi-square density and some hack would have to be used to tackle it but I am sure it can somehow be done. I thought it would be interesting to share it with friends. I did not think hard but may be we can replace the sum of square gaussian terms with some function of square of resulting gaussian( that is sum of two gaussians). 
 
User avatar
Amin
Topic Author
Posts: 2782
Joined: July 14th, 2002, 3:00 am

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

April 28th, 2021, 10:35 am

Dear friends, I am working on solution of stochastic path integrals of SDEs (that I have been previously calling dt-integrals and dz_integrals) as an extension of the method for the calculation of the solution of densities of SDEs that we recently did in previous few posts. I have done a lot of the work and I hope to post the programs for calculation of above stochastic integrals in a a few(2-3) days. 
 
User avatar
Amin
Topic Author
Posts: 2782
Joined: July 14th, 2002, 3:00 am

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

April 28th, 2021, 1:12 pm

I made this post on my off-topic thread and I am copying it here and want to request friends to help me weather this summer by forcing mind control agencies to be better with me. Please ask American embassy in Pakistan to not vigorously lobby for my persecution. I am copying my post here from off-topic forum where I originally posted it:https://forum.wilmott.com/viewtopic.php?f=15&t=94796&p=865546#p865546

Friends, I have not written on this thread for a long while now but it really does not mean that my persecution has stopped. Persecution activity continues but at the same time it is not extreme either. Water is good at a lot of places but at the same time water is increasingly drugged in some parts of the neighborhoods where I live. I found that five liter nestle water bottle had good water and started using but they have started to catch up and have drugged it in some parts of  neighborhoods close to where I live. For example today when I bought 5 liter water bottle, I believed it to be good as usual only to find out after a few hours that my consciousness was decreasing and I was being extremely sluggish and I realized that it was due to drugged water. Water or any drink that I like gets drugged in all close neighborhoods in a few(1-2) weeks. Though this is not good but this is not extremely drastic as compared to bad times when I had to drive around the city to get good water and I would be able to find good water in remote neighborhoods only when I would be lucky somewhere. Now when I try to get good water I can find it somewhere in 1-2 hours but I am afraid that it would continue to become increasingly difficult since mind control agencies continue to study my water finding patters and continue to slowly work to drug places where I can get good water. Summer is already here and trying to get good water and drinks in scorching heat is going to be very difficult especially when I do not use air conditioning in my car.
Another thing I wanted to tell friends that when President Biden won, I was so sure that I would be using air conditioning in coming summer and would do my morning walk in T-shirt and shorts like normal. I use air-conditioning in my room extremely rarely due to huge amount of sickening gas that comes out of air-conditioner in my room and in my car. In fact for at least five-six years in the past fifteen years, I slept on a cot in the open in the extended balcony of our house to be able to weather the scorching heat of the night when I could not use air-conditioning(And they would many times respond by releasing gases out of the pedestal fan that I would use to cool me in the open balcony). And for past six years, I rarely wear T-shirts and shorts in summers because exposed parts of our body can easily help us get into strong mind control and decreasing consciousness and I would have to wear cotton pants and shirts even though it would be very difficult to go out in scorching heat with all the sweat in cotton clothes. I would love to go back to times when I could wear simple shorts and T-shirts (or other loose clothes at night) and enjoy my morning walk even in summers but this cannot be a possibility until there is a marked decrease in mind control. 
I also want to take this opportunity to request American friends to force mind control crooks to not use gases in my air-conditioning this summer. Though I am sure that even if they decrease use of gases in AC due to pressure from friends, they would still always use some gas of lower intensity that would slow me down but would be harder to detect right away and would be less sickening to be able to be noticed right away.
 
User avatar
Amin
Topic Author
Posts: 2782
Joined: July 14th, 2002, 3:00 am

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

May 3rd, 2021, 9:51 am

Friends, I have been able to calculate the arithmetic sum path integrals(dt-integrals type) for mean reverting diffusions quite well. I changed some parts of my strategy to solve the integrals and that seems to work remarkably well. It especially seems to work very well close to zero even for kappa=2 and higher but performance does deteriorate slightly for high kappas(mean reversion speed) as we get to values close to one and higher. I was just recently able to make the breakthrough as I changed my strategy otherwise, employing methods that I have presented before, I was getting rather poor results. But now results are strikingly precise for a large range of parameters. For cases that do not work with precision, I think I know the reason and I will be fixing them in next 2-3 days and improving the method there. I hope to post the new programs in 3-4 days. 
 
User avatar
Amin
Topic Author
Posts: 2782
Joined: July 14th, 2002, 3:00 am

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

May 7th, 2021, 5:24 pm

Friends, I think I have found a way to write ODE for constant CDF lines of stochastic differential equations. As a start I applied it to brownian motion from origin and was able to retrieve the ODE whose solution is constant CDF lines of brownian motion. I was able to do it with the density of brownian motion. Before I write the details of the method, it is important to recognize that when we write the density of brownian motion, we need to emphasize that brownian motion is a function of t. I took CDF of the brownian motion density by integrating it from -infinity to Z(t) and then differentiated the resulting equation with respect to t and then equated all of that to zero (meaning CDF does not change as we take its derivative with respect to time.) and when I simplified the equation to get rid of the common density and other terms, it was a simple ODE And solving the ODE led to constant CDF lines of brownian motion. Now I will explain that with equations and also walk the friends through mathematica steps I used to find the results. Density of brownian motion emphasizing that Z is dependent on t would be. 

[$]\frac{1}{\sqrt{2 \pi t}} \exp[-.5 \frac{{Z(t)}^2}{t}] [$]


I found the CDF of the density as

[$]\int_{-\infty}^{Z(t)} \frac{1}{\sqrt{2 \pi t}} \exp[-.5 \frac{{Z(t)}^2}{t}] dZ[$]

I differentiated the above equation in mathematica where I had emphasized the dependence of Z[t] on t and equated it to zero as previously suggested. I am writing the equation as

[$]\frac{d}{dt}[\int_{-\infty}^{Z(t)} \frac{1}{\sqrt{2 \pi t}} \exp[-.5 \frac{{Z(t)}^2}{t}] dZ]=0[$]

After a few manipulations, I got the simple ODE as a result of the above equation as

[$]-0.5 Z(t) + t \frac{dZ(t)}{dt} =0[$]

I solved the above simple ODE in mathematica to get the answer which is also solution for constant CDF lines of brownian motion as

[$]{{Z[t]->t^{0.5} C[1]}}[$]

I am sure it would easily be possible to write the normal form of evolution SDEs in general and then solve them to find an ODE that can be cheaply solved with this using numerical methods(to give constant CDF line solutions) as I am afraid that getting analytical results for all of the SDEs evolution on constant CDF lines would be impossible.

I used mathematica for all of this and I am going to walk the friends through mathematica commands used to find the above solution.

First use the command

1.  
D[Integrate[1/Sqrt[2 Pi t] Exp[-.5  (Z[t]^2)/t], {Z[t], -Infinity, Z[t]}], t]

It will give you a relatively complicated solution. I am not writing the solution but I represent it here with 'Solution' in mathematica code.

2.
Simplify[Solution == 0]

The above equation means that we have asked mathematica to equate the result of command 1 with zero and simplify it. (You have to copy the result of command 1 and copy it in command 2 in place of 'solution' words I have written there)

3.  As a result of simplification, you will get the following (I have written it in Latex form now.

[$]\frac{e^{-\frac{0.5 Z(t)^2}{t}} \left( t  Z'(t)-0.5 Z(t)\right)}{\sqrt{t}}=0[$]

from which you easily retrieve the ODE

[$]t Z'(t)-0.5 Z(t)=0[$]

4. I used the mathematica command to solve this ODE as

DSolve[-0.5 Z[t]+t Z'[t] == 0, Z[t], t]

which gave me the solution to constant CDF lines of brownian motion as

[$]{{Z[t] -> t^{0.5} C[1]}}[$]


In next few days, I will spend time trying to write ODEs for various SDEs and solving them numerically and would be pasting the code here. I have placed the stochastic integrals project for a small hold but many things we will learn here would be useful in solution of stochastic integrals as well as thinking about them led me to this. I will be posting a lot of code in next few days.

I have not done this so I could not be perfectly certain but I believe that you do not need to know the universal density of SDE diffusions as they are mostly not known but I think that we could just have a density representation of instantaneous evolution of SDEs over small time steps and then use that with numerical solution of ODEs to find the solution of SDEs along constant CDF lines. I will keep friends posted as I do more work on this.
 
User avatar
Amin
Topic Author
Posts: 2782
Joined: July 14th, 2002, 3:00 am

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

May 7th, 2021, 6:24 pm

I want to take this opportunity to thank friends who tried to protect me from severe persecution. If it were not for their efforts, I would never have been able to do any decent research in mathematics. Really a lot of time the persecution activity markedly decreased due to pressure from friends and I was more able to concentrate on my research. Though my persecution still continues and sometime there is marked increase and especially sometimes it becomes difficult to get good water but this is very little as compared to bad times I have seen previously. And my persecution would have been very extreme  if it were not for pressure from good people and their efforts to decrease my persecution. I sincerely thank and I am seriously indebted to all friends and good people who tried to help end my persecution or decrease it. I do hope that persecution and victimization of a large number of totally innocent intelligent people with mind control finally ends and all these good and innocent victims are allowed to freely make their contribution to our society in a good way. I want to thank all good people who stood against mind control  and would request them to continue their efforts until this ill ends from our societies. There are a large number of mind control victims who are targeted only because some bad people simply do not want the intelligent victims to find expression for themselves in a good and positive way which is something extremely cruel of these bad people.
 
User avatar
Amin
Topic Author
Posts: 2782
Joined: July 14th, 2002, 3:00 am

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

May 8th, 2021, 9:06 am

What I have been calling Ito-Taylor expansion (and I also used the same method for monte carlo simulation) is basically a Taylor expansion with repeated application of Ito change of variable formula. Many friends would have read my paper on ODEs in which I used the very same Taylor expansion method(Iterated Integrals) for ODEs. We all know that ODEs work very well with Taylor expansion method and we can take a large number of time steps and apply the series expansion at each step to finally calculate the solution of ODE at terminal time. Barring monte carlo which is different, my attempts at multistep Taylor expansions of SDEs never worked well despite several years of attempts. One step analytic Taylor expansion of SDEs could work well for several years but when I would try the same expansion with smaller multiple time steps(on the normal Z-grid) accuracy would always decrease in my experiments. It was only when I approximately solved the Fokker-planck equation using derivatives of the PDE that I was able to get a tree that branches once at origin and continues to advance each branch in time on the Z-grid. I always wondered why we could not use analytic Taylor expansion with small multiple steps successfully on a Z-grid with small time steps. But I think that I know the reason behind it now. When Taylor coefficients are evaluated at starting time zero, they work well but when they are evaluated in forward time, they do not work perfectly well unlike in ODEs where they work very well. I have some ideas again and I want to spend some time on them again using stochastic integral formulas that adapt to changing variance(due to drift) and also looking at the fact that forward coefficients are themselves stochastic and there evolution as coefficients of hermite polynomials have to be seen as stochastic evolution of product of two random variables. If I can get it to work, it may be the simplest way to solve the SDEs and their stochastic integrals. 
 
User avatar
Amin
Topic Author
Posts: 2782
Joined: July 14th, 2002, 3:00 am

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

May 9th, 2021, 7:58 am

Friends, I was doing research on ODE like formulas for constant CDF lines of the SDE densities. I realized that only transition SDE is probably not enough to calculate the constant CDF lines and I wanted to explore ways to find a complete density for small steps that takes into account the previous density. I was able to find a simple series formula that could update the entire density for short steps using convolution. It would not work for long steps since I am only using first hermite though I really think that it could be possible to extend this series formula to two hermite polynomials in volatility with a bit of hack( I will share that if that works in my experiments since I still cannot be sure.). Here is the derivation of the new formula.

 We are in Bessel/Lamperti format and we are taking a constant  volatility [$]\sigma[$]. I write the SDE as

[$]dy(t)= \mu(y) dt + \sigma dz(t)[$]

In general when we expand the SDE using Taylor, we get multiple terms in drift that can have several powers of t and similarly for volatility so I write the above equation in discrete form and use symbols  [$]\mu[$] and [$]\sigma[$]  to represent sum of all the terms and I have absorbed t or dt in the same symbol representation since it will help keep focus on the main derivation so I write the discrete form of SDE as

[$]y_2=y_1+ \mu(y_1,t)  + \sigma(t) Z [$]
Here [$]y_2[$] denotes the SDE variable at time [$]t_2[$] and [$]y_1[$] is the starting density variable of SDE at time [$]t_1[$].

We can write the evolution equation for the density at next step using convolution as

[$]p(y_2) =\int_{-\infty}^{\infty} f(y_1) g(\frac{(y_2-y_1- \mu(y_1,t))}{ \sigma(t) }) dy_1[$]

here [$]f(y_1)[$] is the original density at time [$]t_1[$] and [$]g(.)[$] is a normal density.
We write the above integral equation with normal density explicitly as

[$]p(y_2) =\int_{-\infty}^{\infty} f(y_1)  \frac{1}{\sigma \sqrt{2 \pi}} \exp(-.5 \big [\frac{(y_2-y_1- \mu(y_1,t))}{ \sigma(t) } \big ]^2 ) dy_1[$]

Expanding the terms inside the normal exponential, we get

[$]p(y_2) =\int_{-\infty}^{\infty} f(y_1)  \frac{1}{\sigma \sqrt{2 \pi}} \exp(-.5 \big [\frac{( {y_2}^2- 2 y_2(y_1+\mu(y_1,t))+ {(\mu(y_1,t)+{y_1})}^2 ) }{ {\sigma(t)}^2 } \big ] ) dy_1[$]

and further re-arrangement of terms gives

[$]p(y_2) =\frac{1}{\sigma \sqrt{2 \pi}}   \exp(-.5 \big [\frac{ {y_2}^2}{ {\sigma(t)}^2 } \big ]) \int_{-\infty}^{\infty} f(y_1)  \exp(-.5 \big [\frac{(- 2 y_2(y_1+\mu(y_1,t)) }{ {\sigma(t)}^2 } \big ]) \exp(-.5 \big [\frac{( {(\mu(y_1,t)+{y_1})}^2 ) }{ {\sigma(t)}^2 } \big ] ) dy_1[$]

We can expand the integral as expectation of a series evaluated with respect to density  [$]f(y_1)[$]. Below whenever I write expectation symbol E, it means that expectation is taken over the density [$]f(y_1)[$], so we can write the above integral as an expansion as

[$]p(y_2) =\frac{1}{\sigma(t) \sqrt{2 \pi}}   \exp(-.5 \big [\frac{ {y_2}^2}{ {\sigma(t)}^2 } \big ]) E \Big [\exp(-.5 \big [\frac{( {(\mu(y_1,t)+{y_1})}^2 ) }{ {\sigma(t)}^2 } \big ] )\Big ] *(1 + y_2 E \big [(\frac{(- 2 (y_1+\mu(y_1,t)) }{ {\sigma(t)}^2 }) \big ] [$]
[$]+ \frac{ {y_2}^2}{2}  E \big [(\frac{(- 2 (y_1+\mu(y_1,t)) }{ {\sigma(t)}^2 }) \big ]^2 + \frac{ {y_2}^3}{6}  E \big [(\frac{(- 2 (y_1+\mu(y_1,t)) }{ {\sigma(t)}^2 }) \big ]^3 + . . .) [$]

please notice that I have taken [$]\frac{1}{\sigma(t) \sqrt{2 \pi}}[$]  out of integral since even if we have expanded [$]\sigma(t)[$] to higher order terms in t, its dependence on [$]y_1[$] is totally negligible in Bessel coordinates unless we are extremely close to zero.
So this way we can easily expand the resulting density of SDE to next steps but we will have to keep time steps reasonable. I have not shown this to keep exposition simple, but the powers of time terms would also explicitly appear in the expansion and can be easily differentiated to work with the ODEs etc. We may have to take a lot of terms in the expansion depending upon what the parameters are but it should not be difficult at all.
So eventually we get an expansion for the density of  [$]y_2[$]  in the form

[$]p(y_2) =C_0  \exp(-.5 \big [\frac{ {y_2}^2}{ {\sigma(t)}^2 } \big ])  *(1 +C_1 y_2  +C_2 \frac{ {y_2}^2}{2}  + C_3 \frac{ {y_2}^3}{6} + . . .) [$]
 where all the constants [$]C_0[$],[$]C_1[$],[$]C_2[$] and further are evaluated as expectations over the original density of [$]y_1[$] at starting time.