SERVING THE QUANTITATIVE FINANCE COMMUNITY

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

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

February 11th, 2020, 3:18 pm

I was thinking about all the possible explanations behind second and third order corrections to the Ito-hermite SDE evolution algorithm. Ok I have some ideas but please take them with a pinch of salt. I believe that realized volatility along time for every point on SD-fraction grid is proportional to path integral of the drift along the SD-lines or constant Z-lines on which the density of w is based. Instantaneous volatility is almost constant and it is path integral of drift that is affecting the variance evolution along the SD-fraction or constant Z-lines of w which is the SDE variable in our case in bessel coordinates. In one dimensional case, w is evolving along one axis as a straight line which expands as a function of time. when w takes just one derivative with respect to Z, we are doing fine and we do not have to make any correction. However whenever w is convex or concave with respect to Z, and has substantial second or third derivatives, we have to make some correction to path integral of drift of w along SD-fractions to get appropriate volatility. We have to remove the effect of higher derivatives with respect to Z in order to find the path integrals of drift along a straight expanding line? 
I will be doing more experiments in coming days on these lines and share with friends.
 
User avatar
Amin
Topic Author
Posts: 2320
Joined: July 14th, 2002, 3:00 am

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

February 12th, 2020, 9:43 am

Some comments before I share my program. In order to update the density we have to add existing density and drift  independently along the orthogonal coordinates in a squared fashion and then take square root. It would probably also take care of instability as well.
But in current version, I just played with legacy version and updated it by subtracting squared sum along orthogonal coordinates. It will be followed by a final version in a few days since we have almost taken our work to conclusion. This version is also pretty good but it is not final at all as I just hastily prepared it. Ok here is the new code.
function [] = ItoHermiteMethodWilmott08MeanCorrected02A()
%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 coordinates.
%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;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.

Tt=16*8;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year; 
OrderA=2;  %Keep at two.
OrderM=2;  %Keep at two.
%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;

x0=.250;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.850;   % volatility power.                                                                                                                                                                                                                                                                           
kappa=1.4850;   %mean reversion parameter.
theta=.1;   %mean reversion target

%you can specify any general mu1 and mu2 and beta1 and beta2.
%The present program is not perfect but quite reasonably good.
%if there is only one term in drift with beta=0, there will be minor
%errors. This case is yet not covered.

mu1=+1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
sigma0=.850;%Volatility value
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

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



 dt2=dt^2/2;
 dt3=dt^3/6.0;
 dtsqrt=sqrt(dt);
 dtonehalf=dt.^1.5;
 
wnStart=1;%


tic

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

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

    wMu0dt(wnStart:Nn)=GG(wnStart:Nn)*dt + dGG(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...
        .5*dGG2(wnStart:Nn).*QQ(wnStart:Nn).*dt2 + ...
        ddGGDD(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
        .5*ddGG2QQ(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
        .5*d2dGGDD(wnStart:Nn).*QQ(wnStart:Nn).*dt3 + ...
        .25*d2dGG2QQ(wnStart:Nn).*QQ(wnStart:Nn).*dt3;
    %    c0(wnStart:Nn)=-1*d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0;
    c1(wnStart:Nn)=sigma0.*dtsqrt+dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*dtonehalf;
    c2(wnStart:Nn)=d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0;
    dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn);% + ...
    dw2(wnStart:Nn)=c1(wnStart:Nn).^2.0.*Z(wnStart:Nn).^2;
    w(isnan(w)==1)=0;
    wMu0dt(isnan(wMu0dt)==1)=0;  
    wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%In this code block, we use second hermite representation of the density
%which is used to solve for mean reversion.
%It is the value of Sigma1 below whcih is used later in the main SDE one
%time step update formula.

    if(tt>1)
        [wMean,ZMean,Sigma00,Sigma11,Sigma22,Sigma1,Sigma2,Sigma3,Sigma0,w1New,w2New] = CalculateHermiteRepresentationOfDensityOrderThree(w,Z,ZProb,wnStart,Nn,dNn,NnMidh,NnMidl);
        [wMu0dtMean,ZMeanM,Sigma00M,Sigma11M,Sigma22M,Sigma1M,Sigma2M,Sigma3M,Sigma0M,wMu0dt1New,wMu0dt2New] = CalculateHermiteRepresentationOfDensityOrderThree(wMu0dt,Z,ZProb,wnStart,Nn,dNn,NnMidh,NnMidl);

%        Sigma00
%        Sigma11
%        Sigma22
%        Sigma3
%        tt
     % str=input('Look at vols');
    else
        Sigma1(wnStart:Nn)=0.0;
        Sigma2(wnStart:Nn)=0.0;
        Sigma3(wnStart:Nn)=0.0;
        Sigma1M(wnStart:Nn)=0.0;
        Sigma2M(wnStart:Nn)=0.0;
        Sigma11=0.0;
        Sigma11M=0.0;
        Sigma22=0.0;
    end
   
   
%This is the main SDE variable one step update equation. 


 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)))- ...
            sqrt(1)*(dt).*sqrt(Sigma11.^2+Sigma11M.^2).*(Z(wnStart:Nn).^2-1)+ ...
            (1/3)*(dt).*sqrt(Sigma2(wnStart:Nn).^2+Sigma2M(wnStart:Nn).^2).*(Z(wnStart:Nn).^3-3*Z(wnStart:Nn));
      
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%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.
    CTheta=1.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+2)+3*w(NnMidh+3)-3/2*w(NnMidh+4)+1/3*w(NnMidh+5))/dNn;
    D2wh=(35/12*w(NnMidh+2) -26/3*w(NnMidh+3)+ 19/2*w(NnMidh+4)-14/3*w(NnMidh+5)+11/12*w(NnMidh+6))/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);



             
        w1(1:Nn-1)=w(1:Nn-1);
        w2(1:Nn-1)=w(2:Nn);
        w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
     
        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
str=input('Analytic Values calculated; Press a key to start monte carlo');


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.
%Pre-compute the time and power exponent values in small multi-dimensional arrays
for k = 0 : (OrderA+1)
    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));

            end
        end
    end
end
%Below call the program that calculates the Ito-Taylor coeffs.
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);





%Yq=ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);
%YqCoeff=Yq/(1-gamma);%These coordinates are divided by (1-gamma)
 %Calculate the  density with monte carlo below.
 %rng(85109634, 'twister')
 rng(29079137, 'twister')
 paths=250000;
 YY(1:paths)=x0;  %Original process monte carlo.
% YY2dt(1:paths)=0.0;
 %YY4dt(1:paths)=0.0;
 X(1:paths)=x0^(1-gamma)/(1-gamma);
 Random1(1:paths)=0;
 YYMean(1:TtM)=0;
 for tt=1:TtM
 %    t1=tt*dtM;
 %    t0=(tt-1)*dtM;
     Random1=randn(size(Random1));
     HermiteP1(1,1:paths)=1;
     HermiteP1(2,1:paths)=Random1(1:paths);
     HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
     HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
     HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;
     YY0=YY;
     YY1=YY;
  %   YYPrev=YY;
  %   XX0=X;
     
     for k = 0 : OrderM
         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;
                     if(l4<=5)   
                         
  %                    X(1:paths)=X(1:paths) + YqCoeff(l1,l2,l3,l4) * ...
  %                   mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
  %                   ((1-gamma)*XX0(1:paths)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* ...
  %                   HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);    
                 
                     YY(1:paths)=YY(1:paths) + YCoeff(l1,l2,l3,l4) * ...
                     mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
                     YY0(1:paths).^Fp(l1,l2,l3,l4) .* HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);
%                     YY(YY>0)=YY1(YY>0);
%                     YY(YY<0)=0;
                 %Above is original coordinates monte carlo equation.
                     end
                 end
             end
         end
     end
     %YY(YY<0)=0;
     YYMean(tt)=YYMean(tt)+sum(YY(1:paths))/paths;
 end
 YY(YY<0)=0;
 disp('Original process average from monte carlo');
 sum(YY(:))/paths %origianl coordinates monte carlo average.
 disp('Original process average from our simulation');
 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');
 theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
 
 
 MaxCutOff=30;
 NoOfBins=500*gamma^2/.25*sigma0;
 %[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
 [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');
 title('Graphical Comparison of Ito-Hermite Method and Monte Carlo Method Density')
 legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
 
 str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
 
  plot((1:Tt),YYMean(1:Tt),'g',(1:Tt),yyMean(1:Tt),'r',(1:Tt),yyMean0(1:Tt),'b');
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [wMean,ZMean,Sigma00,Sigma11,Sigma22,Sigma1,Sigma2,Sigma3,Sigma0,w1,w2] = CalculateHermiteRepresentationOfDensityOrderThree(w,Z,Prob,wnStart,Nn,dNn,NnMidh,NnMidl)

%This program find hermite representation of the SDE variable w to 
%tjird Order as a function of three hermite polynomials.
%This program is not perfectly precise and has very small errors so you
%will have to be careful using it everywhere.
%Here wnStart is the starting grid point of the SDE vector
%mostly wnStart=1.
%Please note we are representing SDE Variable w as a function of hermite 
%polynomials. It is not the density which can however be found by change of variable 
%derivative of w with respect to Z using standard gaussian variable.
%Nn is the index of last density grid point.
%So w(wnStart:Nn) represents the SDE density variable grid.
%  Z(wnStart:Nn) represents the associated value of standard Gaussian on
%  the grid. w and Z can be associated by matching appropriate CDF points
%  on the grid.
%I suppose the grid is uniform in Z and dNn is the spacing between and 
% Z(n) and Z(n+1) two points on the grid.
%Prob is the probability mass associated with each point on the grid.

wMean=sum(w(wnStart:Nn).*Prob(wnStart:Nn))  %Calculate the mean. This is in general 
%not on the discrete grid of w. 

%Below calculate the associated value of ZMean that maps wMean to ZMean.
%Also the grid point index nnMean which is associated with the largest
% w(nn) such that w(nn)<=wMean.
ZMean=(w(NnMidh)+w(NnMidl))/2.0;
nnMean=NnMidl;
for nn=wnStart:Nn-1
    if( (w(nn+1)>wMean) && (w(nn)<=wMean))
        ZMean=Z(nn)+(wMean-w(nn))/(w(nn+1)-w(nn))*(Z(nn+1)-Z(nn));
        nnMean=nn;
    end
end

%Below Calculate the vol associated with each point on the grid
for nn=wnStart:Nn
    Sigma0(nn)=(w(nn)-wMean)./(Z(nn)-ZMean);  
 end

%Below calculate the associated Sigma00 which is a single volatility associated with
%the first orthogonal hermite polynomial. We chose this volatility equal to
%volatility at mean since we are expanding around the mean.
%If mean has zero volatility, we will have to find volatility at mean by
%interpolating the volatility at one higher grid point and one lower grid
%point. Volatility at mean must not be zero (unless we have a delta density).
Sigma00=Sigma0(nnMean)*(Z(nnMean+1)-ZMean)./(Z(nnMean+1)-Z(nnMean))+Sigma0(nnMean+1)*(ZMean-Z(nnMean))./(Z(nnMean+1)-Z(nnMean));

%Below calculate the associated first derivative of Sigma0 with respect to
%Z on the grid.
dSigma0dZ(wnStart)=(-3*Sigma0(wnStart)+4*Sigma0(wnStart+1)-1*Sigma0(wnStart+2))/(2*dNn);
dSigma0dZ(wnStart+1:Nn-1)=(-1*Sigma0(wnStart:Nn-2)+1*Sigma0(wnStart+2:Nn))/(2*dNn);
dSigma0dZ(Nn)=(3*Sigma0(Nn)-4*Sigma0(Nn-1)+1*Sigma0(Nn-2))/(2*dNn);

%Following is the operation related to equation(4) in post 860.     
SSigma1(nnMean)=0;
for nn=nnMean+1:Nn
        SSigma1(nn)=SSigma1(nn-1)+(.5*dSigma0dZ(nn-1)+.5*dSigma0dZ(nn))*dNn;
end
for nn=nnMean-1:-1:wnStart
        SSigma1(nn)=SSigma1(nn+1)-(.5*dSigma0dZ(nn+1)+.5*dSigma0dZ(nn))*dNn;
end

%Below is the averaging
Sigma1(wnStart:nnMean-1)=SSigma1(wnStart:nnMean-1)./(Z(wnStart:nnMean-1)-Z(nnMean));
Sigma1(nnMean+1:Nn)=SSigma1(nnMean+1:Nn)./(Z(nnMean+1:Nn)-Z(nnMean));
Sigma1(nnMean)=(Sigma1(nnMean-1)+Sigma1(nnMean+1))/2.0;


Sigma11=Sigma1(nnMean)*(Z(nnMean+1)-ZMean)./(Z(nnMean+1)-Z(nnMean))+Sigma1(nnMean+1)*(ZMean-Z(nnMean))./(Z(nnMean+1)-Z(nnMean));



%Now below we represent the SDE variable w in terms of first two hermite
%polynomials and this associated value of w produced by hermite polynomials
%coefficients is called w1.
   
 w1(wnStart:Nn)=wMean-Sigma00*ZMean-Sigma1(wnStart:Nn).*(ZMean.^2-1)+ ...
 Sigma00*Z(wnStart:Nn)+ ...
 Sigma1(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
 


dSigma1dZ(wnStart)=(-3*Sigma1(wnStart)+4*Sigma1(wnStart+1)-1*Sigma1(wnStart+2))/(2*dNn);
dSigma1dZ(wnStart+1:Nn-1)=(-1*Sigma1(wnStart:Nn-2)+1*Sigma1(wnStart+2:Nn))/(2*dNn);
dSigma1dZ(Nn)=(3*Sigma1(Nn)-4*Sigma1(Nn-1)+1*Sigma1(Nn-2))/(2*dNn);

% Sigma2(NnMidh)=dSigma1dZ(NnMidh);
% for nn=NnMidh+1:Nn
%     Sigma2(nn)=(Sigma2(nn-1)*Z(nn-1)+(.5*dSigma1dZ(nn-1)+.5*dSigma1dZ(nn))*dNn)./Z(nn);
% end
% Sigma2(NnMidl)=dSigma1dZ(NnMidl);
% for nn=NnMidl-1:-1:wnStart
%     Sigma2(nn)=(Sigma2(nn+1)*Z(nn+1)-(.5*dSigma1dZ(nn+1)+.5*dSigma1dZ(nn))*dNn)./Z(nn);
% end

SSigma2(nnMean)=0;
for nn=nnMean+1:Nn
        SSigma2(nn)=SSigma2(nn-1)+(.5*dSigma1dZ(nn-1)+.5*dSigma1dZ(nn))*dNn;
end
for nn=nnMean-1:-1:wnStart
        SSigma2(nn)=SSigma2(nn+1)-(.5*dSigma1dZ(nn+1)+.5*dSigma1dZ(nn))*dNn;
end

Sigma2(wnStart:nnMean-1)=SSigma2(wnStart:nnMean-1)./(Z(wnStart:nnMean-1)-Z(nnMean));
Sigma2(nnMean+1:Nn)=SSigma2(nnMean+1:Nn)./(Z(nnMean+1:Nn)-Z(nnMean));
Sigma2(nnMean)=(Sigma2(nnMean-1)+Sigma2(nnMean+1))/2.0;

%Sigma22=Sigma2(nnMean);
Sigma22=Sigma2(nnMean)*(Z(nnMean+1)-ZMean)./(Z(nnMean+1)-Z(nnMean))+Sigma2(nnMean+1)*(ZMean-Z(nnMean))./(Z(nnMean+1)-Z(nnMean));

dSigma2dZ(wnStart)=(-3*Sigma2(wnStart)+4*Sigma2(wnStart+1)-1*Sigma2(wnStart+2))/(2*dNn);
dSigma2dZ(wnStart+1:Nn-1)=(-1*Sigma2(wnStart:Nn-2)+1*Sigma2(wnStart+2:Nn))/(2*dNn);
dSigma2dZ(Nn)=(3*Sigma2(Nn)-4*Sigma2(Nn-1)+1*Sigma2(Nn-2))/(2*dNn);


SSigma3(nnMean)=0;
for nn=nnMean+1:Nn
        SSigma3(nn)=SSigma3(nn-1)+(.5*dSigma2dZ(nn-1)+.5*dSigma2dZ(nn))*dNn;
end
for nn=nnMean-1:-1:wnStart
        SSigma3(nn)=SSigma3(nn+1)-(.5*dSigma2dZ(nn+1)+.5*dSigma2dZ(nn))*dNn;
end

Sigma3(wnStart:nnMean-1)=SSigma3(wnStart:nnMean-1)./(Z(wnStart:nnMean-1)-Z(nnMean));
Sigma3(nnMean+1:Nn)=SSigma3(nnMean+1:Nn)./(Z(nnMean+1:Nn)-Z(nnMean));
Sigma3(nnMean)=(Sigma3(nnMean-1)+Sigma3(nnMean+1))/2.0;





w2(wnStart:Nn)=wMean-Sigma00*ZMean-Sigma11.*(ZMean.^2-1)+ ...
 3*Sigma2(wnStart:Nn).*(ZMean.^3-3*ZMean) - ...
 3*Sigma2(wnStart:Nn).*Z(wnStart:Nn).*(ZMean.^2-1) + ... 
 Sigma00*Z(wnStart:Nn)+ ...
 Sigma11.*(Z(wnStart:Nn).^2-1)+ ...   
 0*Sigma2(wnStart:Nn).*(ZMean.^3-3*ZMean) + ...
 Sigma2(wnStart:Nn).*(Z(wnStart:Nn).^3-3*Z(wnStart:Nn));
 



%Below we plot a graphical comparison of the original SDE variable w and
%the variable w1 reproduced from the coefficients of the hermite
%polynomials
% plot(Z(wnStart:Nn),w(wnStart:Nn),'g',Z(wnStart:Nn),w1(wnStart:Nn),'r',Z(wnStart:Nn),w2(wnStart:Nn),'b')
%Sigma00
% str=input('Look at the graphical comparison of w and w1 as a function of standard gaussian Z');
end

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

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

February 17th, 2020, 1:28 pm

Dear friends, I have been thinking about the precise 'cause and effect' that we need to know to get our latest Ito-hermite algorithm working in which we evolve an expanding and drifting SD fractions grid. We divide this expanding and drifting SD fractions grid into subdivisions where center of each subdivision is associated/mapped on a fixed value of standard normal random variable. This SD fractions grid shares the fixed probability mass (even when it evolves) with the underlying standard normal distribution.
Please let me write the evolution SDE in Bessel transformed coordinates to first order for a general two drift SDE. If the original SDE is given as
[$]dy(t)=\mu_1 y^{\beta_1} dt + \mu_2 y^{\beta_2} dt + \sigma y^{\gamma} dz(t)[$] 
Then in transformed coordinates [$]w=\frac{y^{1-\gamma}}{(1-\gamma)}[$] our SDE becomes to first order of expansion
[$]dw(t)=\mu_1 y^{\beta_1-\gamma} dt + \mu_2 y^{\beta_2-\gamma} dt -.5 \gamma {\sigma}^2 y^{\gamma-1} dt  + \sigma dz(t)[$]
Our simple algorithm worked very well for CEV type noises and linear drift SDEs. The reason for these lucky cases as I understand goes as the following. For a proper and perfect algorithm, we needed to calculate the evolution for the entire SD grid subdivision. In particular, We needed to calculate the grid expansion effect integrated over the entire SD grid subdivisions. What we did was to calculate a one step projection effect applied to center(which is a single point) of each subdivision and used ad-hoc techniques that would somehow mimick the integrated grid expansion effect of entire local SD grid subdivision to improve the originally poor results. Luckily in bessel coordinates diffusion effect is very strongly linear so diffusion worked very well in our point evolution algorithm and this is one reason why we could not use the algorithm in original coordinates since in order to handle the SDEs in original coordinates, we would need to integrate the diffusion effect within each SD grid subdivisions. In the w coordinates the drift value would be proportional to [$] \mu_1 y^{\beta_1-\gamma} dt[$] so whenever [$]\beta_1[$] would be very close to [$]\gamma[$] the drift effect would be very linear and in such cases our point evolution algorithm worked perfectly well. 
Especially in case of non-linear drift that has non-linear variance effect in it, a right algorithm must not miss the drift effect over the entire subdivisions and we would need to integrate the drift effect over the SD grid subdivisions and use that to modify the point evolution algorithm.
writing our SDE evolution curve as a function of z, we can write as
[$]dw(z)=\mu_1 y(z)^{\beta_1-\gamma} dt + \mu_2 y(z)^{\beta_2-\gamma} dt -.5 \gamma {\sigma}^2 y(z)^{\gamma-1} dt  + \sigma dz(t)[$]
In our framework, We need to calculate for for each subdivision the drift effect as below. Here [$]z_n[$] is the value of z associated with the center of n-th subdivision

[$]\int_{z_n-.5\Delta z}^{z_n+.5\Delta z} (\mu_1 y(z)^{\beta_1-\gamma} dt + \mu_2 y(z)^{\beta_2-\gamma} dt-.5 \gamma {\sigma}^2 y^{\gamma-1} dt ) dz[$]

We can easily make this correction by taking a Taylor series expansion of the drift (with respect to [$]z_n[$]  around the center of the subdivision. We might have to take some numerical derivatives here). The first term would be of course the point value at the center which we were using in our point evolution algorithm. 2nd linear term would probably cancel around the center of the subdivision when we extend integration on both sides away from the center. However it is the third term with 2nd order convexity effect that we would need to integrate over each subdivision and then add to point drift to get modified drift that we would use in the update equation for the particular SD grid subdivision. We can appreciate that refining the grid in our point evolution algorithm framework would not have helped at all until we explicitly added the 2nd order derivative terms over each of the subdivisions.
Once it works for drift in bessel transformed case, I suppose we would be able to use similar technique for density evolution in original coordinates by applying it to diffusion term as well.

I will be coming back in another two or three days with a modified ito-hermite density update algorithm with the above corrections and present it here. 
 
User avatar
Amin
Topic Author
Posts: 2320
Joined: July 14th, 2002, 3:00 am

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

February 18th, 2020, 3:10 pm

I have been working how to calculate derivatives of SDE variable w with respect to Z analytically. Earlier I thought of following the numerical derivatives route but from my past bitter experience these derivatives are extremely bad at the boundaries and 2nd derivative when calculated numerically almost always makes the evolution equation unstable due to inaccuracies especially at the boundary. Therefore I decided to see if the derivatives [$]\frac{\partial w}{\partial Z}[$] and  [$]\frac{\partial^2 w}{\partial Z^2}[$] could be calculated analytically. If the SDE of w in bessel coordinates is given as
[$]dw(t)=\mu(w) dt + \sigma dz(t)[$]
Then I beileve that we could have heuristically written the evolution equation for the derivative [$]\frac{\partial w}{\partial Z}[$] along the nth SD-fraction as
[$]\frac{\partial w_n(t+1)}{\partial Z_n}=\frac{\partial w_n(t)}{\partial Z_n}+ \frac{\partial \mu(w_n)}{\partial w_n} \frac{\partial w_n(t)}{\partial Z_n} \Delta t + \frac{(\sqrt{(w_n-Mean(w))^2 + (\sigma \sqrt{\Delta t} Z_n)^2} - (w_n-Mean(w))}{Z_n}[$]
The above equation is based on knowledge learnt from our main update equation I have been using to add volatility. We will have to be careful about signs and I have not added anything about signs in the above. And also we will have the initial condition at t=0, [$]\frac{\partial w_n(t=0)}{\partial Z_n}=0[$]
and also 
evolution of second order derivative equation could be written as
[$]\frac{\partial^2 w_n(t+1)}{\partial {Z_n}^2}=\frac{\partial^2 w_n(t)}{\partial {Z_n}^2}+ \frac{\partial^2 \mu(w_n)}{\partial {w_n}^2} (\frac{\partial w_n(t)}{\partial Z_n})^2 \Delta t + \frac{\partial \mu(w_n)}{\partial {w_n}} (\frac{\partial^2 w_n(t)}{\partial {Z_n}^2}) \Delta t[$]

At each update step, we would use these derivatives to calculate the drift correction that is needed in order to integrate the drift over the entire SD-fractions as opposed to point calculations and then use the corrected value of w in future evolution and derivatives calculation.
I have not tested these equations but would be testing them hopefully tomorrow or possibly later tonight. But decided to share with friends. There could however possibly be some errors but I hope these equations are right.
 
User avatar
Amin
Topic Author
Posts: 2320
Joined: July 14th, 2002, 3:00 am

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

February 19th, 2020, 8:07 pm

Ok friends, here is the new progress. Integrating over the subdivisions does not seem to help at all. But when I calculated [$]\frac{\partial \mu(w)}{\partial z}[$] and  [$]\frac{\partial^2 \mu(w)}{\partial z^2}[$] analytically, and tried playing to find an optimal value that could give a good fit to the densities, I found that subtracting the following term greatly improves the fit of the density almost everywhere. Below has to be added to the main SDE variable update equation.
[$]-\frac{1}{(\sqrt{2} {\pi}^2)} \big [ \frac{\partial^2 \mu(w)}{\partial w^2} {(\frac{\partial w}{\partial Z})}^2 - \frac{\partial \mu(w)}{\partial w} {(\frac{\partial^2 w}{\partial Z^2})} \big](Z^2-1)[$]
Shape of the density is universally captured almost everywhere to great accuracy but there is very slight inaccuracy in the tail. My analytical calculations of [$]\frac{\partial \mu(w)}{\partial z}[$] and  [$]\frac{\partial^2 \mu(w)}{\partial z^2}[$] are only approximate and I will try to improve them in coming days. I am posting a program with the above changes included. I have disabled mean correction but when you would enable mean correction, you would have completely exact fit to the density.
function [] = ItoHermiteMethodWilmott08ANewGood()
%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 coordinates.
%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;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.

Tt=128;%16*16;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year; 
OrderA=2;  %Keep at two.
OrderM=2;  %Keep at two.
%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;

x0=.1250;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.850;   % volatility power.                                                                                                                                                                                                                                                                           
kappa=1.4850;   %mean reversion parameter.
theta=.1250;   %mean reversion target

%you can specify any general mu1 and mu2 and beta1 and beta2.
%The present program is not perfect but quite reasonably good.
%if there is only one term in drift with beta=0, there will be minor
%errors. This case is yet not covered.

mu1=+1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
sigma0=.80;%Volatility value
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

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



 dt2=dt^2/2;
 dt3=dt^3/6.0;
 dtsqrt=sqrt(dt);
 dtonehalf=dt.^1.5;
 
wnStart=1;%
d2wdZ2(1:Nn)=0;
dwdZ0(1:Nn)=0;
dwdZ(1:Nn)=0;
dwdZQ(1:Nn)=0;

tic

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

    theta1=mu1;
    theta2=mu2;
    theta3=.5*sigma0^2*(-gamma);
    omega1=beta1-gamma;
    omega2=beta2-gamma;
    omega3=gamma-1;
    
    GG(wnStart:Nn)=theta1*yy(wnStart:Nn).^omega1+theta2*yy(wnStart:Nn).^omega2+theta3*yy(wnStart:Nn).^omega3;
    
    dGGdw(wnStart:Nn)=omega1* theta1 *((1-gamma)* w(wnStart:Nn)).^(-1+omega1/(1-gamma))+ ...
        omega2* theta2* ((1-gamma)* w(wnStart:Nn)).^(-1+omega2/(1-gamma))+ ...
        omega3* theta3* ((1-gamma)* w(wnStart:Nn)).^(-1+omega3/(1-gamma));
    
    
    d2GGdw2(wnStart:Nn)=(1 - gamma)* omega1* (-1 + omega1/(1 - gamma))* theta1* ((1 - gamma)* ...
        w(wnStart:Nn)).^(-2 + omega1/(1 - gamma)) + ...
        (1 - gamma)* omega2* (-1 + omega2/(1 - gamma))* theta2* ((1 - gamma)* w(wnStart:Nn)).^(-2 + omega2/(1 - gamma)) + ...
        (1 - gamma)* omega3* (-1 + omega3/(1 - gamma))* theta3* ((1 - gamma)* w(wnStart:Nn)).^(-2 + omega3/(1 - gamma));
    
    
    dGG(wnStart:Nn)=omega1*theta1*yy(wnStart:Nn).^(omega1-1)+omega2*theta2*yy(wnStart:Nn).^(omega2-1) + ...
        +omega3*theta3*yy(wnStart:Nn).^(omega3-1);
    
    ddGGdw(wnStart:Nn)=(-1+omega1)* omega1* theta1* ((1-gamma)* w(wnStart:Nn)).^(-1+(-1+omega1)/(1-gamma))+ ...
        (-1+omega2)* omega2* theta2* ((1-gamma)* w(wnStart:Nn)).^(-1+(-1+omega2)/(1-gamma))+ ...
        (-1+omega3)* omega3* theta3* ((1-gamma)* w(wnStart:Nn)).^(-1+(-1+omega3)/(1-gamma));
    
     d2dGGdw2(wnStart:Nn)=(1-gamma)* (-1+(-1+omega1)/(1-gamma))* (-1+omega1)* omega1 * ...
         theta1* ((1-gamma)* w(wnStart:Nn)).^(-2+(-1+omega1)/(1-gamma))+ ...
         (1-gamma)* (-1+(-1+omega2)/(1-gamma))* (-1+omega2)* omega2 *theta2* ((1-gamma)* w(wnStart:Nn)).^(-2+(-1+omega2)/(1-gamma))+ ...
         (1-gamma)* (-1+(-1+omega3)/(1-gamma))* (-1+omega3)* omega3 *theta3* ((1-gamma)* w(wnStart:Nn)).^(-2+(-1+omega3)/(1-gamma));
     
    DD(wnStart:Nn)=mu1*yy(wnStart:Nn).^beta1+mu2*yy(wnStart:Nn).^beta2;
    
    dDDdw(wnStart:Nn)=beta1* mu1* ((1 - gamma)* w(wnStart:Nn)).^(-1 + beta1/(1 - gamma)) + ... 
        beta2* mu2* ((1 - gamma)* w(wnStart:Nn)).^(-1 + beta2/(1 - gamma));
    d2DDdw2(wnStart:Nn)=beta1* (-1 + beta1/(1 - gamma))* (1 - gamma)* mu1* ((1 - gamma)* w(wnStart:Nn)).^(-2 + beta1/(1 - gamma)) + ...
        beta2* (-1 + beta2/(1 - gamma))* (1 - gamma)* mu2* ((1 - gamma)* w(wnStart:Nn)).^(-2 + beta2/(1 - gamma));
    
    dGG2(wnStart:Nn)=theta1*omega1*(omega1-1).*yy(wnStart:Nn).^(omega1-2)+ ...
        theta2*omega2*(omega2-1).*yy(wnStart:Nn).^(omega2-2)+ ...
        theta3*omega3*(omega3-1).*yy(wnStart:Nn).^(omega3-2);
    
    ddGG2dw(wnStart:Nn)=(-2+omega1)* (-1+omega1)* omega1* theta1 *((1-gamma)* w(wnStart:Nn)).^(-1+(-2+omega1)/(1-gamma))+ ...
        (-2+omega2)* (-1+omega2)* omega2* theta2* ((1-gamma)* w(wnStart:Nn)).^(-1+(-2+omega2)/(1-gamma))+ ...
        (-2+omega3)* (-1+omega3)* omega3* theta3* ((1-gamma)* w(wnStart:Nn)).^(-1+(-2+omega3)/(1-gamma));
    
    d2dGG2dw2(wnStart:Nn)=(1 - gamma)* (-1 + (-2 + omega1)/(1 - gamma))* (-2 + omega1)* ...
        (-1 + omega1)* omega1* theta1* ((1 - gamma)* w(wnStart:Nn)).^(-2 + (-2 + omega1)/(1 - gamma)) + ...
        (1 - gamma)* (-1 + (-2 + omega2)/(1 - gamma))* (-2 + omega2)* (-1 + omega2)* omega2 * ...
        theta2* ((1 - gamma)* w(wnStart:Nn)).^(-2 + (-2 + omega2)/(1 - gamma)) + ...
        (1 - gamma)*(-1 + (-2 + omega3)/(1 - gamma))* (-2 + omega3)* (-1 + omega3)* omega3 * ...
        theta3* ((1 - gamma)* w(wnStart:Nn)).^(-2 + (-2 + omega3)/(1 - gamma));
    
    QQ(wnStart:Nn)=sigma0^2*yy(wnStart:Nn).^(2*gamma);
    
    dQQdw(wnStart:Nn)=2* gamma* sigma0^2* ((1 - gamma)* w(wnStart:Nn)).^(-1 + (2 *gamma)/(1 - gamma));
    
    d2QQdw2(wnStart:Nn)=2* (1 - gamma)* gamma* (-1 + (2* gamma)/(1 - gamma)) * ...
        sigma0^2* ((1 - gamma)* w(wnStart:Nn)).^(-2 + (2* gamma)/(1 - gamma));
    
    ddGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*yy(wnStart:Nn).^(omega1+beta1-2)+ ...
        mu1*omega2*theta2*(omega2+beta1-1).*yy(wnStart:Nn).^(omega2+beta1-2)+ ...
        mu1*omega3*theta3*(omega3+beta1-1).*yy(wnStart:Nn).^(omega3+beta1-2)+ ...
        mu2*omega1*theta1*(omega1+beta2-1).*yy(wnStart:Nn).^(omega1+beta2-2)+ ...
        mu2*omega2*theta2*(omega2+beta2-1).*yy(wnStart:Nn).^(omega2+beta2-2)+ ...
        mu2*omega3*theta3*(omega3+beta2-1).*yy(wnStart:Nn).^(omega3+beta2-2);
    
    d2dGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*(omega1+beta1-2).*yy(wnStart:Nn).^(omega1+beta1-3)+ ...
        mu1*omega2*theta2*(omega2+beta1-1).*(omega2+beta1-2).*yy(wnStart:Nn).^(omega2+beta1-3)+ ...
        mu1*omega3*theta3*(omega3+beta1-1).*(omega3+beta1-2).*yy(wnStart:Nn).^(omega3+beta1-3)+ ...
        mu2*omega1*theta1*(omega1+beta2-1).*(omega1+beta2-2).*yy(wnStart:Nn).^(omega1+beta2-3)+ ...
        mu2*omega2*theta2*(omega2+beta2-1).*(omega2+beta2-2).*yy(wnStart:Nn).^(omega2+beta2-3)+ ...
        mu2*omega3*theta3*(omega3+beta2-1).*(omega3+beta2-2).*yy(wnStart:Nn).^(omega3+beta2-3);

    ddGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*yy(wnStart:Nn).^(omega1+2*gamma-3)+ ...
        sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*yy(wnStart:Nn).^(omega2+2*gamma-3)+ ...
        sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*yy(wnStart:Nn).^(omega3+2*gamma-3);
        
    d2dGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*(omega1+2*gamma-3).*yy(wnStart:Nn).^(omega1+2*gamma-4)+ ...
        sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*(omega2+2*gamma-3).*yy(wnStart:Nn).^(omega2+2*gamma-4)+ ...
        sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*(omega3+2*gamma-3).*yy(wnStart:Nn).^(omega3+2*gamma-4);
    
    VV(wnStart:Nn)=sigma0.*yy(wnStart:Nn).^gamma;
    d_dGG_VV_dyy(wnStart:Nn)= ...
        sigma0.*omega1.*theta1.*(omega1+gamma-1).*yy(wnStart:Nn).^(omega1+gamma-2) + ...
        sigma0.*omega2.*theta2.*(omega2+gamma-1).*yy(wnStart:Nn).^(omega2+gamma-2) + ...
        sigma0.*omega3.*theta3.*(omega3+gamma-1).*yy(wnStart:Nn).^(omega3+gamma-2);    

    wMu0dt(wnStart:Nn)=GG(wnStart:Nn)*dt + dGG(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...
        .5*dGG2(wnStart:Nn).*QQ(wnStart:Nn).*dt2 + ...
        ddGGDD(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
        .5*ddGG2QQ(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
        .5*d2dGGDD(wnStart:Nn).*QQ(wnStart:Nn).*dt3 + ...
        .25*d2dGG2QQ(wnStart:Nn).*QQ(wnStart:Nn).*dt3;
    dwMu0dtdw(wnStart:Nn)=dGGdw(wnStart:Nn)*dt + ddGGdw(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...
        dGG(wnStart:Nn).*dDDdw(wnStart:Nn).*dt2+ ...
        .5*ddGG2dw(wnStart:Nn).*QQ(wnStart:Nn).*dt2 + ...
        .5*dGG2(wnStart:Nn).*dQQdw(wnStart:Nn).*dt2;
    d2wMu0dtdw2(wnStart:Nn)=d2GGdw2(wnStart:Nn)*dt + d2dGGdw2(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...
        2*ddGGdw(wnStart:Nn).*dDDdw(wnStart:Nn).*dt2+ ...
        dGG(wnStart:Nn).*d2DDdw2(wnStart:Nn).*dt2+ ...
        .5*d2dGG2dw2(wnStart:Nn).*QQ(wnStart:Nn).*dt2 + ...
        .5*2*ddGG2dw(wnStart:Nn).*dQQdw(wnStart:Nn).*dt2 + ...
        .5*dGG2(wnStart:Nn).*d2QQdw2(wnStart:Nn).*dt2;
    
    
    
    %    c0(wnStart:Nn)=-1*d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0;
    c1(wnStart:Nn)=sigma0.*dtsqrt+dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*dtonehalf;
    c2(wnStart:Nn)=d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0;
    dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn);% + ...
    dw2(wnStart:Nn)=c1(wnStart:Nn).^2.0.*Z(wnStart:Nn).^2;
    w(isnan(w)==1)=0;
    wPrev(wnStart:Nn)=w(wnStart:Nn);
    wMu0dt(isnan(wMu0dt)==1)=0;  
    wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%This is the main SDE variable one step update equation. 


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)))- ...
          1/sqrt(2)/(22/7)^2*(d2wMu0dtdw2(wnStart:Nn).*dwdZ(wnStart:Nn).^2-dwMu0dtdw(wnStart:Nn).*d2wdZ2(wnStart:Nn)).*(Z(wnStart:Nn).^2-1);
wMeanNew=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));

d2wdZ2(wnStart:Nn)=d2wdZ2(wnStart:Nn)+ d2wMu0dtdw2(wnStart:Nn).*dwdZ(wnStart:Nn).^2+dwMu0dtdw(wnStart:Nn).*d2wdZ2(wnStart:Nn);

dwdZ0(wnStart:Nn)=dwdZ0(wnStart:Nn)+ dwMu0dtdw(wnStart:Nn).*dwdZ0(wnStart:Nn)+ ...
    ((sign(wPrev(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
          sqrt(abs(sign(wPrev(wnStart:Nn)-wMeanPrev).*(wPrev(wnStart:Nn)-wMeanPrev).^2+ ...
          sign(dw(wnStart:Nn)).*dw2(wnStart:Nn))))-(wPrev(wnStart:Nn)-wMeanPrev))./Z(wnStart:Nn);
      
dwdZQ(wnStart:Nn)=((w(wnStart:Nn)-wMeanNew))./Z(wnStart:Nn);      
 
   dwdZN(wnStart+1:Nn-1)=(w(wnStart+2:Nn)-w(wnStart:Nn-2))/(2*dNn);
    dwdZN(wnStart)= InterpolateOrderN(4,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),0,dwdZN(wnStart+1),dwdZN(wnStart+2),dwdZN(wnStart+3),dwdZN(wnStart+4),0);
    dwdZN(Nn)= InterpolateOrderN(4,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),0,dwdZN(Nn-1),dwdZN(Nn-2),dwdZN(Nn-3),dwdZN(Nn-4),0);

 dwdZ(wnStart:Nn)=(.5*dwdZ0(wnStart:Nn)+.5*dwdZQ(wnStart:Nn));
    
    
    
%    plot(Z(wnStart+1:Nn-1),dwdZ(wnStart+1:Nn-1),'k',Z(wnStart+1:Nn-1),dwdZ0(wnStart+1:Nn-1),'r',Z(wnStart+1:Nn-1),dwdZN(wnStart+1:Nn-1),'g',Z(wnStart+1:Nn-1),dwdZQ(wnStart+1:Nn-1),'b'); 
    
    
%    str=input('Look at dwdz');
    
     
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
%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.
    CTheta=1.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+2)+3*w(NnMidh+3)-3/2*w(NnMidh+4)+1/3*w(NnMidh+5))/dNn;
    D2wh=(35/12*w(NnMidh+2) -26/3*w(NnMidh+3)+ 19/2*w(NnMidh+4)-14/3*w(NnMidh+5)+11/12*w(NnMidh+6))/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);



             
        w1(1:Nn-1)=w(1:Nn-1);
        w2(1:Nn-1)=w(2:Nn);
        w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
     
        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
str=input('Analytic Values calculated; Press a key to start monte carlo');


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.
%Pre-compute the time and power exponent values in small multi-dimensional arrays
for k = 0 : (OrderA+1)
    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));

            end
        end
    end
end
%Below call the program that calculates the Ito-Taylor coeffs.
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);





%Yq=ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);
%YqCoeff=Yq/(1-gamma);%These coordinates are divided by (1-gamma)
 %Calculate the  density with monte carlo below.
 %rng(85109634, 'twister')
 rng(29079137, 'twister')
 paths=250000;
 YY(1:paths)=x0;  %Original process monte carlo.
% YY2dt(1:paths)=0.0;
 %YY4dt(1:paths)=0.0;
 X(1:paths)=x0^(1-gamma)/(1-gamma);
 Random1(1:paths)=0;
 YYMean(1:TtM)=0;
 for tt=1:TtM
 %    t1=tt*dtM;
 %    t0=(tt-1)*dtM;
     Random1=randn(size(Random1));
     HermiteP1(1,1:paths)=1;
     HermiteP1(2,1:paths)=Random1(1:paths);
     HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
     HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
     HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;
     YY0=YY;
     YY1=YY;
  %   YYPrev=YY;
  %   XX0=X;
     
     for k = 0 : OrderM
         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;
                     if(l4<=5)   
                         
  %                    X(1:paths)=X(1:paths) + YqCoeff(l1,l2,l3,l4) * ...
  %                   mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
  %                   ((1-gamma)*XX0(1:paths)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* ...
  %                   HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);    
                 
                     YY(1:paths)=YY(1:paths) + YCoeff(l1,l2,l3,l4) * ...
                     mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
                     YY0(1:paths).^Fp(l1,l2,l3,l4) .* HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);
%                     YY(YY>0)=YY1(YY>0);
%                     YY(YY<0)=0;
                 %Above is original coordinates monte carlo equation.
                     end
                 end
             end
         end
     end
     %YY(YY<0)=0;
     YYMean(tt)=YYMean(tt)+sum(YY(1:paths))/paths;
 end
 YY(YY<0)=0;
 disp('Original process average from monte carlo');
 sum(YY(:))/paths %origianl coordinates monte carlo average.
 disp('Original process average from our simulation');
 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');
 theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
 
 
 MaxCutOff=30;
 NoOfBins=round(500*gamma^2*1*sigma0);
 %[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
 [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');
 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
The program requires this new dependency function
function [y0] = InterpolateOrderN(N,x0,x1,x2,x3,x4,x5,y1,y2,y3,y4,y5)


if(N==5)
y0=(x0-x2)*(x0-x3)*(x0-x4)*(x0-x5)/((x1-x2).*(x1-x3).*(x1-x4).*(x1-x5))*y1+ ...
   (x0-x1)*(x0-x3)*(x0-x4)*(x0-x5)/((x2-x1).*(x2-x3).*(x2-x4).*(x2-x5))*y2+ ...
   (x0-x1)*(x0-x2)*(x0-x4)*(x0-x5)/((x3-x1).*(x3-x2).*(x3-x4).*(x3-x5))*y3+ ...
   (x0-x1)*(x0-x2)*(x0-x3)*(x0-x5)/((x4-x1).*(x4-x2).*(x4-x3).*(x4-x5))*y4+ ...
   (x0-x1)*(x0-x2)*(x0-x3)*(x0-x4)/((x5-x1).*(x5-x2).*(x5-x3).*(x5-x4))*y5;
end
if(N==4)
y0=(x0-x2)*(x0-x3)*(x0-x4)/((x1-x2).*(x1-x3).*(x1-x4))*y1+ ...
   (x0-x1)*(x0-x3)*(x0-x4)/((x2-x1).*(x2-x3).*(x2-x4))*y2+ ...
   (x0-x1)*(x0-x2)*(x0-x4)/((x3-x1).*(x3-x2).*(x3-x4))*y3+ ...
   (x0-x1)*(x0-x2)*(x0-x3)/((x4-x1).*(x4-x2).*(x4-x3))*y4;
end
if(N==3)
y0=(x0-x2)*(x0-x3)/((x1-x2).*(x1-x3))*y1+ ...
   (x0-x1)*(x0-x3)/((x2-x1).*(x2-x3))*y2+ ...
   (x0-x1)*(x0-x2)/((x3-x1).*(x3-x2))*y3;
end
if(N==2)
y0=(x0-x2)/((x1-x2))*y1+ ...
   (x0-x1)/((x2-x1))*y2;
end

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

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

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

February 19th, 2020, 9:23 pm

Here is the graph output comparison of true density with Ito-Hermite method density with parameters as in the posted program.

[$]\Delta t=.125/16, Tt=128,T=1, x_0=.1250, \beta_1=0.0, \beta_2=1.0, \kappa=1.4850,\theta=.125,\sigma=.8,\gamma=.85[$]



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

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

February 21st, 2020, 4:18 am

While I will continue to work for a genuinely better analytical approximation of derivatives of w with respect to Z. Here is something very simple that improves the above derivatives quite a bit.

please replace the line of code
dwdZ(wnStart:Nn)=(.5*dwdZ0(wnStart:Nn)+.5*dwdZQ(wnStart:Nn));
with the following two lines of code
dwdZ(wnStart:NnMidl)=(.4*dwdZ0(wnStart:NnMidl)+.6*dwdZQ(wnStart:NnMidl));
dwdZ(NnMidh:Nn)=(.6*dwdZ0(NnMidh:Nn)+.4*dwdZQ(NnMidh:Nn));

and you will have somewhat better results.
 
User avatar
Amin
Topic Author
Posts: 2320
Joined: July 14th, 2002, 3:00 am

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

February 24th, 2020, 3:38 pm

Friends here is the progress on new work. The above algorithm works fine for evolving the densities of SDEs from Z=-5 to Z=1.0 (even in this region I made a slight improvement); However Z>1.0 takes a distortion in tail proportional to total effective volatility of the SDE. I have been working on ways to correct the right tail. I have had very good progress and I hope withing 2-3 days, I will have an algorithm that will find the density of SDEs to great accuracy with calculated mean of the evolved density within a few basis points of the true mean even for high volatilities. I hope to post the new program within a few days.
 
User avatar
Amin
Topic Author
Posts: 2320
Joined: July 14th, 2002, 3:00 am

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

February 28th, 2020, 6:03 pm

Dear friends it will still be several days before a completely worked out algorithm for the evolution of SDEs using Ito-Hermite method could be available. In the meantime, I decided to share some thoughts with the friends about the intuition of how drift contributes to variance in transformed Bessel/Lamperti coordinates. Let us suppose we have a general SDE with two drift terms given as
[$]dy(t)=\mu_1 y(t)^{\beta_1} dt + \mu_2 y(t)^{\beta_2} dt + \sigma y(t)^{\gamma} dz(t)[$]

To avoid confusion, I want to mention that original SDE variable is [$]y(t)[$] and its transformed version is given as [$]w(t)=\frac{y(t)^{(1-\gamma)}}{(1-\gamma)}[$]
and evolution of [$]w(t)[$] is given as
 [$]dw(t)=\mu_1 y(t)^{\beta_1-\gamma} dt + \mu_2 y(t)^{\beta_2-\gamma} dt -.5 \gamma {\sigma}^2  y(t)^{\gamma-1} dt+ \sigma dz(t)[$]

Please note that noise term in the transformed SDE of w(t) is linear but it is the drift that is making a contribution to variance and it is the higher order contribution effect of drift that makes the transformed SDE density different from normal density. The first three dt terms in the above transformed SDE are what I refer to drift in my comments below and denote as [$]\mu(t)[$] in following. Here are a few observations

1. In our framework of SD fractions based on a particular value of standard normal Z, if the drift of any SD fraction is different from the drift of the mean SD-fraction, it will make a contribution to local variance of the SDE. Only when drift of all SD-fractions are identically the same as the drift of the mean SD-fraction that there will be no contribution from the drift to local variance of the SDE. Since the diffusion term is linear of first order all non-linearity/effects of higher order in the variance of density are caused by different local drift along the density or along each SD-fraction(with unique value of Z associated with each SD-fraction/subdivision) in our framework.

2. To conceptually understand the effect of drift on variance, let us suppose that [$]\frac{\partial \mu(t)}{\partial Z}=constant[$] and other higher derivatives of [$]\mu(t)[$] with respect to Z along SD-fractions are equal to zero, such a drift will make an equal and constant first order contribution to local variance all along the SDE (or in our case make a constant contribution to variance of each SD-fraction associated with different particular value of Z). And since the contribution of diffusion term is also linear of first order, and the drift term also contributes to diffusion/variance only in first linear order, higher order hermite polynomials will be zero and our density of w(t) will be very close to scaled density of Z, the normal random variable.

3. since diffusion term is linear of first order in  SDE of w(t) and does not directly affect the higher order variance, the higher order variance effects in the density of w(t) are caused by higher order derivatives of drift with respect  to Z along the SD-fractions. These derivatives are given as [$]\frac{\partial^2 \mu(t)}{\partial Z^2}[$], [$]\frac{\partial^3 \mu(t)}{\partial Z^3}[$] and so on. If we could easily calculate these derivatives, we could calculate the effect of drift by adding the derivatives of drift (with respect to Z) to the existing density after multiplying each higher order derivative of drift with hermite polynomial of the same degree(and accounting for appropriate scaling coefficients). 

Please also notice that it was second order derivatives of drift with respect to Z that were able to make a better fit to the density in my last program when added after multiplication with 2nd order hermite polynomial.
 
User avatar
Amin
Topic Author
Posts: 2320
Joined: July 14th, 2002, 3:00 am

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

March 1st, 2020, 11:44 am

Dear friends, my work is still not complete and still may take a few days. But here is the line of research I have taken to correct the SDE densities and the ideas I am currently working on. We added a term with second derivative of drift with respect to Z and it worked extremely well but the tail remained problematic. Now, I believe there need to be possibly one or both first and third order correction terms with corresponding first and third order derivatives of drift with respect to Z and multiplied with appropriate hermite polynomials. And this is what I am currently working with to fix the program. I hope to be able to complete the work in a few days. 
Previously I had taken a different approach but it did not seem to work well into time after a few years so I had to abandon that.
 
User avatar
Amin
Topic Author
Posts: 2320
Joined: July 14th, 2002, 3:00 am

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

March 2nd, 2020, 5:33 pm

I have not been able to work very well for past two days but I was still able to determine that it is the third derivative correction that has to be further after the second derivative correction(as given in the previous program). After the third derivative correction, we can get very exact density if the mean is known and mean correction could be applied. I would be posting this program in two to three days after testing it well. Though shape of the density and fit is quite excellent for most SDEs even if the mean correction is not applied but the density mean in such cases could be slightly off from true mean. After posting the new program in 2-3 days, I will continue work on precision calculation of derivatives of drift with Z,  and exact mean for every SDE. There is still a long list of things that still remain to be done like, for example, introducing transition probabilities framework on the SD-fractions grid and work on arithmetic path integrals of SDEs using the SD-fractions and transition probabilities together. I would also want to extend the work to two dimensional SDEs for asset pricing with stochastic volatility. So I hope to add a lot of material to this thread in coming weeks and months.
 
User avatar
Amin
Topic Author
Posts: 2320
Joined: July 14th, 2002, 3:00 am

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

March 4th, 2020, 8:14 pm

Friends, it will still be a few days before I could post a final program. But there is some good progress and I am posting an informal update here. I am writing the term that seems to work well to correct the density though I expect some further minor changes in this term but I am sure the final correction will be some variation on the same theme I am posting here now. For friends who want to play with the program on their own here is the term you need to add to main density update equation.

---------------------------------------------------------------------------------------------------------------
-.225*yyMean.^(-2+2*gamma)/(22/7)^2.*sign(CTerm(wnStart:Nn)).*sqrt(abs(CTerm(wnStart:Nn))).*H1Zero0(wnStart:Nn).*sqrt(dt)

In the above term CTerm is defined as
CTerm(wnStart:Nn)=(d2wMu0dtdw2(wnStart:Nn).*dwdZ(wnStart:Nn).^2 - dwMu0dtdw(wnStart:Nn).*d2wdZ2(wnStart:Nn));
and H1Zero0(wnStart:Nn) is defined as

H1Zero0(1:Nn)=0;
H1Zero0(NnMidh:Nn)=Z(NnMidh:Nn);

and yyMean is defined as 

yyMean=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn));

------------------------------------------------------------------------------------------------------------

Please note that H1Zero0 is first hermite polynomial defined only in positive half domain and zero in negative domain. May be after skew is added, the first hermite polynomial here could be defined both for positive and negative domains. But I am very sure that the above is the right term and authentic final correction would be a slight variation on the same theme.  If you add the above term and then apply(enable) mean correction, you will be able to simulate the density almost really perfectly.
 
User avatar
Amin
Topic Author
Posts: 2320
Joined: July 14th, 2002, 3:00 am

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

March 5th, 2020, 5:11 am

While I will continue to work for a genuinely better analytical approximation of derivatives of w with respect to Z. Here is something very simple that improves the above derivatives quite a bit.

please replace the line of code
dwdZ(wnStart:Nn)=(.5*dwdZ0(wnStart:Nn)+.5*dwdZQ(wnStart:Nn));
with the following two lines of code
dwdZ(wnStart:NnMidl)=(.4*dwdZ0(wnStart:NnMidl)+.6*dwdZQ(wnStart:NnMidl));
dwdZ(NnMidh:Nn)=(.6*dwdZ0(NnMidh:Nn)+.4*dwdZQ(NnMidh:Nn));

and you will have somewhat better results.
Sorry that I failed to update with some changes. Right after the above change, please add these lines of code to your program as this will make dwdZ nice and smooth in the middle.

dwdZ(NnMidl-2)=6/7*dwdZ(NnMidl-3)+1/7*dwdZ(NnMidh+3);
 dwdZ(NnMidl-1)=5/7*dwdZ(NnMidl-3)+2/7*dwdZ(NnMidh+3);
 dwdZ(NnMidl)=  4/7*dwdZ(NnMidl-3)+3/7*dwdZ(NnMidh+3);
 dwdZ(NnMidh)=  3/7*dwdZ(NnMidl-3)+4/7*dwdZ(NnMidh+3);
 dwdZ(NnMidh+1)=2/7*dwdZ(NnMidl-3)+5/7*dwdZ(NnMidh+3);
 dwdZ(NnMidh+2)=1/7*dwdZ(NnMidl-3)+6/7*dwdZ(NnMidh+3);


I will come up with a reworked and better version of derivatives  dwdZ and d2wdZ2 in a few days.
 
User avatar
Amin
Topic Author
Posts: 2320
Joined: July 14th, 2002, 3:00 am

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

March 7th, 2020, 11:56 am

Post deleted. I will come back with a reworked program in a few days.
 
User avatar
Amin
Topic Author
Posts: 2320
Joined: July 14th, 2002, 3:00 am

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

March 8th, 2020, 11:04 pm

I was thinking of trying something on the lines of the following equation for the main update equation

[$] w(t+\Delta t)=\mu_w +[ \mu(t) \Delta t - .5 \frac{\partial^2 \mu(t) \Delta t}{\partial Z^2}  H_2(Z)] + \sqrt{(([w(t)-.5 \frac{\partial^2 w(t)}{\partial Z^2} H_2(Z)]-\mu_w )^2 + (\Delta w)^2)}[$]
[$]+ .5 \sqrt{(( \frac{\partial^2 w(t)}{\partial Z^2})^2 + ( \frac{\partial^2 \mu(t) \Delta t}{\partial Z^2})^2)}  H_2(Z) [$]

This is just a rough hypothesis at the moment but I will try this and its variants tomorrow. I decided  to share it as a food for thought for friends so they could try their own variants if they find something appropriate in this idea. I will share the results with friends how it goes tomorrow.
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