SERVING THE QUANTITATIVE FINANCE COMMUNITY

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

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

July 27th, 2020, 4:50 pm

Please note that I posted a program in last post on previous page just now. This post refers to interpolation used in that program.
With straightforward Lagrange interpolation, friends can easily choose any grid they wish. Also higher kappa values which become unstable seem to become somewhat better when we decrease the time step size. For example, I was able to get a program run with kappa=4.0 and gamma=.75 when I halved the time step. I was not being able to get a stable run for these parameter values with the original time stepping as in program I uploaded in previous post. Here is the graph of the program I mentioned. Parameters are on graph.

Image

Here is another program run with large mean-reversion that became possible when I halved the time step.

Image
Last edited by Amin on July 27th, 2020, 5:12 pm, edited 1 time in total.
 
User avatar
Amin
Topic Author
Posts: 2536
Joined: July 14th, 2002, 3:00 am

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

July 27th, 2020, 5:09 pm

Just another thing I noticed. It seems we have to double C33 constant in program when we halve the time step size. It (constant C33) might possibly have nothing to do with dt though it is approximately at right value for given step size in the program posted. But anyway please see  if making it double helps when you halve the time step size.
 
User avatar
Amin
Topic Author
Posts: 2536
Joined: July 14th, 2002, 3:00 am

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

July 27th, 2020, 6:33 pm

This is a new modified program. I have not tested it extensively but it looks good and therefore I decided to post it. It works with slightly shorter stepping.
function [TrueMean,ItoHermiteMean,y_w,fy_w,wnStart] = FPERevisited044Wilmott16()

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

%I have not directly simulated the SDE but simulated the transformed 
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.
%In this version, the mean correction has been disabled. If you are
%simulating standard mean-reverting SV type SDEs, please enable the mean 
%correction by uncommenting the appropriate line in the body of code below.

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


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

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

x0=1;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.75;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=4.0;%.950;   %mean reversion parameter.
theta=.250;%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;
                
w(1:Nn)=x0^(1-gamma)/(1-gamma);

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

Z
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);

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

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

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

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

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

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

for k = 0 : (OrderA)
    for m = 0:k
        l4 = k - m + 1;
        for n = 0 : m
            l3 = m - n + 1;
            for j = 0:n
                l2 = n - j + 1;
                l1 = j + 1;
                %Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
                Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                
                YCoeff0(l1,l2,l3,l4) =YCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
                YqCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
            end
        end
    end
end

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

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

for tt=1:Tt
    
     yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
    
    [wMu0dt,dwMu0dtdw,d2wMu0dtdw2,d3wMu0dtdw3,c1,c2,c3,c4] = CalculateDriftAndVolOrder4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);

    wPrev=w;
    dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn) ;% ...
    dw2(wnStart:Nn)=dw(wnStart:Nn).^2;
    w(isnan(w)==1)=0;
    wMu0dt(isnan(wMu0dt)==1)=0;  
    wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
   
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



d2wdZ2(wnStart:Nn)=d2wdZ2(wnStart:Nn)+ d2wMu0dtdw2(wnStart:Nn).*dwdZ(wnStart:Nn).^2+dwMu0dtdw(wnStart:Nn).*d2wdZ2(wnStart:Nn)+.33*(Z(wnStart:Nn)).*sigma0^2*dt^2;;%+ ...


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


if (tt>16)
d2wdZ2=d2wdZ2A;

      dwdZ(NnMidh+3)=InterpolateOrderN8(8,Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),dwdZ(NnMidh+4),dwdZ(NnMidh+5),dwdZ(NnMidh+6),dwdZ(NnMidh+7),dwdZ(NnMidl-3),dwdZ(NnMidl-4),dwdZ(NnMidl-5),dwdZ(NnMidl-6));  
      dwdZ(NnMidh+2)=InterpolateOrderN8(8,Z(NnMidh+2),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),dwdZ(NnMidh+4),dwdZ(NnMidh+5),dwdZ(NnMidh+6),dwdZ(NnMidh+7),dwdZ(NnMidl-3),dwdZ(NnMidl-4),dwdZ(NnMidl-5),dwdZ(NnMidl-6));  
      dwdZ(NnMidh+1)=InterpolateOrderN8(8,Z(NnMidh+1),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),dwdZ(NnMidh+4),dwdZ(NnMidh+5),dwdZ(NnMidh+6),dwdZ(NnMidh+7),dwdZ(NnMidl-3),dwdZ(NnMidl-4),dwdZ(NnMidl-5),dwdZ(NnMidl-6));  
      dwdZ(NnMidh)=InterpolateOrderN8(8,Z(NnMidh),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),dwdZ(NnMidh+4),dwdZ(NnMidh+5),dwdZ(NnMidh+6),dwdZ(NnMidh+7),dwdZ(NnMidl-3),dwdZ(NnMidl-4),dwdZ(NnMidl-5),dwdZ(NnMidl-6));  
      dwdZ(NnMidl)=InterpolateOrderN8(8,Z(NnMidl),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),dwdZ(NnMidh+4),dwdZ(NnMidh+5),dwdZ(NnMidh+6),dwdZ(NnMidh+7),dwdZ(NnMidl-3),dwdZ(NnMidl-4),dwdZ(NnMidl-5),dwdZ(NnMidl-6));  
      dwdZ(NnMidl-1)=InterpolateOrderN8(8,Z(NnMidl-1),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),dwdZ(NnMidh+4),dwdZ(NnMidh+5),dwdZ(NnMidh+6),dwdZ(NnMidh+7),dwdZ(NnMidl-3),dwdZ(NnMidl-4),dwdZ(NnMidl-5),dwdZ(NnMidl-6));  
      dwdZ(NnMidl-2)=InterpolateOrderN8(8,Z(NnMidl-2),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),dwdZ(NnMidh+4),dwdZ(NnMidh+5),dwdZ(NnMidh+6),dwdZ(NnMidh+7),dwdZ(NnMidl-3),dwdZ(NnMidl-4),dwdZ(NnMidl-5),dwdZ(NnMidl-6));  
%       
%       dwdZ(NnMidh+2)=InterpolateOrderN8(8,Z(NnMidh+2),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),dwdZ(NnMidh+3),dwdZ(NnMidh+4),dwdZ(NnMidh+5),dwdZ(NnMidh+6),dwdZ(NnMidl-2),dwdZ(NnMidl-3),dwdZ(NnMidl-4),dwdZ(NnMidl-5));
%       dwdZ(NnMidh+1)=InterpolateOrderN8(8,Z(NnMidh+1),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),dwdZ(NnMidh+3),dwdZ(NnMidh+4),dwdZ(NnMidh+5),dwdZ(NnMidh+6),dwdZ(NnMidl-2),dwdZ(NnMidl-3),dwdZ(NnMidl-4),dwdZ(NnMidl-5));
%       dwdZ(NnMidh)=InterpolateOrderN8(8,Z(NnMidh),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),dwdZ(NnMidh+3),dwdZ(NnMidh+4),dwdZ(NnMidh+5),dwdZ(NnMidh+6),dwdZ(NnMidl-2),dwdZ(NnMidl-3),dwdZ(NnMidl-4),dwdZ(NnMidl-5));
%       dwdZ(NnMidl)=InterpolateOrderN8(8,Z(NnMidl),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),dwdZ(NnMidh+3),dwdZ(NnMidh+4),dwdZ(NnMidh+5),dwdZ(NnMidh+6),dwdZ(NnMidl-2),dwdZ(NnMidl-3),dwdZ(NnMidl-4),dwdZ(NnMidl-5));
%       dwdZ(NnMidl-1)=InterpolateOrderN8(8,Z(NnMidl-1),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),dwdZ(NnMidh+3),dwdZ(NnMidh+4),dwdZ(NnMidh+5),dwdZ(NnMidh+6),dwdZ(NnMidl-2),dwdZ(NnMidl-3),dwdZ(NnMidl-4),dwdZ(NnMidl-5));


d2wdZ2(NnMidh+4)=InterpolateOrderN8(8,Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidh+8),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),Z(NnMidl-7),d2wdZ2(NnMidh+5),d2wdZ2(NnMidh+6),d2wdZ2(NnMidh+7),d2wdZ2(NnMidh+8),d2wdZ2(NnMidl-4),d2wdZ2(NnMidl-5),d2wdZ2(NnMidl-6),d2wdZ2(NnMidl-7));  
d2wdZ2(NnMidh+3)=InterpolateOrderN8(8,Z(NnMidh+3),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidh+8),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),Z(NnMidl-7),d2wdZ2(NnMidh+5),d2wdZ2(NnMidh+6),d2wdZ2(NnMidh+7),d2wdZ2(NnMidh+8),d2wdZ2(NnMidl-4),d2wdZ2(NnMidl-5),d2wdZ2(NnMidl-6),d2wdZ2(NnMidl-7));  
d2wdZ2(NnMidh+2)=InterpolateOrderN8(8,Z(NnMidh+2),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidh+8),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),Z(NnMidl-7),d2wdZ2(NnMidh+5),d2wdZ2(NnMidh+6),d2wdZ2(NnMidh+7),d2wdZ2(NnMidh+8),d2wdZ2(NnMidl-4),d2wdZ2(NnMidl-5),d2wdZ2(NnMidl-6),d2wdZ2(NnMidl-7));  
d2wdZ2(NnMidh+1)=InterpolateOrderN8(8,Z(NnMidh+1),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidh+8),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),Z(NnMidl-7),d2wdZ2(NnMidh+5),d2wdZ2(NnMidh+6),d2wdZ2(NnMidh+7),d2wdZ2(NnMidh+8),d2wdZ2(NnMidl-4),d2wdZ2(NnMidl-5),d2wdZ2(NnMidl-6),d2wdZ2(NnMidl-7));  
d2wdZ2(NnMidh)=InterpolateOrderN8(8,Z(NnMidh),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidh+8),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),Z(NnMidl-7),d2wdZ2(NnMidh+5),d2wdZ2(NnMidh+6),d2wdZ2(NnMidh+7),d2wdZ2(NnMidh+8),d2wdZ2(NnMidl-4),d2wdZ2(NnMidl-5),d2wdZ2(NnMidl-6),d2wdZ2(NnMidl-7));  
d2wdZ2(NnMidl)=InterpolateOrderN8(8,Z(NnMidl),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidh+8),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),Z(NnMidl-7),d2wdZ2(NnMidh+5),d2wdZ2(NnMidh+6),d2wdZ2(NnMidh+7),d2wdZ2(NnMidh+8),d2wdZ2(NnMidl-4),d2wdZ2(NnMidl-5),d2wdZ2(NnMidl-6),d2wdZ2(NnMidl-7));  
d2wdZ2(NnMidl-1)=InterpolateOrderN8(8,Z(NnMidl-1),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidh+8),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),Z(NnMidl-7),d2wdZ2(NnMidh+5),d2wdZ2(NnMidh+6),d2wdZ2(NnMidh+7),d2wdZ2(NnMidh+8),d2wdZ2(NnMidl-4),d2wdZ2(NnMidl-5),d2wdZ2(NnMidl-6),d2wdZ2(NnMidl-7));  
d2wdZ2(NnMidl-2)=InterpolateOrderN8(8,Z(NnMidl-2),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidh+8),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),Z(NnMidl-7),d2wdZ2(NnMidh+5),d2wdZ2(NnMidh+6),d2wdZ2(NnMidh+7),d2wdZ2(NnMidh+8),d2wdZ2(NnMidl-4),d2wdZ2(NnMidl-5),d2wdZ2(NnMidl-6),d2wdZ2(NnMidl-7));  
d2wdZ2(NnMidl-3)=InterpolateOrderN8(8,Z(NnMidl-3),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7),Z(NnMidh+8),Z(NnMidl-4),Z(NnMidl-5),Z(NnMidl-6),Z(NnMidl-7),d2wdZ2(NnMidh+5),d2wdZ2(NnMidh+6),d2wdZ2(NnMidh+7),d2wdZ2(NnMidh+8),d2wdZ2(NnMidl-4),d2wdZ2(NnMidl-5),d2wdZ2(NnMidl-6),d2wdZ2(NnMidl-7));  

% 
%       d2wdZ2(NnMidh+2)=InterpolateOrderN8(8,Z(NnMidh+2),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),d2wdZ2(NnMidh+3),d2wdZ2(NnMidh+4),d2wdZ2(NnMidh+5),d2wdZ2(NnMidh+6),d2wdZ2(NnMidl-2),d2wdZ2(NnMidl-3),d2wdZ2(NnMidl-4),d2wdZ2(NnMidl-5));
%       d2wdZ2(NnMidh+1)=InterpolateOrderN8(8,Z(NnMidh+1),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),d2wdZ2(NnMidh+3),d2wdZ2(NnMidh+4),d2wdZ2(NnMidh+5),d2wdZ2(NnMidh+6),d2wdZ2(NnMidl-2),d2wdZ2(NnMidl-3),d2wdZ2(NnMidl-4),d2wdZ2(NnMidl-5));
%       d2wdZ2(NnMidh)=InterpolateOrderN8(8,Z(NnMidh),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),d2wdZ2(NnMidh+3),d2wdZ2(NnMidh+4),d2wdZ2(NnMidh+5),d2wdZ2(NnMidh+6),d2wdZ2(NnMidl-2),d2wdZ2(NnMidl-3),d2wdZ2(NnMidl-4),d2wdZ2(NnMidl-5));
%       d2wdZ2(NnMidl)=InterpolateOrderN8(8,Z(NnMidl),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),d2wdZ2(NnMidh+3),d2wdZ2(NnMidh+4),d2wdZ2(NnMidh+5),d2wdZ2(NnMidh+6),d2wdZ2(NnMidl-2),d2wdZ2(NnMidl-3),d2wdZ2(NnMidl-4),d2wdZ2(NnMidl-5));
%       d2wdZ2(NnMidl-1)=InterpolateOrderN8(8,Z(NnMidl-1),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),d2wdZ2(NnMidh+3),d2wdZ2(NnMidh+4),d2wdZ2(NnMidh+5),d2wdZ2(NnMidh+6),d2wdZ2(NnMidl-2),d2wdZ2(NnMidl-3),d2wdZ2(NnMidl-4),d2wdZ2(NnMidl-5));

%d2wdZ2(21:29)=d2wdZ2(20)*(30-(21:29))/10+d2wdZ2(30)*((21:29)-20)/10;

d2wdZ2(1:15)=d2wdZ2(1:15).*exp(-32*kappa^2*(1-gamma).*dNn*(15-(1:15)));
%d2wdZ2(1:11)=d2wdZ2(1:11).*exp(-32*kappa^2*(1-gamma).*dNn*(11-(1:11)));
%Change 1:7/25/2020: I have very heavily weighted the second derivative to
%the left of density peak with an exponential and it allows far larger
%range of parameters with stable execution without sacrificing much
%accuracy. Of course, it is just ad hoc.
end
%If the calculation of the program blows, it means that above calculated
%derivative is unstable or bad. Plese replace with your favourite method of
%calculating this second derivative.


%plot((wnStart:Nn),d2wdZ2,'r',(wnStart:Nn),d2wdZ2A,'g');
%tt
%str=input('Look at comparison of d2wdZ2. Interpolation starts after tt=32');


dZdw(wnStart:Nn)=1.0./dwdZ(wnStart:Nn);
H2(wnStart:Nn)=(Z(wnStart:Nn).^2-1);


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

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

  %wS in point interpolated to the left of the grid.
  %wE in point interpolated to the right of the grid.
  
  wAdjacentB(wnStart:Nn)=0;
  wAdjacentF(wnStart:Nn)=0;
  ZAdjacentB(wnStart:Nn)=0;
  ZAdjacentF(wnStart:Nn)=0;
  dwdZB(wnStart:Nn)=0;
  dwdZF(wnStart:Nn)=0;
  
  
  wAdjacentB(wnStart+1:Nn)=w(wnStart:Nn-1);
  wAdjacentB(wnStart)=wS;
  ZAdjacentB(wnStart+1:Nn)=Z(wnStart:Nn-1);
  ZAdjacentB(wnStart)=Z(wnStart)-dNn;
  dwdZB(wnStart:Nn)=(w(wnStart:Nn)-wAdjacentB(wnStart:Nn))/(dNn);
  
  wAdjacentF(wnStart:Nn-1)=w(wnStart+1:Nn);
  wAdjacentF(Nn)=wE;
  ZAdjacentF(wnStart:Nn-1)=Z(wnStart+1:Nn);
  ZAdjacentF(Nn)=Z(Nn)+dNn;
  dwdZF(wnStart:Nn)=(w(wnStart:Nn)-wAdjacentF(wnStart:Nn))/(-dNn);
  
  DeltaZwDeltaZ(wnStart:Nn)=(ZAdjacentF(wnStart:Nn).*wAdjacentF(wnStart:Nn)-ZAdjacentB(wnStart:Nn).*wAdjacentB(wnStart:Nn))/(2*dNn);
  Term0(wnStart:Nn)=(.5*wAdjacentB(wnStart:Nn)+.5*wAdjacentF(wnStart:Nn)-DeltaZwDeltaZ(wnStart:Nn));
  
  C11(wnStart:Nn)=.707; %dwdZ(wnStart:Nn);%.707000;
  C22=-8*dt;%
  C22=-sqrt(2)/pi/16;%
  C33=2*4*dt*sqrt(2*2*2);%
  C33=sqrt(2)/4/pi;%
%  C22=.5*sigma0*sqrt(dt);%
%  C33=sigma0*sqrt(dt);%
 
  
  w(wnStart:Nn)=w(wnStart:Nn)+wMu0dt(wnStart:Nn)+((Term0(wnStart:Nn)+C11(wnStart:Nn).*d2wdZ2(wnStart:Nn))+ ...
            sign(-Term0(wnStart:Nn)-C11(wnStart:Nn).*d2wdZ2(wnStart:Nn)+dw(wnStart:Nn)).* ...
            sqrt(abs(sign(-Term0(wnStart:Nn)-C11(wnStart:Nn).*d2wdZ2(wnStart:Nn)).* ... 
            (-Term0(wnStart:Nn)-C11(wnStart:Nn).*d2wdZ2(wnStart:Nn)).^2+ ...
            sign(+dw(wnStart:Nn)).*dw2(wnStart:Nn))))+ ...
            -C33.*dwMu0dtdw(wnStart:Nn).*dwdZ(wnStart:Nn).^2 + ...
            +C22.*3*Z(wnStart:Nn).^1.*dZdw(wnStart:Nn).*d2wdZ2(wnStart:Nn)*sigma0^2.*dt + ...
            +0*(1-1/sqrt(3))*dwMu0dtdw(wnStart:Nn).*sigma0^2*dt - ...
            0*.125/2*.125*d2wdZ2(wnStart:Nn).*(w(wnStart:Nn)-Term0(wnStart:Nn))*dt;
        
        %Change 2:7/25/2020:%I have added a last term from quadratic variation between drift
        %and w terms. I am not very sure if we need it or not.

  
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Below are instability correction equations. There is instability in the 
% SDE representation equation roughly between -.5 SD to +.5 SD in normal terms.
%Instability correction Algorithm takes two ends that are here given as
%lower end=w(NnMidl-2) and upper end=w(NnMidh+3).
%In the block below the interpolation starts from lower end.
%1st equation is
%wll(NnMidl-1)=w(NnMidl-2)+(D1wl+CTheta*D2wl*dNn/1+.5*DeltaD1w)*dNn;
%Second equation is: wll(NnMidl)=wll(NnMidl-1)+(D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
%And so on.
%There is a sixth equation that is not written below(since it takes us to 
%the upper interpolation point which is already known) which is 
%wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
%Sum of these six equations add up to the distance between upper and lower
%interpolation points.
%In the sum there are six instances of first derivative which is already 
%calculated at the boundary, six instances of second derivative which is also
%already calculated at lower boundary and 18.0 instances of DeltaD1w which is 
%unknown. This equation calculates the value of one unit of DeltaD1w addition
%DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
% I hope to find a better analytical method to tackle this issue. Friends
% are welcome to play with slightly changed related alternatives to the 
%above interpolation.
%Change 7/27/2020
%I have switched to Lagrange interpolation by taking four points on both 
%sides of instability and using 8th order interpolation of five points in
%mid between the eight points on sides. This works far better and you can
%now use nuemrical 2nd derivative without any modification for a large 
%range of parameters. It truly improves the fit to true density when it 
%works. 

% Good interpolation below
      w(NnMidh+2)=InterpolateOrderN8(8,Z(NnMidh+2),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidl-2),w(NnMidl-3),w(NnMidl-4),w(NnMidl-5));
      w(NnMidh+1)=InterpolateOrderN8(8,Z(NnMidh+1),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidl-2),w(NnMidl-3),w(NnMidl-4),w(NnMidl-5));
      w(NnMidh)=InterpolateOrderN8(8,Z(NnMidh),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidl-2),w(NnMidl-3),w(NnMidl-4),w(NnMidl-5));
      w(NnMidl)=InterpolateOrderN8(8,Z(NnMidl),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidl-2),w(NnMidl-3),w(NnMidl-4),w(NnMidl-5));
      w(NnMidl-1)=InterpolateOrderN8(8,Z(NnMidl-1),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidl-2),w(NnMidl-3),w(NnMidl-4),w(NnMidl-5));
           
           
      
    
% 
% 
% 
%     CTheta=0;
%     D1wl=(11/6*w(NnMidl-2)-3*w(NnMidl-3)+3/2*w(NnMidl-4)-1/3*w(NnMidl-5))/dNn;
%     D2wl=(35/12*w(NnMidl-2) -26/3*w(NnMidl-3)+ 19/2*w(NnMidl-4)-14/3*w(NnMidl-5)+11/12*w(NnMidl-6))/dNn.^2;
%     DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
%     
%     wll(NnMidl-1)=w(NnMidl-2)+  (D1wl+CTheta*D2wl*dNn+.5*DeltaD1w)*dNn;
%     wll(NnMidl)=wll(NnMidl-1)+  (D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
%     wll(NnMidh)=wll(NnMidl)+    (D1wl+CTheta*D2wl*dNn+2.5*DeltaD1w)*dNn;
%     wll(NnMidh+1)=wll(NnMidh)+  (D1wl+CTheta*D2wl*dNn+3.5*DeltaD1w)*dNn;
%     wll(NnMidh+2)=wll(NnMidh+1)+(D1wl+CTheta*D2wl*dNn+4.5*DeltaD1w)*dNn;
%     %wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
% %The last equation is not used since point w(NnMidh+3)is given but it is used in 
% %calculation of distance between two interpolation points. This distance is
% %used to determine the value of DeltaD1w.
% 
%     D1wh=(-11/6*w(NnMidh+3)+3*w(NnMidh+4)-3/2*w(NnMidh+5)+1/3*w(NnMidh+6))/dNn;
%     D2wh=(35/12*w(NnMidh+3) -26/3*w(NnMidh+4)+ 19/2*w(NnMidh+5)-14/3*w(NnMidh+6)+11/12*w(NnMidh+7))/dNn.^2;
%     DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wh*dNn-6.0*CTheta*D2wh*dNn.^2.0/1.0)/18.0/dNn;
% %     DeltaD1w=(w(NnMidh+2)-w(NnMidl-1)-4*D1wh*dNn-16.0*D2wh*dNn.^2.0/1.0)/6.0/dNn;
% 
% 
%     whh(NnMidh+2)=w(NnMidh+3)-(D1wh+CTheta*D2wh*dNn+.5*DeltaD1w)*dNn;
%     whh(NnMidh+1)=w(NnMidh+2)-(D1wh+CTheta*D2wh*dNn+1.5* DeltaD1w)*dNn;
%     whh(NnMidh)=whh(NnMidh+1)-(D1wh+CTheta*D2wh*dNn+2.5* DeltaD1w)*dNn;
%     whh(NnMidl)=whh(NnMidh)-(D1wh+CTheta*D2wh*dNn+3.5* DeltaD1w)*dNn;
%     whh(NnMidl-1)=whh(NnMidl)-(D1wh+CTheta*D2wh*dNn+4.5* DeltaD1w)*dNn;
% 
%     %whh(NnMidl-2)=whh(NnMidl-1)-(D1wh+CTheta*D2wh*dNn/1+5.5* DeltaD1w)*dNn;
% 
% 
%         w(NnMidh+2)=.5*whh(NnMidh+2)+.5*wll(NnMidh+2);
%         w(NnMidh+1)=.5*whh(NnMidh+1)+.5*wll(NnMidh+1);
%         w(NnMidh)=.5*whh(NnMidh)+.5*wll(NnMidh);
%         w(NnMidl)=.5*whh(NnMidl)+.5*wll(NnMidl);
%         w(NnMidl-1)=.5*whh(NnMidl-1)+.5*wll(NnMidl-1);


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

        
        w1(1:Nn-1)=w(1:Nn-1);
        w1(Nn)=w(Nn);
        w2(1:Nn-1)=w(2:Nn);
        w2(Nn)=wE;     
        
w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
     %Change 3:7/25/2020: I have improved zero correction in above.
        w(w<0)=0.0;
        for nn=1:Nn
            if(w(nn)<=0)
                wnStart=nn+1;
           end
        end
        
        %%1st order mean correction. We know the mean in original
        %%coordinates so I shift the density into original coordinates,
        %%apply the mean correction and then transform again into Lamperti
        %%coordinates. Algebra solves two equations in two unknowns.
        %%Equations are Sum_0^N{(Y_w(wnStart:Nn)-Y0)*W0}=1 and the second
        %%equation is Sum_0^N{Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0}=u0
        %%Two unknows are Y0 and W0. u0 is known mean.
        tt;
        u0=theta+(x0-theta)*exp(-kappa*(tt*dt)); %analytic mean of the density
        
        %If you are not using stochastic volatility, replace above with
        %true mean otherwise results would become garbage. 

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

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

toc

ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 
 disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
 TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average


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

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

%Uncomment for fourth order monte carlo

%   
%   YY(1:paths)=YY(1:paths) + ...
%       (YCoeff0(1,1,2,1).*YY(1:paths).^Fp(1,1,2,1)+YCoeff0(1,2,1,1).*YY(1:paths).^Fp(1,2,1,1)+ ...
%       YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1))*dtM + ...
%       (YCoeff0(1,1,3,1).*YY(1:paths).^Fp(1,1,3,1)+YCoeff0(1,2,2,1).*YY(1:paths).^Fp(1,2,2,1)+ ...
%       YCoeff0(2,1,2,1).*YY(1:paths).^Fp(2,1,2,1)+YCoeff0(1,3,1,1).*YY(1:paths).^Fp(1,3,1,1)+ ...
%       YCoeff0(2,2,1,1).*YY(1:paths).^Fp(2,2,1,1)+YCoeff0(3,1,1,1).*YY(1:paths).^Fp(3,1,1,1))*dtM^2 + ...
%       (YCoeff0(1,1,4,1).*YY(1:paths).^Fp(1,1,4,1)+YCoeff0(1,2,3,1).*YY(1:paths).^Fp(1,2,3,1)+ ...
%       YCoeff0(2,1,3,1).*YY(1:paths).^Fp(2,1,3,1)+YCoeff0(1,3,2,1).*YY(1:paths).^Fp(1,3,2,1)+ ...
%       YCoeff0(2,2,2,1).*YY(1:paths).^Fp(2,2,2,1)+YCoeff0(3,1,2,1).*YY(1:paths).^Fp(3,1,2,1)+ ...
%       YCoeff0(1,4,1,1).*YY(1:paths).^Fp(1,4,1,1)+YCoeff0(2,3,1,1).*YY(1:paths).^Fp(2,3,1,1)+ ...
%       YCoeff0(3,2,1,1).*YY(1:paths).^Fp(3,2,1,1)+YCoeff0(4,1,1,1).*YY(1:paths).^Fp(4,1,1,1))*dtM^3 + ...
%        (YCoeff0(1,1,5,1).*YY(1:paths).^Fp(1,1,5,1)+YCoeff0(1,2,4,1).*yy(1:paths).^Fp(1,2,4,1)+ ...
%     YCoeff0(2,1,4,1).*YY(1:paths).^Fp(2,1,4,1)+YCoeff0(2,1,1,1).*yy(1:paths).^Fp(2,1,1,1)+ ...
%      YCoeff0(2,2,3,1).*YY(1:paths).^Fp(2,2,3,1)+ ...
%      YCoeff0(3,1,3,1).*YY(1:paths).^Fp(3,1,3,1)+YCoeff0(1,4,2,1).*YY(1:paths).^Fp(1,4,2,1)+ ...
%      YCoeff0(2,3,2,1).*YY(1:paths).^Fp(2,3,2,1)+YCoeff0(3,2,2,1).*YY(1:paths).^Fp(3,2,2,1)+ ...
%      YCoeff0(4,1,2,1).*YY(1:paths).^Fp(4,1,2,1)+YCoeff0(1,5,1,1).*YY(1:paths).^Fp(1,5,1,1)+ ...
%      YCoeff0(2,4,1,1).*YY(1:paths).^Fp(2,4,1,1)+ ...
%      YCoeff0(3,3,1,1).*YY(1:paths).^Fp(3,3,1,1)+YCoeff0(4,2,1,1).*YY(1:paths).^Fp(4,2,1,1)+ ...
%      YCoeff0(5,1,1,1).*YY(1:paths).^Fp(5,1,1,1))*dtM^4+ ...
%       ((YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
%     (YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
%     YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5+ ...
%     (YCoeff0(1,1,3,2).*YY(1:paths).^Fp(1,1,3,2)+YCoeff0(1,2,2,2).*YY(1:paths).^Fp(1,2,2,2)+ ...
%     YCoeff0(2,1,2,2).*YY(1:paths).^Fp(2,1,2,2)+YCoeff0(1,3,1,2).*YY(1:paths).^Fp(1,3,1,2)+ ...
%     YCoeff0(2,2,1,2).*YY(1:paths).^Fp(2,2,1,2)+YCoeff0(3,1,1,2).*YY(1:paths).^Fp(3,1,1,2)).*dtM^2.5+ ...
%     (YCoeff0(1,1,4,2).*YY(1:paths).^Fp(1,1,4,2)+YCoeff0(1,2,3,2).*YY(1:paths).^Fp(1,2,3,2)+ ...
%     YCoeff0(2,1,3,2).*YY(1:paths).^Fp(2,1,3,2)+YCoeff0(1,3,2,2).*YY(1:paths).^Fp(1,3,2,2)+ ...
%     YCoeff0(2,2,2,2).*YY(1:paths).^Fp(2,2,2,2)+ YCoeff0(3,1,2,2).*YY(1:paths).^Fp(3,1,2,2)+ ...
%     YCoeff0(1,4,1,2).*YY(1:paths).^Fp(1,4,1,2)+YCoeff0(2,3,1,2).*YY(1:paths).^Fp(2,3,1,2)+ ...
%     YCoeff0(3,2,1,2).*YY(1:paths).^Fp(3,2,1,2)+YCoeff0(4,1,1,2).*YY(1:paths).^Fp(4,1,1,2)).*dtM^3.5) .*HermiteP1(2,1:paths) + ...
%     ((YCoeff0(1,1,1,3).*YY(1:paths).^Fp(1,1,1,3) *dtM) + ...
%     (YCoeff0(1,1,2,3).*YY(1:paths).^Fp(1,1,2,3)+YCoeff0(1,2,1,3).*YY(1:paths).^Fp(1,2,1,3)+ ...
%     YCoeff0(2,1,1,3).*YY(1:paths).^Fp(2,1,1,3)).*dtM^2+ ...
%     (YCoeff0(1,1,3,3).*YY(1:paths).^Fp(1,1,3,3)+YCoeff0(1,2,2,3).*YY(1:paths).^Fp(1,2,2,3)+ ...
%     YCoeff0(2,1,2,3).*YY(1:paths).^Fp(2,1,2,3) + YCoeff0(1,3,1,3).*YY(1:paths).^Fp(1,3,1,3)+ ...
%     YCoeff0(2,2,1,3).*YY(1:paths).^Fp(2,2,1,3)+YCoeff0(3,1,1,3).*YY(1:paths).^Fp(3,1,1,3)).*dtM^3).*HermiteP1(3,1:paths) + ...
%     ((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )+ ...
%     (YCoeff0(1,1,2,4).*YY(1:paths).^Fp(1,1,2,4)+YCoeff0(1,2,1,4).*YY(1:paths).^Fp(1,2,1,4)+ ...
%     YCoeff0(2,1,1,4).*YY(1:paths).^Fp(2,1,1,4))*dtM^2.5).*HermiteP1(4,1:paths) + ...
%     (YCoeff0(1,1,1,5).*YY(1:paths).^Fp(1,1,1,5)*dtM^2.0).*HermiteP1(5,1:paths);
%     
 
     YYMean(tt)=YYMean(tt)+sum(YY(1:paths))/paths;
 end
 YY(YY<0)=0;
 disp('Original process average from monte carlo');
 MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
 disp('Original process average from our simulation');
 ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 
 disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
 TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
 
 
 MaxCutOff=30;
 NoOfBins=round(300*gamma^2*4*sigma0);%Decrease the number of bins if the graph is too 
 [YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
 plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 %plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');
 
 title(sprintf('Ito-Hermite Method VS Monte Carlo Density, x0 = %.2f,theta=%.2f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f', x0,theta,kappa,gamma,sigma0,T));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
 %title('text',sprintf('x0 = %f', x0),sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T))
 %xlabel(a)
 %title('Graphical Comparison of Ito-Hermite Method and Monte Carlo Method Density')
 legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
 
 str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
 

end
Here is the output from the program.
Original process average from monte carlo
MCMean =
   0.249506263266265
Original process average from our simulation
ItoHermiteMean =
   0.249432851976764
true Mean only applicble to standard SV mean reverting type models otherwise disregard
TrueMean =
   0.250251596970927
IndexMax =
   676

Here is the output graph.

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

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

July 28th, 2020, 9:24 am

Here are notes from today's work.

1. I interpolated first and second derivatives in my last posted program yesterday and that was actually sharply deteriorating the results. Please just interpolate the SDE variable only using Lagrange interpolation and use numerical derivatives just as such and decrease the step size if the program creates problems.

2. I have used in the program variable Term0 which is given as
[$].5 w_{n+1}+.5 w_{n-1} - .5 \frac{Z_{n+1} w_{n+1} - Z_{n-1} w_{n-1}}{\Delta Z}[$]

but this is exactly identical with [$]Z\frac{dw}{dZ}[$] and actually
[$]Z\frac{dw}{dZ}= .5 w_{n+1}+.5 w_{n-1} - .5 \frac{Z_{n+1} w_{n+1} - Z_{n-1} w_{n-1}}{\Delta Z}[$]
so in the current program, I have just used Term0=[$]Z\frac{dw}{dZ}[$]

3. I have used the following suggestive rough  values of constants used in the program
  C11(wnStart:Nn)=.707;
  C22=-1/pi/8;
  C33=1.2/4/pi;

4.  C11 is coefficient of 2nd derivative and my wild guess is that somehow it might be getting a division with [$]\sqrt{2}[$] due to the coefficient in normalization of second hermite polynomial.
C33 is very very close to .1 which has the possibility of being equal to [$]\frac{\Delta Z}{2}[$] because in our program  [$]\Delta Z=.2[$]

5. Though new interpolation with Lagrange polynomials greatly improved the program, I still think there is need for improvement and a better interpolation that preserves derivatives would further increase accuracy.

6. I have continued to miss a third derivative term in all of these programs since even second derivatives just recently became stable. In many cases, adding this term would further improve accuracy of the program.

Here is my new improved program with derivatives interpolations deleted and new constants
function [TrueMean,ItoHermiteMean,y_w,fy_w,wnStart] = FPERevisited044Wilmott16A()

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

%I have not directly simulated the SDE but simulated the transformed 
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.
%In this version, the mean correction has been disabled. If you are
%simulating standard mean-reverting SV type SDEs, please enable the mean 
%correction by uncommenting the appropriate line in the body of code below.

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


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

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

x0=.10;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.75;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=2.0;%.950;   %mean reversion parameter.
theta=.10;%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;
                
w(1:Nn)=x0^(1-gamma)/(1-gamma);

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

Z
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);

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

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

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

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

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

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

for k = 0 : (OrderA)
    for m = 0:k
        l4 = k - m + 1;
        for n = 0 : m
            l3 = m - n + 1;
            for j = 0:n
                l2 = n - j + 1;
                l1 = j + 1;
                %Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
                Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                
                YCoeff0(l1,l2,l3,l4) =YCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
                YqCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
            end
        end
    end
end

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

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

for tt=1:Tt
    
     yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
    
    [wMu0dt,dwMu0dtdw,d2wMu0dtdw2,d3wMu0dtdw3,c1,c2,c3,c4] = CalculateDriftAndVolOrder4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);

    wPrev=w;
    dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn) ;% ...
    dw2(wnStart:Nn)=dw(wnStart:Nn).^2;
    w(isnan(w)==1)=0;
    wMu0dt(isnan(wMu0dt)==1)=0;  
    wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
   
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



d2wdZ2(wnStart:Nn)=d2wdZ2(wnStart:Nn)+ d2wMu0dtdw2(wnStart:Nn).*dwdZ(wnStart:Nn).^2+dwMu0dtdw(wnStart:Nn).*d2wdZ2(wnStart:Nn)+.33*(Z(wnStart:Nn)).*sigma0^2*dt^2;;%+ ...


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


if (tt>16)
d2wdZ2=d2wdZ2A;

d2wdZ2(1:8)=d2wdZ2(1:8).*exp(-32*kappa^2*(1-gamma).*dNn*(8-(1:8)));
%Change 1:7/25/2020: I have very heavily weighted the second derivative to
%the left of density peak with an exponential and it allows far larger
%range of parameters with stable execution without sacrificing much
%accuracy. Of course, it is just ad hoc.
end
%If the calculation of the program blows, it means that above calculated
%derivative is unstable or bad. Plese replace with your favourite method of
%calculating this second derivative.


%plot((wnStart:Nn),d2wdZ2,'r',(wnStart:Nn),d2wdZ2A,'g');
%tt
%str=input('Look at comparison of d2wdZ2. Interpolation starts after tt=32');


dZdw(wnStart:Nn)=1.0./dwdZ(wnStart:Nn);
H2(wnStart:Nn)=(Z(wnStart:Nn).^2-1);


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

  
%   [wS] = InterpolateOrderN6(6,Z(1)-dNn,Z(1),Z(2),Z(3),Z(4),Z(5),Z(6),w(1),w(2),w(3),w(4),w(5),w(6));
%   [wE] = InterpolateOrderN6(6,Z(Nn)+dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));
% 
%   %wS in point interpolated to the left of the grid.
%   %wE in point interpolated to the right of the grid.
%   
%   wAdjacentB(wnStart:Nn)=0;
%   wAdjacentF(wnStart:Nn)=0;
%   ZAdjacentB(wnStart:Nn)=0;
%   ZAdjacentF(wnStart:Nn)=0;
%   dwdZB(wnStart:Nn)=0;
%   dwdZF(wnStart:Nn)=0;
%   
%   
%   wAdjacentB(wnStart+1:Nn)=w(wnStart:Nn-1);
%   wAdjacentB(wnStart)=wS;
%   ZAdjacentB(wnStart+1:Nn)=Z(wnStart:Nn-1);
%   ZAdjacentB(wnStart)=Z(wnStart)-dNn;
%   dwdZB(wnStart:Nn)=(w(wnStart:Nn)-wAdjacentB(wnStart:Nn))/(dNn);
%   
%   wAdjacentF(wnStart:Nn-1)=w(wnStart+1:Nn);
%   wAdjacentF(Nn)=wE;
%   ZAdjacentF(wnStart:Nn-1)=Z(wnStart+1:Nn);
%   ZAdjacentF(Nn)=Z(Nn)+dNn;
%   dwdZF(wnStart:Nn)=(w(wnStart:Nn)-wAdjacentF(wnStart:Nn))/(-dNn);
%   
%   DeltaZwDeltaZ(wnStart:Nn)=(ZAdjacentF(wnStart:Nn).*wAdjacentF(wnStart:Nn)-ZAdjacentB(wnStart:Nn).*wAdjacentB(wnStart:Nn))/(2*dNn);
%   Term0(wnStart:Nn)=(.5*wAdjacentB(wnStart:Nn)+.5*wAdjacentF(wnStart:Nn)-DeltaZwDeltaZ(wnStart:Nn));
  
    %if(tt>32)
  %Term0(wnStart:Nn)=-(Z(wnStart:Nn).*dwdZ(wnStart:Nn));
  Term0(wnStart:Nn)=-(Z(wnStart:Nn).*dwdZ(wnStart:Nn));
  %end
  C11(wnStart:Nn)=.707; %dwdZ(wnStart:Nn);%.707000;
  
  C22=-1/pi/8;%

  C33=1.2/4/pi;%

  
  w(wnStart:Nn)=w(wnStart:Nn)+wMu0dt(wnStart:Nn)+((Term0(wnStart:Nn)+C11(wnStart:Nn).*d2wdZ2(wnStart:Nn))+ ...
            sign(-Term0(wnStart:Nn)-C11(wnStart:Nn).*d2wdZ2(wnStart:Nn)+dw(wnStart:Nn)).* ...
            sqrt(abs(sign(-Term0(wnStart:Nn)-C11(wnStart:Nn).*d2wdZ2(wnStart:Nn)).* ... 
            (-Term0(wnStart:Nn)-C11(wnStart:Nn).*d2wdZ2(wnStart:Nn)).^2+ ...
            sign(+dw(wnStart:Nn)).*dw2(wnStart:Nn))))+ ...
            -C33.*dwMu0dtdw(wnStart:Nn).*dwdZ(wnStart:Nn).^2 + ...
            +C22.*3*Z(wnStart:Nn).^1.*dZdw(wnStart:Nn).*d2wdZ2(wnStart:Nn)*sigma0^2.*dt + ...
            +0*(1-1/sqrt(3))*dwMu0dtdw(wnStart:Nn).*sigma0^2*dt - ...
            0*.125/2*.125*d2wdZ2(wnStart:Nn).*(w(wnStart:Nn)-Term0(wnStart:Nn))*dt;
        
        %Change 2:7/25/2020:%I have added a last term from quadratic variation between drift
        %and w terms. I am not very sure if we need it or not.

  
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Below are instability correction equations. There is instability in the 
% SDE representation equation roughly between -.5 SD to +.5 SD in normal terms.
%Instability correction Algorithm takes two ends that are here given as
%lower end=w(NnMidl-2) and upper end=w(NnMidh+3).
%In the block below the interpolation starts from lower end.
%1st equation is
%wll(NnMidl-1)=w(NnMidl-2)+(D1wl+CTheta*D2wl*dNn/1+.5*DeltaD1w)*dNn;
%Second equation is: wll(NnMidl)=wll(NnMidl-1)+(D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
%And so on.
%There is a sixth equation that is not written below(since it takes us to 
%the upper interpolation point which is already known) which is 
%wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
%Sum of these six equations add up to the distance between upper and lower
%interpolation points.
%In the sum there are six instances of first derivative which is already 
%calculated at the boundary, six instances of second derivative which is also
%already calculated at lower boundary and 18.0 instances of DeltaD1w which is 
%unknown. This equation calculates the value of one unit of DeltaD1w addition
%DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
% I hope to find a better analytical method to tackle this issue. Friends
% are welcome to play with slightly changed related alternatives to the 
%above interpolation.
%Change 7/27/2020
%I have switched to Lagrange interpolation by taking four points on both 
%sides of instability and using 8th order interpolation of five points in
%mid between the eight points on sides. This works far better and you can
%now use nuemrical 2nd derivative without any modification for a large 
%range of parameters. It truly improves the fit to true density when it 
%works. 

% Good interpolation below
      w(NnMidh+2)=InterpolateOrderN8(8,Z(NnMidh+2),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidl-2),w(NnMidl-3),w(NnMidl-4),w(NnMidl-5));
      w(NnMidh+1)=InterpolateOrderN8(8,Z(NnMidh+1),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidl-2),w(NnMidl-3),w(NnMidl-4),w(NnMidl-5));
      w(NnMidh)=InterpolateOrderN8(8,Z(NnMidh),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidl-2),w(NnMidl-3),w(NnMidl-4),w(NnMidl-5));
      w(NnMidl)=InterpolateOrderN8(8,Z(NnMidl),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidl-2),w(NnMidl-3),w(NnMidl-4),w(NnMidl-5));
      w(NnMidl-1)=InterpolateOrderN8(8,Z(NnMidl-1),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidl-2),w(NnMidl-3),w(NnMidl-4),w(NnMidl-5));
           
           
      
    
% 
% 
% 
%     CTheta=0;
%     D1wl=(11/6*w(NnMidl-2)-3*w(NnMidl-3)+3/2*w(NnMidl-4)-1/3*w(NnMidl-5))/dNn;
%     D2wl=(35/12*w(NnMidl-2) -26/3*w(NnMidl-3)+ 19/2*w(NnMidl-4)-14/3*w(NnMidl-5)+11/12*w(NnMidl-6))/dNn.^2;
%     DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
%     
%     wll(NnMidl-1)=w(NnMidl-2)+  (D1wl+CTheta*D2wl*dNn+.5*DeltaD1w)*dNn;
%     wll(NnMidl)=wll(NnMidl-1)+  (D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
%     wll(NnMidh)=wll(NnMidl)+    (D1wl+CTheta*D2wl*dNn+2.5*DeltaD1w)*dNn;
%     wll(NnMidh+1)=wll(NnMidh)+  (D1wl+CTheta*D2wl*dNn+3.5*DeltaD1w)*dNn;
%     wll(NnMidh+2)=wll(NnMidh+1)+(D1wl+CTheta*D2wl*dNn+4.5*DeltaD1w)*dNn;
%     %wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
% %The last equation is not used since point w(NnMidh+3)is given but it is used in 
% %calculation of distance between two interpolation points. This distance is
% %used to determine the value of DeltaD1w.
% 
%     D1wh=(-11/6*w(NnMidh+3)+3*w(NnMidh+4)-3/2*w(NnMidh+5)+1/3*w(NnMidh+6))/dNn;
%     D2wh=(35/12*w(NnMidh+3) -26/3*w(NnMidh+4)+ 19/2*w(NnMidh+5)-14/3*w(NnMidh+6)+11/12*w(NnMidh+7))/dNn.^2;
%     DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wh*dNn-6.0*CTheta*D2wh*dNn.^2.0/1.0)/18.0/dNn;
% %     DeltaD1w=(w(NnMidh+2)-w(NnMidl-1)-4*D1wh*dNn-16.0*D2wh*dNn.^2.0/1.0)/6.0/dNn;
% 
% 
%     whh(NnMidh+2)=w(NnMidh+3)-(D1wh+CTheta*D2wh*dNn+.5*DeltaD1w)*dNn;
%     whh(NnMidh+1)=w(NnMidh+2)-(D1wh+CTheta*D2wh*dNn+1.5* DeltaD1w)*dNn;
%     whh(NnMidh)=whh(NnMidh+1)-(D1wh+CTheta*D2wh*dNn+2.5* DeltaD1w)*dNn;
%     whh(NnMidl)=whh(NnMidh)-(D1wh+CTheta*D2wh*dNn+3.5* DeltaD1w)*dNn;
%     whh(NnMidl-1)=whh(NnMidl)-(D1wh+CTheta*D2wh*dNn+4.5* DeltaD1w)*dNn;
% 
%     %whh(NnMidl-2)=whh(NnMidl-1)-(D1wh+CTheta*D2wh*dNn/1+5.5* DeltaD1w)*dNn;
% 
% 
%         w(NnMidh+2)=.5*whh(NnMidh+2)+.5*wll(NnMidh+2);
%         w(NnMidh+1)=.5*whh(NnMidh+1)+.5*wll(NnMidh+1);
%         w(NnMidh)=.5*whh(NnMidh)+.5*wll(NnMidh);
%         w(NnMidl)=.5*whh(NnMidl)+.5*wll(NnMidl);
%         w(NnMidl-1)=.5*whh(NnMidl-1)+.5*wll(NnMidl-1);


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

        
        w1(1:Nn-1)=w(1:Nn-1);
        w1(Nn)=w(Nn);
        w2(1:Nn-1)=w(2:Nn);
        w2(Nn)=wE;     
        
w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
     %Change 3:7/25/2020: I have improved zero correction in above.
        w(w<0)=0.0;
        for nn=1:Nn
            if(w(nn)<=0)
                wnStart=nn+1;
           end
        end
        
        %%1st order mean correction. We know the mean in original
        %%coordinates so I shift the density into original coordinates,
        %%apply the mean correction and then transform again into Lamperti
        %%coordinates. Algebra solves two equations in two unknowns.
        %%Equations are Sum_0^N{(Y_w(wnStart:Nn)-Y0)*W0}=1 and the second
        %%equation is Sum_0^N{Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0}=u0
        %%Two unknows are Y0 and W0. u0 is known mean.
        tt;
        u0=theta+(x0-theta)*exp(-kappa*(tt*dt)); %analytic mean of the density
        
        %If you are not using stochastic volatility, replace above with
        %true mean otherwise results would become garbage. 

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

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

toc

ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 
 disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
 TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average


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

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

%Uncomment for fourth order monte carlo

%   
%   YY(1:paths)=YY(1:paths) + ...
%       (YCoeff0(1,1,2,1).*YY(1:paths).^Fp(1,1,2,1)+YCoeff0(1,2,1,1).*YY(1:paths).^Fp(1,2,1,1)+ ...
%       YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1))*dtM + ...
%       (YCoeff0(1,1,3,1).*YY(1:paths).^Fp(1,1,3,1)+YCoeff0(1,2,2,1).*YY(1:paths).^Fp(1,2,2,1)+ ...
%       YCoeff0(2,1,2,1).*YY(1:paths).^Fp(2,1,2,1)+YCoeff0(1,3,1,1).*YY(1:paths).^Fp(1,3,1,1)+ ...
%       YCoeff0(2,2,1,1).*YY(1:paths).^Fp(2,2,1,1)+YCoeff0(3,1,1,1).*YY(1:paths).^Fp(3,1,1,1))*dtM^2 + ...
%       (YCoeff0(1,1,4,1).*YY(1:paths).^Fp(1,1,4,1)+YCoeff0(1,2,3,1).*YY(1:paths).^Fp(1,2,3,1)+ ...
%       YCoeff0(2,1,3,1).*YY(1:paths).^Fp(2,1,3,1)+YCoeff0(1,3,2,1).*YY(1:paths).^Fp(1,3,2,1)+ ...
%       YCoeff0(2,2,2,1).*YY(1:paths).^Fp(2,2,2,1)+YCoeff0(3,1,2,1).*YY(1:paths).^Fp(3,1,2,1)+ ...
%       YCoeff0(1,4,1,1).*YY(1:paths).^Fp(1,4,1,1)+YCoeff0(2,3,1,1).*YY(1:paths).^Fp(2,3,1,1)+ ...
%       YCoeff0(3,2,1,1).*YY(1:paths).^Fp(3,2,1,1)+YCoeff0(4,1,1,1).*YY(1:paths).^Fp(4,1,1,1))*dtM^3 + ...
%        (YCoeff0(1,1,5,1).*YY(1:paths).^Fp(1,1,5,1)+YCoeff0(1,2,4,1).*yy(1:paths).^Fp(1,2,4,1)+ ...
%     YCoeff0(2,1,4,1).*YY(1:paths).^Fp(2,1,4,1)+YCoeff0(2,1,1,1).*yy(1:paths).^Fp(2,1,1,1)+ ...
%      YCoeff0(2,2,3,1).*YY(1:paths).^Fp(2,2,3,1)+ ...
%      YCoeff0(3,1,3,1).*YY(1:paths).^Fp(3,1,3,1)+YCoeff0(1,4,2,1).*YY(1:paths).^Fp(1,4,2,1)+ ...
%      YCoeff0(2,3,2,1).*YY(1:paths).^Fp(2,3,2,1)+YCoeff0(3,2,2,1).*YY(1:paths).^Fp(3,2,2,1)+ ...
%      YCoeff0(4,1,2,1).*YY(1:paths).^Fp(4,1,2,1)+YCoeff0(1,5,1,1).*YY(1:paths).^Fp(1,5,1,1)+ ...
%      YCoeff0(2,4,1,1).*YY(1:paths).^Fp(2,4,1,1)+ ...
%      YCoeff0(3,3,1,1).*YY(1:paths).^Fp(3,3,1,1)+YCoeff0(4,2,1,1).*YY(1:paths).^Fp(4,2,1,1)+ ...
%      YCoeff0(5,1,1,1).*YY(1:paths).^Fp(5,1,1,1))*dtM^4+ ...
%       ((YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
%     (YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
%     YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5+ ...
%     (YCoeff0(1,1,3,2).*YY(1:paths).^Fp(1,1,3,2)+YCoeff0(1,2,2,2).*YY(1:paths).^Fp(1,2,2,2)+ ...
%     YCoeff0(2,1,2,2).*YY(1:paths).^Fp(2,1,2,2)+YCoeff0(1,3,1,2).*YY(1:paths).^Fp(1,3,1,2)+ ...
%     YCoeff0(2,2,1,2).*YY(1:paths).^Fp(2,2,1,2)+YCoeff0(3,1,1,2).*YY(1:paths).^Fp(3,1,1,2)).*dtM^2.5+ ...
%     (YCoeff0(1,1,4,2).*YY(1:paths).^Fp(1,1,4,2)+YCoeff0(1,2,3,2).*YY(1:paths).^Fp(1,2,3,2)+ ...
%     YCoeff0(2,1,3,2).*YY(1:paths).^Fp(2,1,3,2)+YCoeff0(1,3,2,2).*YY(1:paths).^Fp(1,3,2,2)+ ...
%     YCoeff0(2,2,2,2).*YY(1:paths).^Fp(2,2,2,2)+ YCoeff0(3,1,2,2).*YY(1:paths).^Fp(3,1,2,2)+ ...
%     YCoeff0(1,4,1,2).*YY(1:paths).^Fp(1,4,1,2)+YCoeff0(2,3,1,2).*YY(1:paths).^Fp(2,3,1,2)+ ...
%     YCoeff0(3,2,1,2).*YY(1:paths).^Fp(3,2,1,2)+YCoeff0(4,1,1,2).*YY(1:paths).^Fp(4,1,1,2)).*dtM^3.5) .*HermiteP1(2,1:paths) + ...
%     ((YCoeff0(1,1,1,3).*YY(1:paths).^Fp(1,1,1,3) *dtM) + ...
%     (YCoeff0(1,1,2,3).*YY(1:paths).^Fp(1,1,2,3)+YCoeff0(1,2,1,3).*YY(1:paths).^Fp(1,2,1,3)+ ...
%     YCoeff0(2,1,1,3).*YY(1:paths).^Fp(2,1,1,3)).*dtM^2+ ...
%     (YCoeff0(1,1,3,3).*YY(1:paths).^Fp(1,1,3,3)+YCoeff0(1,2,2,3).*YY(1:paths).^Fp(1,2,2,3)+ ...
%     YCoeff0(2,1,2,3).*YY(1:paths).^Fp(2,1,2,3) + YCoeff0(1,3,1,3).*YY(1:paths).^Fp(1,3,1,3)+ ...
%     YCoeff0(2,2,1,3).*YY(1:paths).^Fp(2,2,1,3)+YCoeff0(3,1,1,3).*YY(1:paths).^Fp(3,1,1,3)).*dtM^3).*HermiteP1(3,1:paths) + ...
%     ((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )+ ...
%     (YCoeff0(1,1,2,4).*YY(1:paths).^Fp(1,1,2,4)+YCoeff0(1,2,1,4).*YY(1:paths).^Fp(1,2,1,4)+ ...
%     YCoeff0(2,1,1,4).*YY(1:paths).^Fp(2,1,1,4))*dtM^2.5).*HermiteP1(4,1:paths) + ...
%     (YCoeff0(1,1,1,5).*YY(1:paths).^Fp(1,1,1,5)*dtM^2.0).*HermiteP1(5,1:paths);
%     
 
     YYMean(tt)=YYMean(tt)+sum(YY(1:paths))/paths;
 end
 YY(YY<0)=0;
 disp('Original process average from monte carlo');
 MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
 disp('Original process average from our simulation');
 ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 
 disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
 TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
 
 
 MaxCutOff=30;
 NoOfBins=round(300*gamma^2*4*sigma0);%Decrease the number of bins if the graph is too 
 [YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
 plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 %plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');
 
 title(sprintf('Ito-Hermite Method VS Monte Carlo Density, x0 = %.2f,theta=%.2f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f', x0,theta,kappa,gamma,sigma0,T));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
 %title('text',sprintf('x0 = %f', x0),sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T))
 %xlabel(a)
 %title('Graphical Comparison of Ito-Hermite Method and Monte Carlo Method Density')
 legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
 
 str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
 

end
Original process average from monte carlo
MCMean =
   0.099666492688799
Original process average from our simulation
ItoHermiteMean =
   0.097592958223771
true Mean only applicble to standard SV mean reverting type models otherwise disregard
TrueMean =
   0.100000000000000
IndexMax =
   609

and has the following graph
Image
 
User avatar
Amin
Topic Author
Posts: 2536
Joined: July 14th, 2002, 3:00 am

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

July 28th, 2020, 12:52 pm

I used spline interpolation function of matlab for interpolation of the 'instability data points' in the mid of the density with vastly improved results. Here I post the new matlab function.
function [TrueMean,ItoHermiteMean,y_w,fy_w,wnStart] = FPERevisited044Wilmott16A()

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

%I have not directly simulated the SDE but simulated the transformed 
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.
%In this version, the mean correction has been disabled. If you are
%simulating standard mean-reverting SV type SDEs, please enable the mean 
%correction by uncommenting the appropriate line in the body of code below.

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


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

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

x0=.10;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.65;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=4.0;%.950;   %mean reversion parameter.
theta=.10;%mean reversion target
sigma0=.75;%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;
                
w(1:Nn)=x0^(1-gamma)/(1-gamma);

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

Z
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);

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

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

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

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

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

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

for k = 0 : (OrderA)
    for m = 0:k
        l4 = k - m + 1;
        for n = 0 : m
            l3 = m - n + 1;
            for j = 0:n
                l2 = n - j + 1;
                l1 = j + 1;
                %Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
                Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                
                YCoeff0(l1,l2,l3,l4) =YCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
                YqCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
            end
        end
    end
end

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

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

for tt=1:Tt
    
     yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
    
    [wMu0dt,dwMu0dtdw,d2wMu0dtdw2,d3wMu0dtdw3,c1,c2,c3,c4] = CalculateDriftAndVolOrder4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);

    wPrev=w;
    dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn) ;% ...
    dw2(wnStart:Nn)=dw(wnStart:Nn).^2;
    w(isnan(w)==1)=0;
    wMu0dt(isnan(wMu0dt)==1)=0;  
    wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
   
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



d2wdZ2(wnStart:Nn)=d2wdZ2(wnStart:Nn)+ d2wMu0dtdw2(wnStart:Nn).*dwdZ(wnStart:Nn).^2+dwMu0dtdw(wnStart:Nn).*d2wdZ2(wnStart:Nn)+.33*(Z(wnStart:Nn)).*sigma0^2*dt^2;;%+ ...


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


if (tt>16)
d2wdZ2=d2wdZ2A;

d2wdZ2(1:8)=d2wdZ2(1:8).*exp(-32*kappa^2*(1-gamma).*dNn*(8-(1:8)));
%Change 1:7/25/2020: I have very heavily weighted the second derivative to
%the left of density peak with an exponential and it allows far larger
%range of parameters with stable execution without sacrificing much
%accuracy. Of course, it is just ad hoc.
end
%If the calculation of the program blows, it means that above calculated
%derivative is unstable or bad. Plese replace with your favourite method of
%calculating this second derivative.


%plot((wnStart:Nn),d2wdZ2,'r',(wnStart:Nn),d2wdZ2A,'g');
%tt
%str=input('Look at comparison of d2wdZ2. Interpolation starts after tt=32');


dZdw(wnStart:Nn)=1.0./dwdZ(wnStart:Nn);
H2(wnStart:Nn)=(Z(wnStart:Nn).^2-1);


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

  
%   [wS] = InterpolateOrderN6(6,Z(1)-dNn,Z(1),Z(2),Z(3),Z(4),Z(5),Z(6),w(1),w(2),w(3),w(4),w(5),w(6));
%   [wE] = InterpolateOrderN6(6,Z(Nn)+dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));
% 
%   %wS in point interpolated to the left of the grid.
%   %wE in point interpolated to the right of the grid.
%   
%   wAdjacentB(wnStart:Nn)=0;
%   wAdjacentF(wnStart:Nn)=0;
%   ZAdjacentB(wnStart:Nn)=0;
%   ZAdjacentF(wnStart:Nn)=0;
%   dwdZB(wnStart:Nn)=0;
%   dwdZF(wnStart:Nn)=0;
%   
%   
%   wAdjacentB(wnStart+1:Nn)=w(wnStart:Nn-1);
%   wAdjacentB(wnStart)=wS;
%   ZAdjacentB(wnStart+1:Nn)=Z(wnStart:Nn-1);
%   ZAdjacentB(wnStart)=Z(wnStart)-dNn;
%   dwdZB(wnStart:Nn)=(w(wnStart:Nn)-wAdjacentB(wnStart:Nn))/(dNn);
%   
%   wAdjacentF(wnStart:Nn-1)=w(wnStart+1:Nn);
%   wAdjacentF(Nn)=wE;
%   ZAdjacentF(wnStart:Nn-1)=Z(wnStart+1:Nn);
%   ZAdjacentF(Nn)=Z(Nn)+dNn;
%   dwdZF(wnStart:Nn)=(w(wnStart:Nn)-wAdjacentF(wnStart:Nn))/(-dNn);
%   
%   DeltaZwDeltaZ(wnStart:Nn)=(ZAdjacentF(wnStart:Nn).*wAdjacentF(wnStart:Nn)-ZAdjacentB(wnStart:Nn).*wAdjacentB(wnStart:Nn))/(2*dNn);
%   Term0(wnStart:Nn)=(.5*wAdjacentB(wnStart:Nn)+.5*wAdjacentF(wnStart:Nn)-DeltaZwDeltaZ(wnStart:Nn));
  
    %if(tt>32)
  %Term0(wnStart:Nn)=-(Z(wnStart:Nn).*dwdZ(wnStart:Nn));
  Term0(wnStart:Nn)=-(Z(wnStart:Nn).*dwdZ(wnStart:Nn));
  %end
  C11(wnStart:Nn)=.707; %dwdZ(wnStart:Nn);%.707000;
  
  C22=-1/pi/8;%

  C33=1.2/4/pi;%

  
  w(wnStart:Nn)=w(wnStart:Nn)+wMu0dt(wnStart:Nn)+((Term0(wnStart:Nn)+C11(wnStart:Nn).*d2wdZ2(wnStart:Nn))+ ...
            sign(-Term0(wnStart:Nn)-C11(wnStart:Nn).*d2wdZ2(wnStart:Nn)+dw(wnStart:Nn)).* ...
            sqrt(abs(sign(-Term0(wnStart:Nn)-C11(wnStart:Nn).*d2wdZ2(wnStart:Nn)).* ... 
            (-Term0(wnStart:Nn)-C11(wnStart:Nn).*d2wdZ2(wnStart:Nn)).^2+ ...
            sign(+dw(wnStart:Nn)).*dw2(wnStart:Nn))))+ ...
            -C33.*dwMu0dtdw(wnStart:Nn).*dwdZ(wnStart:Nn).^2 + ...
            +C22.*3*Z(wnStart:Nn).^1.*dZdw(wnStart:Nn).*d2wdZ2(wnStart:Nn)*sigma0^2.*dt + ...
            +0*(1-1/sqrt(3))*dwMu0dtdw(wnStart:Nn).*sigma0^2*dt - ...
            0*.125/2*.125*d2wdZ2(wnStart:Nn).*(w(wnStart:Nn)-Term0(wnStart:Nn))*dt;
        
        %Change 2:7/25/2020:%I have added a last term from quadratic variation between drift
        %and w terms. I am not very sure if we need it or not.

  
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Below are instability correction equations. There is instability in the 
% SDE representation equation roughly between -.5 SD to +.5 SD in normal terms.
%Instability correction Algorithm takes two ends that are here given as
%lower end=w(NnMidl-2) and upper end=w(NnMidh+3).
%In the block below the interpolation starts from lower end.
%1st equation is
%wll(NnMidl-1)=w(NnMidl-2)+(D1wl+CTheta*D2wl*dNn/1+.5*DeltaD1w)*dNn;
%Second equation is: wll(NnMidl)=wll(NnMidl-1)+(D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
%And so on.
%There is a sixth equation that is not written below(since it takes us to 
%the upper interpolation point which is already known) which is 
%wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
%Sum of these six equations add up to the distance between upper and lower
%interpolation points.
%In the sum there are six instances of first derivative which is already 
%calculated at the boundary, six instances of second derivative which is also
%already calculated at lower boundary and 18.0 instances of DeltaD1w which is 
%unknown. This equation calculates the value of one unit of DeltaD1w addition
%DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
% I hope to find a better analytical method to tackle this issue. Friends
% are welcome to play with slightly changed related alternatives to the 
%above interpolation.
%Change 7/27/2020
%I have switched to Lagrange interpolation by taking four points on both 
%sides of instability and using 8th order interpolation of five points in
%mid between the eight points on sides. This works far better and you can
%now use nuemrical 2nd derivative without any modification for a large 
%range of parameters. It truly improves the fit to true density when it 
%works. 

% % Good interpolation below
%       w(NnMidh+2)=InterpolateOrderN8(8,Z(NnMidh+2),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidl-2),w(NnMidl-3),w(NnMidl-4),w(NnMidl-5));
%       w(NnMidh+1)=InterpolateOrderN8(8,Z(NnMidh+1),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidl-2),w(NnMidl-3),w(NnMidl-4),w(NnMidl-5));
%       w(NnMidh)=InterpolateOrderN8(8,Z(NnMidh),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidl-2),w(NnMidl-3),w(NnMidl-4),w(NnMidl-5));
%       w(NnMidl)=InterpolateOrderN8(8,Z(NnMidl),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidl-2),w(NnMidl-3),w(NnMidl-4),w(NnMidl-5));
%       w(NnMidl-1)=InterpolateOrderN8(8,Z(NnMidl-1),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidl-2),Z(NnMidl-3),Z(NnMidl-4),Z(NnMidl-5),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidl-2),w(NnMidl-3),w(NnMidl-4),w(NnMidl-5));
%            

      
      %yq = makima(x,y,xq)
      
      
      Zi=[Z(NnMidl-5),Z(NnMidl-4),Z(NnMidl-3),Z(NnMidl-2),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6)];
      wi=[w(NnMidl-5),w(NnMidl-4),w(NnMidl-3),w(NnMidl-2),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6)];
      Zq=[Z(NnMidl-1),Z(NnMidl),Z(NnMidh),Z(NnMidh+1),Z(NnMidh+2)];
      wq=spline(Zi,wi,Zq);
      w(NnMidl-1:NnMidh+2)=wq(1:5);
      
      
      
    
% 
% 
% 
%     CTheta=0;
%     D1wl=(11/6*w(NnMidl-2)-3*w(NnMidl-3)+3/2*w(NnMidl-4)-1/3*w(NnMidl-5))/dNn;
%     D2wl=(35/12*w(NnMidl-2) -26/3*w(NnMidl-3)+ 19/2*w(NnMidl-4)-14/3*w(NnMidl-5)+11/12*w(NnMidl-6))/dNn.^2;
%     DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
%     
%     wll(NnMidl-1)=w(NnMidl-2)+  (D1wl+CTheta*D2wl*dNn+.5*DeltaD1w)*dNn;
%     wll(NnMidl)=wll(NnMidl-1)+  (D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
%     wll(NnMidh)=wll(NnMidl)+    (D1wl+CTheta*D2wl*dNn+2.5*DeltaD1w)*dNn;
%     wll(NnMidh+1)=wll(NnMidh)+  (D1wl+CTheta*D2wl*dNn+3.5*DeltaD1w)*dNn;
%     wll(NnMidh+2)=wll(NnMidh+1)+(D1wl+CTheta*D2wl*dNn+4.5*DeltaD1w)*dNn;
%     %wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
% %The last equation is not used since point w(NnMidh+3)is given but it is used in 
% %calculation of distance between two interpolation points. This distance is
% %used to determine the value of DeltaD1w.
% 
%     D1wh=(-11/6*w(NnMidh+3)+3*w(NnMidh+4)-3/2*w(NnMidh+5)+1/3*w(NnMidh+6))/dNn;
%     D2wh=(35/12*w(NnMidh+3) -26/3*w(NnMidh+4)+ 19/2*w(NnMidh+5)-14/3*w(NnMidh+6)+11/12*w(NnMidh+7))/dNn.^2;
%     DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wh*dNn-6.0*CTheta*D2wh*dNn.^2.0/1.0)/18.0/dNn;
% %     DeltaD1w=(w(NnMidh+2)-w(NnMidl-1)-4*D1wh*dNn-16.0*D2wh*dNn.^2.0/1.0)/6.0/dNn;
% 
% 
%     whh(NnMidh+2)=w(NnMidh+3)-(D1wh+CTheta*D2wh*dNn+.5*DeltaD1w)*dNn;
%     whh(NnMidh+1)=w(NnMidh+2)-(D1wh+CTheta*D2wh*dNn+1.5* DeltaD1w)*dNn;
%     whh(NnMidh)=whh(NnMidh+1)-(D1wh+CTheta*D2wh*dNn+2.5* DeltaD1w)*dNn;
%     whh(NnMidl)=whh(NnMidh)-(D1wh+CTheta*D2wh*dNn+3.5* DeltaD1w)*dNn;
%     whh(NnMidl-1)=whh(NnMidl)-(D1wh+CTheta*D2wh*dNn+4.5* DeltaD1w)*dNn;
% 
%     %whh(NnMidl-2)=whh(NnMidl-1)-(D1wh+CTheta*D2wh*dNn/1+5.5* DeltaD1w)*dNn;
% 
% 
%         w(NnMidh+2)=.5*whh(NnMidh+2)+.5*wll(NnMidh+2);
%         w(NnMidh+1)=.5*whh(NnMidh+1)+.5*wll(NnMidh+1);
%         w(NnMidh)=.5*whh(NnMidh)+.5*wll(NnMidh);
%         w(NnMidl)=.5*whh(NnMidl)+.5*wll(NnMidl);
%         w(NnMidl-1)=.5*whh(NnMidl-1)+.5*wll(NnMidl-1);


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

        
        w1(1:Nn-1)=w(1:Nn-1);
        w1(Nn)=w(Nn);
        w2(1:Nn-1)=w(2:Nn);
        w2(Nn)=wE;     
        
w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
     %Change 3:7/25/2020: I have improved zero correction in above.
        w(w<0)=0.0;
        for nn=1:Nn
            if(w(nn)<=0)
                wnStart=nn+1;
           end
        end
        
        %%1st order mean correction. We know the mean in original
        %%coordinates so I shift the density into original coordinates,
        %%apply the mean correction and then transform again into Lamperti
        %%coordinates. Algebra solves two equations in two unknowns.
        %%Equations are Sum_0^N{(Y_w(wnStart:Nn)-Y0)*W0}=1 and the second
        %%equation is Sum_0^N{Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0}=u0
        %%Two unknows are Y0 and W0. u0 is known mean.
        tt;
        u0=theta+(x0-theta)*exp(-kappa*(tt*dt)); %analytic mean of the density
        
        %If you are not using stochastic volatility, replace above with
        %true mean otherwise results would become garbage. 

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

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

toc

ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 
 disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
 TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average


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

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

%Uncomment for fourth order monte carlo

%   
%   YY(1:paths)=YY(1:paths) + ...
%       (YCoeff0(1,1,2,1).*YY(1:paths).^Fp(1,1,2,1)+YCoeff0(1,2,1,1).*YY(1:paths).^Fp(1,2,1,1)+ ...
%       YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1))*dtM + ...
%       (YCoeff0(1,1,3,1).*YY(1:paths).^Fp(1,1,3,1)+YCoeff0(1,2,2,1).*YY(1:paths).^Fp(1,2,2,1)+ ...
%       YCoeff0(2,1,2,1).*YY(1:paths).^Fp(2,1,2,1)+YCoeff0(1,3,1,1).*YY(1:paths).^Fp(1,3,1,1)+ ...
%       YCoeff0(2,2,1,1).*YY(1:paths).^Fp(2,2,1,1)+YCoeff0(3,1,1,1).*YY(1:paths).^Fp(3,1,1,1))*dtM^2 + ...
%       (YCoeff0(1,1,4,1).*YY(1:paths).^Fp(1,1,4,1)+YCoeff0(1,2,3,1).*YY(1:paths).^Fp(1,2,3,1)+ ...
%       YCoeff0(2,1,3,1).*YY(1:paths).^Fp(2,1,3,1)+YCoeff0(1,3,2,1).*YY(1:paths).^Fp(1,3,2,1)+ ...
%       YCoeff0(2,2,2,1).*YY(1:paths).^Fp(2,2,2,1)+YCoeff0(3,1,2,1).*YY(1:paths).^Fp(3,1,2,1)+ ...
%       YCoeff0(1,4,1,1).*YY(1:paths).^Fp(1,4,1,1)+YCoeff0(2,3,1,1).*YY(1:paths).^Fp(2,3,1,1)+ ...
%       YCoeff0(3,2,1,1).*YY(1:paths).^Fp(3,2,1,1)+YCoeff0(4,1,1,1).*YY(1:paths).^Fp(4,1,1,1))*dtM^3 + ...
%        (YCoeff0(1,1,5,1).*YY(1:paths).^Fp(1,1,5,1)+YCoeff0(1,2,4,1).*yy(1:paths).^Fp(1,2,4,1)+ ...
%     YCoeff0(2,1,4,1).*YY(1:paths).^Fp(2,1,4,1)+YCoeff0(2,1,1,1).*yy(1:paths).^Fp(2,1,1,1)+ ...
%      YCoeff0(2,2,3,1).*YY(1:paths).^Fp(2,2,3,1)+ ...
%      YCoeff0(3,1,3,1).*YY(1:paths).^Fp(3,1,3,1)+YCoeff0(1,4,2,1).*YY(1:paths).^Fp(1,4,2,1)+ ...
%      YCoeff0(2,3,2,1).*YY(1:paths).^Fp(2,3,2,1)+YCoeff0(3,2,2,1).*YY(1:paths).^Fp(3,2,2,1)+ ...
%      YCoeff0(4,1,2,1).*YY(1:paths).^Fp(4,1,2,1)+YCoeff0(1,5,1,1).*YY(1:paths).^Fp(1,5,1,1)+ ...
%      YCoeff0(2,4,1,1).*YY(1:paths).^Fp(2,4,1,1)+ ...
%      YCoeff0(3,3,1,1).*YY(1:paths).^Fp(3,3,1,1)+YCoeff0(4,2,1,1).*YY(1:paths).^Fp(4,2,1,1)+ ...
%      YCoeff0(5,1,1,1).*YY(1:paths).^Fp(5,1,1,1))*dtM^4+ ...
%       ((YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
%     (YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
%     YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5+ ...
%     (YCoeff0(1,1,3,2).*YY(1:paths).^Fp(1,1,3,2)+YCoeff0(1,2,2,2).*YY(1:paths).^Fp(1,2,2,2)+ ...
%     YCoeff0(2,1,2,2).*YY(1:paths).^Fp(2,1,2,2)+YCoeff0(1,3,1,2).*YY(1:paths).^Fp(1,3,1,2)+ ...
%     YCoeff0(2,2,1,2).*YY(1:paths).^Fp(2,2,1,2)+YCoeff0(3,1,1,2).*YY(1:paths).^Fp(3,1,1,2)).*dtM^2.5+ ...
%     (YCoeff0(1,1,4,2).*YY(1:paths).^Fp(1,1,4,2)+YCoeff0(1,2,3,2).*YY(1:paths).^Fp(1,2,3,2)+ ...
%     YCoeff0(2,1,3,2).*YY(1:paths).^Fp(2,1,3,2)+YCoeff0(1,3,2,2).*YY(1:paths).^Fp(1,3,2,2)+ ...
%     YCoeff0(2,2,2,2).*YY(1:paths).^Fp(2,2,2,2)+ YCoeff0(3,1,2,2).*YY(1:paths).^Fp(3,1,2,2)+ ...
%     YCoeff0(1,4,1,2).*YY(1:paths).^Fp(1,4,1,2)+YCoeff0(2,3,1,2).*YY(1:paths).^Fp(2,3,1,2)+ ...
%     YCoeff0(3,2,1,2).*YY(1:paths).^Fp(3,2,1,2)+YCoeff0(4,1,1,2).*YY(1:paths).^Fp(4,1,1,2)).*dtM^3.5) .*HermiteP1(2,1:paths) + ...
%     ((YCoeff0(1,1,1,3).*YY(1:paths).^Fp(1,1,1,3) *dtM) + ...
%     (YCoeff0(1,1,2,3).*YY(1:paths).^Fp(1,1,2,3)+YCoeff0(1,2,1,3).*YY(1:paths).^Fp(1,2,1,3)+ ...
%     YCoeff0(2,1,1,3).*YY(1:paths).^Fp(2,1,1,3)).*dtM^2+ ...
%     (YCoeff0(1,1,3,3).*YY(1:paths).^Fp(1,1,3,3)+YCoeff0(1,2,2,3).*YY(1:paths).^Fp(1,2,2,3)+ ...
%     YCoeff0(2,1,2,3).*YY(1:paths).^Fp(2,1,2,3) + YCoeff0(1,3,1,3).*YY(1:paths).^Fp(1,3,1,3)+ ...
%     YCoeff0(2,2,1,3).*YY(1:paths).^Fp(2,2,1,3)+YCoeff0(3,1,1,3).*YY(1:paths).^Fp(3,1,1,3)).*dtM^3).*HermiteP1(3,1:paths) + ...
%     ((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )+ ...
%     (YCoeff0(1,1,2,4).*YY(1:paths).^Fp(1,1,2,4)+YCoeff0(1,2,1,4).*YY(1:paths).^Fp(1,2,1,4)+ ...
%     YCoeff0(2,1,1,4).*YY(1:paths).^Fp(2,1,1,4))*dtM^2.5).*HermiteP1(4,1:paths) + ...
%     (YCoeff0(1,1,1,5).*YY(1:paths).^Fp(1,1,1,5)*dtM^2.0).*HermiteP1(5,1:paths);
%     
 
     YYMean(tt)=YYMean(tt)+sum(YY(1:paths))/paths;
 end
 YY(YY<0)=0;
 disp('Original process average from monte carlo');
 MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
 disp('Original process average from our simulation');
 ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 
 disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
 TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
 
 
 MaxCutOff=30;
 NoOfBins=round(300*gamma^2*4*sigma0);%Decrease the number of bins if the graph is too 
 [YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
 plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 %plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');
 
 title(sprintf('Ito-Hermite Method VS Monte Carlo Density, x0 = %.2f,theta=%.2f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f', x0,theta,kappa,gamma,sigma0,T));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
 %title('text',sprintf('x0 = %f', x0),sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T))
 %xlabel(a)
 %title('Graphical Comparison of Ito-Hermite Method and Monte Carlo Method Density')
 legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
 
 str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
 

end
Here is the result of the function run.

Original process average from monte carlo
MCMean =
   0.099651168377852
Original process average from our simulation
ItoHermiteMean =
   0.098879674266063
true Mean only applicble to standard SV mean reverting type models otherwise disregard
TrueMean =
   0.100000000000000
IndexMax =
   381

And here is the output graph

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

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

July 31st, 2020, 2:35 am

Sorry friends but the part of work we earlier did for FPE related to change of derivative from SDE variable to normal random variable and related removal of probability formulas are only approximate. Since we were not getting pinpoint precision, I decided to numerically verify those formulas and realized they have very small lax. Later, I will be posting my worked out matlab programs for comparison of numerical and analytical first derivative calculations based on change of derivative. I have some ideas about correcting this (slight) error and I will be following them in my program as well.  
 
User avatar
Amin
Topic Author
Posts: 2536
Joined: July 14th, 2002, 3:00 am

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

July 31st, 2020, 9:48 am

I checked the replication of first derivative by analytical method using a comparison between numerical first derivative and analytical first derivative. Here by analytical derivative, I mean the derivative obtained by differentiating the equation for transformed probability distribution that we get multiplying the original (normal) probability distribution with change of probability derivative. There are some light discrepancies but they seem to be within numerical errors. The analytical first order derivative seemed to be quite reasonably close to original derivative and I specifically checked that effect of 2nd derivative was close to what it should be. I had roughly checked the derivatives when I had written the original equations several months ago. 
Since first derivative seemed very fine, I could not explain the anomaly why 2nd derivative(that arises in expansion of first analytical derivative of FPE) is multiplied with .707. When I used .707 for the derivative, there were many instabilities and equation would blow again and again when I kept the coefficient of 2nd derivative at one. But now I was able to get even more nice and more stable solution by keeping the coefficient of 2nd derivative at one and altering some other coefficients(C33). Now I get a far better match to shape of true distribution with a better mean in most cases. I will be posting a new program later tonight or tomorrow with this new change(coefficient of 2nd derivative set equal to one).   
Though the solution was reasonable earlier, with this new change the accuracy is greatly improved and a lot of cases with extreme parameters which were giving bad results and where we could not solve the equation earlier, have greatly improved with much better accuracy. For general easy parameters, not only the shape is perfect, the mean is also in a very narrow range to the true mean.
 
User avatar
Amin
Topic Author
Posts: 2536
Joined: July 14th, 2002, 3:00 am

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

July 31st, 2020, 4:49 pm

I will post the program later but here are the constants I am using

  C11(wnStart:Nn)=1;
  C22=-1/pi/8;
  C33=1.15/pi/16;

I have not checked carefully and I am not sure if it is still needed or not but I am still weighting 2nd derivative with an exponential for smaller values where it does not matter that much and just in order for friends to get my output I am writing the weighting exponential equation here.

d2wdZ2(1:18)=d2wdZ2(1:18).*exp(-16*kappa^2*(1-gamma).*dNn*(18-(1:18)));

I am getting great accuracy and mean in every easy and medium difficulty case is within 0-15 bp of true mean. Though gamma closer to .5 (equal to or less than .65)  near zero can give worse mean values but even those cases have far more accuracy now.

I will be posting a complete program in a few hours or possibly tomorrow.
 
User avatar
Amin
Topic Author
Posts: 2536
Joined: July 14th, 2002, 3:00 am

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

August 1st, 2020, 1:31 pm

I could not work today because my persecution has dramatically increased for past one day. Please protest against my continued persecution and try to help if you can. Please rad here about the latest persecution activities. https://forum.wilmott.com/viewtopic.php?f=15&t=94796&p=859574#p859574
 
User avatar
Amin
Topic Author
Posts: 2536
Joined: July 14th, 2002, 3:00 am

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

August 1st, 2020, 4:51 pm

In the graphs below I have simulated fokker-planck equation of density of following SDE with various values of parameters.

[$]dx(t)=\kappa (\theta-x(t)) + \sigma x(t)^{\gamma} dz(t) [$]
with [$]x(0)=x_0[$]
You have to interpret the parameters given in title of each graph with reference to above SDE. I have solved the fokker-Planck equation related to above SDE in Lamperti/Bessel form.

I am posting graphs from the new program. I have also included effect from third derivative term. Fit to the true density is perfectly excellent. Mean of the density also very closely matches the true mean of the SDE densities. All of these densities are mean reverting SDE densities that I have been working with for a long time. 
Parameters are posted on graph title. FPE new method density mean and true mean is also posted on graph title. These graphs are totally beautiful.
very close to zero, at high volatilities, gamma parameter values around .5 still remain problematic but might work better with decreased step size.
I have made these graphs from the one program without any change and I will be posting the program in a few minutes in next post.

Image

Image

Image

Image

Image

Image

Image

Image

Image

Image

Image

Image

Image

Image

Image

 
Last edited by Amin on August 4th, 2020, 6:46 pm, edited 3 times in total.
 
User avatar
Amin
Topic Author
Posts: 2536
Joined: July 14th, 2002, 3:00 am

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

August 1st, 2020, 5:00 pm

Here is the program I used for graphs in previous post.


function [TrueMean,ItoHermiteMean,y_w,fy_w,wnStart] = FPERevisited044Wilmott16_2ndD()

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

%I have not directly simulated the SDE but simulated the transformed 
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.
%In this version, the mean correction has been disabled. If you are
%simulating standard mean-reverting SV type SDEs, please enable the mean 
%correction by uncommenting the appropriate line in the body of code below.

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


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

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

x0=.04;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.75;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=4.0;%.950;   %mean reversion parameter.
theta=.04;%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;
                
w(1:Nn)=x0^(1-gamma)/(1-gamma);

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

Z
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);

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

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

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

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

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

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

for k = 0 : (OrderA)
    for m = 0:k
        l4 = k - m + 1;
        for n = 0 : m
            l3 = m - n + 1;
            for j = 0:n
                l2 = n - j + 1;
                l1 = j + 1;
                %Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
                Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                
                YCoeff0(l1,l2,l3,l4) =YCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
                YqCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
            end
        end
    end
end

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

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

for tt=1:Tt
    
     yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
    
    [wMu0dt,dwMu0dtdw,d2wMu0dtdw2,d3wMu0dtdw3,c1,c2,c3,c4] = CalculateDriftAndVolOrder4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);

    wPrev=w;
    dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn) ;% ...
    dw2(wnStart:Nn)=dw(wnStart:Nn).^2;
    w(isnan(w)==1)=0;
    wMu0dt(isnan(wMu0dt)==1)=0;  
    wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
   
    yyMedian=.5*yy(NnMidh)+.5*yy(NnMidl);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

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

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

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



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


dZdw(wnStart:Nn)=1.0./dwdZ(wnStart:Nn);

if(tt==1)
w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+ ...
          sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
          sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
          sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)));
end
if(tt>1)
    
  
  C22=1/pi/8;
  C33=.9/pi/16;%  
  C44=1/pi/8;

  
   w(wnStart:Nn)=w(wnStart:Nn)+wMu0dt(wnStart:Nn)+(-(Z(wnStart:Nn).*dwdZ(wnStart:Nn))+d2wdZ2(wnStart:Nn))+ ...
             sign((Z(wnStart:Nn).*dwdZ(wnStart:Nn))-d2wdZ2(wnStart:Nn)+dw(wnStart:Nn)).* ...
             sqrt(abs(sign((Z(wnStart:Nn).*dwdZ(wnStart:Nn))-d2wdZ2(wnStart:Nn)).* ... 
             ((Z(wnStart:Nn).*dwdZ(wnStart:Nn))-d2wdZ2(wnStart:Nn)).^2+ ...
             sign(+dw(wnStart:Nn)).*dw2(wnStart:Nn)))+ ...
             -C33.*dwMu0dtdw(wnStart:Nn).*dwdZ(wnStart:Nn).^2 + ...
             +C22.*-3*Z(wnStart:Nn).^1.*dZdw(wnStart:Nn).*d2wdZ2(wnStart:Nn)*sigma0^2.*dt + ...
             +C44.*(-d3wdZ3(wnStart:Nn).*dZdw(wnStart:Nn)+3*d2wdZ2(wnStart:Nn).^2.*dZdw(wnStart:Nn).^2)*sigma0^2.*dt;
        
        
        %Change 2:7/25/2020:%I have added a last term from quadratic variation between drift
        %and w terms. I am not very sure if we need it or not.

  
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Below are instability correction equations. There is instability in the 
% SDE representation equation roughly between -.5 SD to +.5 SD in normal terms.
%Instability correction Algorithm takes two ends that are here given as
%lower end=w(NnMidl-2) and upper end=w(NnMidh+3).
%In the block below the interpolation starts from lower end.
%1st equation is
%wll(NnMidl-1)=w(NnMidl-2)+(D1wl+CTheta*D2wl*dNn/1+.5*DeltaD1w)*dNn;
%Second equation is: wll(NnMidl)=wll(NnMidl-1)+(D1wl+CTheta*D2wl*dNn+1.5*DeltaD1w)*dNn;
%And so on.
%There is a sixth equation that is not written below(since it takes us to 
%the upper interpolation point which is already known) which is 
%wll(NnMidh+3)=wll(NnMidh+2)+(D1wl+CTheta*D2wl*dNn+5.5*DeltaD1w)*dNn;
%Sum of these six equations add up to the distance between upper and lower
%interpolation points.
%In the sum there are six instances of first derivative which is already 
%calculated at the boundary, six instances of second derivative which is also
%already calculated at lower boundary and 18.0 instances of DeltaD1w which is 
%unknown. This equation calculates the value of one unit of DeltaD1w addition
%DeltaD1w=(w(NnMidh+3)-w(NnMidl-2)-6*D1wl*dNn-6.0*CTheta*D2wl*dNn.^2.0/1.0)/18.0/dNn;
% I hope to find a better analytical method to tackle this issue. Friends
% are welcome to play with slightly changed related alternatives to the 
%above interpolation.
%Change 7/27/2020
%I have switched to Lagrange interpolation by taking four points on both 
%sides of instability and using 8th order interpolation of five points in
%mid between the eight points on sides. This works far better and you can
%now use nuemrical 2nd derivative without any modification for a large 
%range of parameters. It truly improves the fit to true density when it 
%works. 

      
      Zi=[Z(NnMidl-6),Z(NnMidl-5),Z(NnMidl-4),Z(NnMidl-3),Z(NnMidl-2),Z(NnMidh+3),Z(NnMidh+4),Z(NnMidh+5),Z(NnMidh+6),Z(NnMidh+7)];
      wi=[w(NnMidl-6),w(NnMidl-5),w(NnMidl-4),w(NnMidl-3),w(NnMidl-2),w(NnMidh+3),w(NnMidh+4),w(NnMidh+5),w(NnMidh+6),w(NnMidh+7)];
      Zq=[Z(NnMidl-1),Z(NnMidl),Z(NnMidh),Z(NnMidh+1),Z(NnMidh+2)];
      wq=spline(Zi,wi,Zq);
      w(NnMidl-1:NnMidh+2)=wq(1:5);
      
      

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

        
        w1(1:Nn-1)=w(1:Nn-1);
        w1(Nn)=w(Nn);
        w2(1:Nn-1)=w(2:Nn);
        w2(Nn)=wE;     
        
    w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
     %Change 3:7/25/2020: I have improved zero correction in above.
        w(w<0)=0.0;
        for nn=1:Nn
            if(w(nn)<=0)
                wnStart=nn+1;
           end
        end
        
        %%1st order mean correction. We know the mean in original
        %%coordinates so I shift the density into original coordinates,
        %%apply the mean correction and then transform again into Lamperti
        %%coordinates. Algebra solves two equations in two unknowns.
        %%Equations are Sum_0^N{(Y_w(wnStart:Nn)-Y0)*W0}=1 and the second
        %%equation is Sum_0^N{Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0}=u0
        %%Two unknows are Y0 and W0. u0 is known mean.
        tt;
        u0=theta+(x0-theta)*exp(-kappa*(tt*dt)); %analytic mean of the density
        
        %If you are not using stochastic volatility, replace above with
        %true mean otherwise results would become garbage. 

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

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

toc

ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 
 disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
 TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average


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

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

%Uncomment for fourth order monte carlo

%   
%   YY(1:paths)=YY(1:paths) + ...
%       (YCoeff0(1,1,2,1).*YY(1:paths).^Fp(1,1,2,1)+YCoeff0(1,2,1,1).*YY(1:paths).^Fp(1,2,1,1)+ ...
%       YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1))*dtM + ...
%       (YCoeff0(1,1,3,1).*YY(1:paths).^Fp(1,1,3,1)+YCoeff0(1,2,2,1).*YY(1:paths).^Fp(1,2,2,1)+ ...
%       YCoeff0(2,1,2,1).*YY(1:paths).^Fp(2,1,2,1)+YCoeff0(1,3,1,1).*YY(1:paths).^Fp(1,3,1,1)+ ...
%       YCoeff0(2,2,1,1).*YY(1:paths).^Fp(2,2,1,1)+YCoeff0(3,1,1,1).*YY(1:paths).^Fp(3,1,1,1))*dtM^2 + ...
%       (YCoeff0(1,1,4,1).*YY(1:paths).^Fp(1,1,4,1)+YCoeff0(1,2,3,1).*YY(1:paths).^Fp(1,2,3,1)+ ...
%       YCoeff0(2,1,3,1).*YY(1:paths).^Fp(2,1,3,1)+YCoeff0(1,3,2,1).*YY(1:paths).^Fp(1,3,2,1)+ ...
%       YCoeff0(2,2,2,1).*YY(1:paths).^Fp(2,2,2,1)+YCoeff0(3,1,2,1).*YY(1:paths).^Fp(3,1,2,1)+ ...
%       YCoeff0(1,4,1,1).*YY(1:paths).^Fp(1,4,1,1)+YCoeff0(2,3,1,1).*YY(1:paths).^Fp(2,3,1,1)+ ...
%       YCoeff0(3,2,1,1).*YY(1:paths).^Fp(3,2,1,1)+YCoeff0(4,1,1,1).*YY(1:paths).^Fp(4,1,1,1))*dtM^3 + ...
%        (YCoeff0(1,1,5,1).*YY(1:paths).^Fp(1,1,5,1)+YCoeff0(1,2,4,1).*yy(1:paths).^Fp(1,2,4,1)+ ...
%     YCoeff0(2,1,4,1).*YY(1:paths).^Fp(2,1,4,1)+YCoeff0(2,1,1,1).*yy(1:paths).^Fp(2,1,1,1)+ ...
%      YCoeff0(2,2,3,1).*YY(1:paths).^Fp(2,2,3,1)+ ...
%      YCoeff0(3,1,3,1).*YY(1:paths).^Fp(3,1,3,1)+YCoeff0(1,4,2,1).*YY(1:paths).^Fp(1,4,2,1)+ ...
%      YCoeff0(2,3,2,1).*YY(1:paths).^Fp(2,3,2,1)+YCoeff0(3,2,2,1).*YY(1:paths).^Fp(3,2,2,1)+ ...
%      YCoeff0(4,1,2,1).*YY(1:paths).^Fp(4,1,2,1)+YCoeff0(1,5,1,1).*YY(1:paths).^Fp(1,5,1,1)+ ...
%      YCoeff0(2,4,1,1).*YY(1:paths).^Fp(2,4,1,1)+ ...
%      YCoeff0(3,3,1,1).*YY(1:paths).^Fp(3,3,1,1)+YCoeff0(4,2,1,1).*YY(1:paths).^Fp(4,2,1,1)+ ...
%      YCoeff0(5,1,1,1).*YY(1:paths).^Fp(5,1,1,1))*dtM^4+ ...
%       ((YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
%     (YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
%     YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5+ ...
%     (YCoeff0(1,1,3,2).*YY(1:paths).^Fp(1,1,3,2)+YCoeff0(1,2,2,2).*YY(1:paths).^Fp(1,2,2,2)+ ...
%     YCoeff0(2,1,2,2).*YY(1:paths).^Fp(2,1,2,2)+YCoeff0(1,3,1,2).*YY(1:paths).^Fp(1,3,1,2)+ ...
%     YCoeff0(2,2,1,2).*YY(1:paths).^Fp(2,2,1,2)+YCoeff0(3,1,1,2).*YY(1:paths).^Fp(3,1,1,2)).*dtM^2.5+ ...
%     (YCoeff0(1,1,4,2).*YY(1:paths).^Fp(1,1,4,2)+YCoeff0(1,2,3,2).*YY(1:paths).^Fp(1,2,3,2)+ ...
%     YCoeff0(2,1,3,2).*YY(1:paths).^Fp(2,1,3,2)+YCoeff0(1,3,2,2).*YY(1:paths).^Fp(1,3,2,2)+ ...
%     YCoeff0(2,2,2,2).*YY(1:paths).^Fp(2,2,2,2)+ YCoeff0(3,1,2,2).*YY(1:paths).^Fp(3,1,2,2)+ ...
%     YCoeff0(1,4,1,2).*YY(1:paths).^Fp(1,4,1,2)+YCoeff0(2,3,1,2).*YY(1:paths).^Fp(2,3,1,2)+ ...
%     YCoeff0(3,2,1,2).*YY(1:paths).^Fp(3,2,1,2)+YCoeff0(4,1,1,2).*YY(1:paths).^Fp(4,1,1,2)).*dtM^3.5) .*HermiteP1(2,1:paths) + ...
%     ((YCoeff0(1,1,1,3).*YY(1:paths).^Fp(1,1,1,3) *dtM) + ...
%     (YCoeff0(1,1,2,3).*YY(1:paths).^Fp(1,1,2,3)+YCoeff0(1,2,1,3).*YY(1:paths).^Fp(1,2,1,3)+ ...
%     YCoeff0(2,1,1,3).*YY(1:paths).^Fp(2,1,1,3)).*dtM^2+ ...
%     (YCoeff0(1,1,3,3).*YY(1:paths).^Fp(1,1,3,3)+YCoeff0(1,2,2,3).*YY(1:paths).^Fp(1,2,2,3)+ ...
%     YCoeff0(2,1,2,3).*YY(1:paths).^Fp(2,1,2,3) + YCoeff0(1,3,1,3).*YY(1:paths).^Fp(1,3,1,3)+ ...
%     YCoeff0(2,2,1,3).*YY(1:paths).^Fp(2,2,1,3)+YCoeff0(3,1,1,3).*YY(1:paths).^Fp(3,1,1,3)).*dtM^3).*HermiteP1(3,1:paths) + ...
%     ((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )+ ...
%     (YCoeff0(1,1,2,4).*YY(1:paths).^Fp(1,1,2,4)+YCoeff0(1,2,1,4).*YY(1:paths).^Fp(1,2,1,4)+ ...
%     YCoeff0(2,1,1,4).*YY(1:paths).^Fp(2,1,1,4))*dtM^2.5).*HermiteP1(4,1:paths) + ...
%     (YCoeff0(1,1,1,5).*YY(1:paths).^Fp(1,1,1,5)*dtM^2.0).*HermiteP1(5,1:paths);
%     
 
     YYMean(tt)=YYMean(tt)+sum(YY(1:paths))/paths;
 end
 YY(YY<0)=0;
 disp('Original process average from monte carlo');
 MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
 disp('Original process average from our simulation');
 ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 
 disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
 TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
 
 
 MaxCutOff=30;
 NoOfBins=round(300*gamma^2*4*sigma0/sqrt(yyMedian)/(1+kappa));%Decrease the number of bins if the graph is too 
 [YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
 plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 %plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');
 
 title(sprintf('x0 = %.4f,theta=%.2f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,Mean=%.5f,TrueMean=%.5f', x0,theta,kappa,gamma,sigma0,T,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
 %title('text',sprintf('x0 = %f', x0),sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T))
 %xlabel(a)
 %title('Graphical Comparison of Ito-Hermite Method and Monte Carlo Method Density')
 legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
 
 str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
 

end
This program needs the following new function to calculate first three derivatives on equi-spaced grid.
function [dwdZ,d2wdZ2,d3wdZ3] = First3Derivatives2ndOrderEqSpaced(wnStart,Nn,dNn,w,Z)

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

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

 [wS] = InterpolateOrderN6(5,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5));
 [wE] = InterpolateOrderN6(5,Z(Nn)+dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));
 
 
 Zi=[Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),Z(wnStart+8),Z(wnStart+9)];
      wi=[w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7),w(wnStart+8),w(wnStart+9)];
%      Zi=[Z(Nn-6),Z(Nn-5),Z(Nn-4),Z(Nn-3),Z(Nn-2),Z(Nn-1),Z(Nn)];
%      wi=[w(Nn-6),w(Nn-5),w(Nn-4),w(Nn-3),w(Nn-2),w(Nn-1),w(Nn)];

      Zi=[Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6)];
      wi=[w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6)];
      Zq=Z(wnStart)-dNn;
      wS=spline(Zi,wi,Zq);
 
      Zq=Z(wnStart)-2*dNn;
      wS2=spline(Zi,wi,Zq);
 
 
 
%      Zi=[Z(Nn-9),Z(Nn-8),Z(Nn-7),Z(Nn-6),Z(Nn-5),Z(Nn-4),Z(Nn-3),Z(Nn-2),Z(Nn-1),Z(Nn)];
%      wi=[w(Nn-9),w(Nn-8),w(Nn-7),w(Nn-6),w(Nn-5),w(Nn-4),w(Nn-3),w(Nn-2),w(Nn-1),w(Nn)];
      Zi=[Z(Nn-6),Z(Nn-5),Z(Nn-4),Z(Nn-3),Z(Nn-2),Z(Nn-1),Z(Nn)];
      wi=[w(Nn-6),w(Nn-5),w(Nn-4),w(Nn-3),w(Nn-2),w(Nn-1),w(Nn)];
      Zq=Z(Nn)+dNn;
      wE=spline(Zi,wi,Zq);
     
 
      Zq=Z(Nn)+2*dNn;
      wE2=spline(Zi,wi,Zq);
 
 
 
 ZS=Z(wnStart)-dNn;
 ZE=Z(Nn)+dNn;

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

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

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

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

August 2nd, 2020, 8:45 am

I want to say that the program I posted above is not final. I will be posting another program with slight fine-tuning of parameters. 
So I expect a few minor changes in a day or two.
Later, I will be posting another program with a grid of N=100 just so that friends can choose any grid size they wish. Please note that the present grid has a size of N=50.
Edit: I edited the above post and removed reference to sign of 3rd derivative.
 
User avatar
Amin
Topic Author
Posts: 2536
Joined: July 14th, 2002, 3:00 am

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

August 3rd, 2020, 5:01 pm

I was trying to post on some linkedin groups asking friends to download the program in previous posts but my posts do not appear on some groups even after I see the notification that my post was successful. It seems that mind control and allied crooks have given money to their horses who want to present my research as their own. Why are they stopping my research from reaching other people?
 
User avatar
Amin
Topic Author
Posts: 2536
Joined: July 14th, 2002, 3:00 am

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

August 3rd, 2020, 6:29 pm

Lot of friends in several countries can recall how crooks in US defence and their allied apparatus were telling my linkedin connections when I did my monte carlo simulation research and approached my connections about my research that I was a very bad guy and "horses" of US defence crooks would present the same (my) research to them. Well, these crooks never grow up and they want to do the same when I have presented solution to fokker-planck equation.
I will try to post latex equations for the solution of Fokker-planck equation later today. But in the meantime, I want to request friends to protest against my persecution. Here I want to share some threads with friends especially my application to human rights watch for help.
My Letter To HRW About My Human Rights Abuse 
Muslims, mind control, Muslim-Jewish relationship and American defense
My Open Letter to Prime Minister of Pakistan, Imran Khan, to Intervene In Order to Protect the Interests of Pakistan
 
User avatar
Amin
Topic Author
Posts: 2536
Joined: July 14th, 2002, 3:00 am

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

August 5th, 2020, 1:04 pm

Post in progress.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

We can write the above equation after noting that some terms are getting constant multiplier coefficients not immediately known. So we re-write equation (17) as
 
[$].5 \Big[X(t_1)-X(t_0) - \int_{t_0}^{t_1} \mu(X) ds - \int_{t_0}^{t_1} 2 \sigma(X) \frac{\partial \sigma(X)}{\partial X}  ds [$]
[$]  + ( Z (\frac{\partial X}{\partial Z}) -  \frac{\partial^2 X}{\partial Z^2})\Big]^2 [$] 
[$]=  +{\Big[ ( Z (\frac{\partial X}{\partial Z}) -  \frac{\partial^2 X}{\partial Z^2})\Big]}^2 [$]
[$]+ .5 {\sigma(X)}^2 (Z^2 +3C_2 Z  {(\frac{\partial Z}{\partial X})} \frac{\partial^2 X}{\partial Z^2} +C_4 {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3}) dt[$] 
[$] - C_3\frac{\partial \mu(X)}{\partial X}  {(\frac{\partial X}{\partial Z})}^2 dt[$]
[$]+C_5({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) {(\frac{\partial X}{\partial Z})}^2 dt[$]  Equation(18)

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

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

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

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

Please pardon any error with signs. I quickly wrote the post but I will carefully check it again and fix any errors in a day or two.
ABOUT WILMOTT

PW by JB

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


Twitter LinkedIn Instagram

JOBS BOARD

JOBS BOARD

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


GZIP: On