SERVING THE QUANTITATIVE FINANCE COMMUNITY

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

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

August 25th, 2019, 5:40 pm

My second order correction to evolution of density of non-linear mean-reverting SDEs is only approximate but is a huge improvement on first order and gives a very good idea about structure of the optimal correction. Please go ahead and try to think of a better correction but please do not forget to give me credit for this work.

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)
where mu1=kappa*theta and mu2=-kappa with  beta1=0 and beta2=1;
I have not directly simulated the SDE but simulated the transformed  Bessel 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. The present program is meant for mean-reversion SDEs that take two terms in non-linear drift. It will work for simpler SDEs but for that you will  have to turn SecondOC variable, second order correction variable to zero.
For diffusions going into zero, you may see instabilities, I have not corrected for that yet in this initial version. Heston is still not good as substantial density goes into zero at practically realistic vol. Heston when density does not hit zero will work perfectly.
I try volatility-power=gamma above .65 or .70; gamma=.95 is proxy for lognormal, the best and coveted SV model but not used due to lack of tractibility.
Heston and even smaller gamma will work far better when x0=1 and theta=1 as in Interest rate LMM like models with local volatility.
Very close to zero, heston density simulation does not work very well with this program. Whenever substantial  density goes into zero, performance of the current program deteriorates.
 
User avatar
Amin
Topic Author
Posts: 2122
Joined: July 14th, 2002, 3:00 am

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

August 26th, 2019, 5:34 pm

I have removed the previous program that I posted yesterday as it was very suboptimal compared to the corrected version. I am posting a new good program now that works excellently. Now Heston is also working with my new program in many many cases even on higher volatilities. I am sure friends would love this program. Though there is need for more theoretical research, I doubt if the program below can altered significantly now since it is already working very well. Here is the only line of code with the second order correction that I newly added to the  existing program 

SecondOC(wnStart:Nn)=-.0027*((tt-1)*dt).*dt./gamma^(2*gamma-1).*sigma0^4./(w(wnStart:Nn)*(1-gamma)).*(Z(wnStart:Nn).^2-1);

Here is the new  code. Most diffusions going into zero are working perfectly fine now. I will be posting some graphs produced with this code tomorrow. 
function [] = ItoTaylorDensityMRRecombining14F2Wilmott()
%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)
%where mu1=kappa*theta and mu2=-kappa with
%beta1=0 and beta2=1;
%I have not directly simulated the SDE but simulated the transformed 
%Bessel 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.
%The present program is meant for mean-reversion SDEs that take two terms
%in non-linear drift. It will work for simpler SDEs but for that you will 
%have to turn SecondOC variable, second order correction variable to zero.
%
%For diffusions going into zero, you may see instabilities, I have not
%corrected for that yet in this initial version.
%Heston is still not good as substantial density goes into zero at 
%practically realistic vol.
%Heston when density does not hit zero will work perfectly.
%I try volatility power/gamma above .65 0r .70
%gamma=.95 is proxy for lognormal, the best and coveted SV model but
%not used due to lack of tractibility.
%HEston and even smaller gamma will work far better when
%x0=1 and theta=1 as in Interest rate LMM like models with local volatility.
%Very close to zero, heston does not work very well. Whenever substantial 
%density goes into zero, model performance deteriorates.
dt=.125/32;   % Simulation time interval.%For diffusions close to zero
             %decrease dt for accuracy.
Tt=32*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/4;%Monte carlo time interval size dtM.
TtM=4*8;%Monte carlo number of simulation intervals.

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

x0=.25;   % starting value of SDE
beta1=0.0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.95;   % volatility power.
kappa=.950;   %mean reversion parameter.
theta=.10;   %mean reversion target

mu1=1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
sigma0=.90;%Volatility value
alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;
y(1:Nn)=x0;                
w(1:Nn)=x0^(1-gamma)/(1-gamma);
Z(1:Nn)=(((1:Nn)-3.5)*dNn-NnMid);


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;
%    dw0(wnStart:Nn)=(sigma0.*Z(wnStart:Nn).^1.*dtsqrt+ ...
%        dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*Z(wnStart:Nn).^1.*dtonehalf + ... %;%+ ...
%        1*d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0.*(Z(wnStart:Nn).^2-1));
    dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn) + ...
        c2(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
   dw2(wnStart:Nn)=c1(wnStart:Nn).^2.0.*Z(wnStart:Nn).^2 + ...
        c2(wnStart:Nn).^2.*(Z(wnStart:Nn).^4-2*Z(wnStart:Nn).^2+1);
   
    w(isnan(w)==1)=0;
    
    wMu0dt(isnan(wMu0dt)==1)=0;
    
    wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
 
    %Below is second order correction.Uses second hermite polynomial
    SecondOC(wnStart:Nn)=-0.027*((tt-1)*dt).* ...
        (dt)./gamma^(2*gamma-1).* ...
        sigma0^4./(w(wnStart:Nn)*(1-gamma)).*(Z(wnStart:Nn).^2-1);%sqrt(dt)changed to dt.
    
 %   plot((wnStart:Nn),SecondOC(wnStart:Nn),'g')
    
 %   str=input('Look at plots');
    
    
    w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+  ...
        SecondOC(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)));   



%Calculation of dt-integral starts    
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+2)-w(NnMidl-2)-5*D1wl*dNn-5.0*D2wl*dNn.^2.0/1.0)/12.50/dNn;
    
    wll(NnMidl-1)=w(NnMidl-2)+(.5*(D1wl+D2wl*dNn/1)+.5*(D1wl+D2wl*dNn/1+DeltaD1w))*dNn;
    wll(NnMidl)=wll(NnMidl-1)+(.5*(D1wl+D2wl*dNn/1+1.0*DeltaD1w)+.5*(D1wl+D2wl*dNn/1+2.0*DeltaD1w))*dNn;
    wll(NnMidh)=wll(NnMidl)+(.5*(D1wl+D2wl*dNn/1+2.0*DeltaD1w)+.5*(D1wl+D2wl*dNn/1+3.0*DeltaD1w))*dNn;
    wll(NnMidh+1)=wll(NnMidh)+(.5*(D1wl+D2wl*dNn/1+3.0*DeltaD1w)+.5*(D1wl+D2wl*dNn/1+4.0*DeltaD1w))*dNn;

%         
    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+2)-w(NnMidl-2)-5*D1wh*dNn-5.0*D2wh*dNn.^2.0/1.0)/12.50/dNn;
% 

    whh(NnMidh+1)=w(NnMidh+2)-(.5*(D1wh+D2wh*dNn/1)+.5*(D1wh+D2wh*dNn/1+DeltaD1w))*dNn;
    whh(NnMidh)=whh(NnMidh+1)-(.5*(D1wh+D2wh*dNn/1+1*DeltaD1w)+.5*(D1wh+D2wh*dNn/1+2.0*DeltaD1w))*dNn;
    whh(NnMidl)=whh(NnMidh)-(.5*(D1wh+D2wh*dNn/1+2.0*DeltaD1w)+.5*(D1wh+D2wh*dNn/1+3.0*DeltaD1w))*dNn;
    whh(NnMidl-1)=whh(NnMidl)-(.5*(D1wh+D2wh*dNn/1+3.0*DeltaD1w)+.5*(D1wh+D2wh*dNn/1+4.0*DeltaD1w))*dNn;

    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
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')
 paths=200000;
 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;
 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;
     YY0=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);
                 %Above is original coordinates monte carlo equation.
                     end
                 end
             end
         end
     end
 end
 YY(YY<0)=0;
 X(X<0)=0;
 
 sum(YY(:))/paths %origianl coordinates monte carlo average.
 sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
 
 BinSize=.0075/4;%Please change this bin size accordingly. When data range is too small decrease it.
 %When density range is too large increase it.Roughly, When starting point of SDE is
 %1, set at .0075*2; when starting point is x0=.1, set at .0075/4.0;
 MaxCutOff=30;
 [YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
 plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 str=input('red line is density of SDE from orthogonal incrments method, green is monte carlo.');
 
 
end
Dependencies of the code have already been uploaded in previous posts on this thread.
 
User avatar
Amin
Topic Author
Posts: 2122
Joined: July 14th, 2002, 3:00 am

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

August 28th, 2019, 6:09 pm

Post Deleted.
Last edited by Amin on September 1st, 2019, 2:50 pm, edited 1 time in total.
 
User avatar
Amin
Topic Author
Posts: 2122
Joined: July 14th, 2002, 3:00 am

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

August 29th, 2019, 5:21 am

In my previous post I talked about the second term of taylor expansion of the existing distribution variable. The total Taylor expansion approximately would go as

[$]X_0+\frac{\partial X}{\partial Z} H_1(Z)+ .5 \frac{{\partial}^2 X}{{\partial Z}^2} H_2(Z) +\frac{1}{6} \frac{{\partial}^3 X}{{\partial Z}^3} H_3(Z)+ H.O.T [$]

I believe the solution term we calculated in previous post is only an approximate proxy for [$]\frac{{\partial}^2 X}{{\partial Z}^2}[$] that takes into account one step evolution of the distribution variable and works well only because of the high relative linearity of the Lamperti transformed variable. And the These calculations would not be so straightforward for non-transformed SDE variables. I am sure once there is more research and there is a better understanding of these phenomenon, we would be evolving difficult SDEs without taking their Lamperti transform.
The first term in above taylor expansion is linear variance that I have calculated by adding increments in a squared fashion and multiplies with Z which is the first hermite polynomial. Friends would have noticed that I did not add the second term of the Taylor expansion in a squared fashion to first term and added it separately to first term because it is a different distinct term in Taylor expansion. Obviously higher order terms, if somewhere needed, would have to be added separately and not in a squared fashion to previous terms. So linear variance and higher order variance is not to be mixed up.
 
User avatar
Amin
Topic Author
Posts: 2122
Joined: July 14th, 2002, 3:00 am

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

August 31st, 2019, 3:47 pm

I am posting a new and improved version of the program. I have changed the methodology for dealing with second derivative altogether and the results are extremely good. My program posted in previous post was just alright and had some problems since long term after six or seven years became very poor and there were also suboptimalities in between for smaller maturities. This program is simply excellent and I have tried it for ten years in the code I am posting with superb results. It is a great improvement after I made changes to my methodology. Here is the code. I will explain the changes and new methodology further with equations in next post in just a few hours. For ten years, for example the difference in mean of the density in the code posted is just three basis points with quite high volatility and mean reversion. Please enjoy the new program. 

function [] = ItoTaylorDensityMRRecombining14F4Wilmott()
%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)
%where mu1=kappa*theta and mu2=-kappa with
%beta1=0 and beta2=1;
%I have not directly simulated the SDE but simulated the transformed 
%Bessel 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.
%The present program is meant for mean-reversion SDEs that take two terms
%in non-linear drift. It will work for simpler SDEs but for that you will 
%have to turn SecondOC variable, second order correction variable to zero.
%
%For diffusions going into zero, you may see instabilities, I have not
%corrected for that yet in this initial version.
%Heston is still not good as substantial density goes into zero at 
%practically realistic vol.
%Heston when density does not hit zero will work perfectly.
%I try volatility power/gamma above .65 0r .70
%gamma=.95 is proxy for lognormal, the best and coveted SV model but
%not used due to lack of tractibility.
%HEston and even smaller gamma will work far better when
%x0=1 and theta=1 as in Interest rate LMM like models with local volatility.
%Very close to zero, heston does not work very well. Whenever substantial 
%density goes into zero, model performance deteriorates.
dt=.125/32;   % Simulation time interval.%For diffusions close to zero
             %decrease dt for accuracy.
Tt=32*80;     % 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/4;%Monte carlo time interval size dtM.
TtM=4*80;%Monte carlo number of simulation intervals.

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

x0=.25;   % starting value of SDE
beta1=0.0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.95;   % volatility power.
kappa=.90;   %mean reversion parameter.
theta=.10;   %mean reversion target

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;
y(1:Nn)=x0;                
w(1:Nn)=x0^(1-gamma)/(1-gamma);
Z(1:Nn)=(((1:Nn)-3.5)*dNn-NnMid);
H4(1:Nn)=Z(1:Nn).^4-6*Z(1:Nn)+3;

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;
SecondOC(1:Nn)=0;
SecondD(1:Nn)=0.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;
    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;
%    dw0(wnStart:Nn)=(sigma0.*Z(wnStart:Nn).^1.*dtsqrt+ ...
%        dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*Z(wnStart:Nn).^1.*dtonehalf + ... %;%+ ...
%        1*d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0.*(Z(wnStart:Nn).^2-1));
    dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn);% + ...
    d2w(wnStart:Nn)=c2(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
   dw2(wnStart:Nn)=c1(wnStart:Nn).^2.0.*Z(wnStart:Nn).^2;% + ...
   d2w2(wnStart:Nn)= c2(wnStart:Nn).^2.*(Z(wnStart:Nn).^4-2*Z(wnStart:Nn).^2+1);
   
    w(isnan(w)==1)=0;
    
    wMu0dt(isnan(wMu0dt)==1)=0;
    
    wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
    w2MeanPrev=sum(ZProb(wnStart:Nn).*(w(wnStart:Nn)-wMeanPrev).^2);
    w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)-  ...
        1/sqrt(2*pi)*.5*(sign((w(wnStart:Nn)-wMeanPrev).^2-w2MeanPrev+d2w(wnStart:Nn))).*...
        sqrt(abs(sign((w(wnStart:Nn)-wMeanPrev).^2-w2MeanPrev).*((w(wnStart:Nn)-wMeanPrev).^2-w2MeanPrev).^2+ ...
        sign(d2w(wnStart:Nn)).*d2w2(wnStart:Nn))).*sigma0^2*dt/2+ ...
        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)));  
   

%Calculation of dt-integral starts    
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+2)-w(NnMidl-2)-5*D1wl*dNn-5.0*D2wl*dNn.^2.0/1.0)/12.50/dNn;
    
    wll(NnMidl-1)=w(NnMidl-2)+(.5*(D1wl+D2wl*dNn/1)+.5*(D1wl+D2wl*dNn/1+DeltaD1w))*dNn;
    wll(NnMidl)=wll(NnMidl-1)+(.5*(D1wl+D2wl*dNn/1+1.0*DeltaD1w)+.5*(D1wl+D2wl*dNn/1+2.0*DeltaD1w))*dNn;
    wll(NnMidh)=wll(NnMidl)+(.5*(D1wl+D2wl*dNn/1+2.0*DeltaD1w)+.5*(D1wl+D2wl*dNn/1+3.0*DeltaD1w))*dNn;
    wll(NnMidh+1)=wll(NnMidh)+(.5*(D1wl+D2wl*dNn/1+3.0*DeltaD1w)+.5*(D1wl+D2wl*dNn/1+4.0*DeltaD1w))*dNn;

%         
    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+2)-w(NnMidl-2)-5*D1wh*dNn-5.0*D2wh*dNn.^2.0/1.0)/12.50/dNn;
% 

    whh(NnMidh+1)=w(NnMidh+2)-(.5*(D1wh+D2wh*dNn/1)+.5*(D1wh+D2wh*dNn/1+DeltaD1w))*dNn;
    whh(NnMidh)=whh(NnMidh+1)-(.5*(D1wh+D2wh*dNn/1+1*DeltaD1w)+.5*(D1wh+D2wh*dNn/1+2.0*DeltaD1w))*dNn;
    whh(NnMidl)=whh(NnMidh)-(.5*(D1wh+D2wh*dNn/1+2.0*DeltaD1w)+.5*(D1wh+D2wh*dNn/1+3.0*DeltaD1w))*dNn;
    whh(NnMidl-1)=whh(NnMidl)-(.5*(D1wh+D2wh*dNn/1+3.0*DeltaD1w)+.5*(D1wh+D2wh*dNn/1+4.0*DeltaD1w))*dNn;

    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
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')
 paths=200000;
 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;
 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;
     YY0=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);
                 %Above is original coordinates monte carlo equation.
                     end
                 end
             end
         end
     end
 end
 YY(YY<0)=0;
 X(X<0)=0;
 
 sum(YY(:))/paths %origianl coordinates monte carlo average.
 sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
 
 BinSize=.0075/4;%Please change this bin size accordingly. When data range is too small decrease it.
 %When density range is too large increase it.Roughly, When starting point of SDE is
 %1, set at .0075*2; when starting point is x0=.1, set at .0075/4.0;
 MaxCutOff=30;
 [YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
 plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 str=input('red line is density of SDE from orthogonal incrments method, green is monte carlo.');
 
 
end
The above program was tried for volatility exponents upwards of .65. It will not work for heston due to other technical reasons in the program and I will fix them in a few days and assure friends that same "main algorithm" will work perfectly fine for heston as well. 
 
User avatar
Amin
Topic Author
Posts: 2122
Joined: July 14th, 2002, 3:00 am

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

September 1st, 2019, 2:34 pm

Since my work is very close to completion, I would like to explain my method in a more formal way.
Let us consider a non-linear SDE of the form 
[$]dX(t)=\mu_1 X^{\beta_1}(t) dt +\mu_2 X^{\beta_2}(t) dt + \sigma X^{\gamma} dz(t)[$]
In my work I took a Lamperti transform of the above SDE and defined the new transformed variable as [$]w(t)=\frac{X(t)^{(1-\gamma)}}{(1-\gamma)}[$]. I did not take [$]\sigma[$] in the denominator of the transformed variable.
Applying Ito change of variable formula on the transformed variable, we get
[$]dw(t)=d[\frac{X(t)^{(1-\gamma)}}{1-\gamma}]=\mu_1 X(t)^{(\beta_1-\gamma)} dt+ \mu_2 X(t)^{(\beta_2-\gamma)} dt[$]
[$] -.5 \gamma {\sigma}^2 X(t)^{(\gamma-1)} dt +\sigma dz(t)[$] 
Integral form of the above SDE would be written as
[$]w(t)=\frac{X(t)^{(1-\gamma)}}{1-\gamma}=w(0)+\int_0^t \mu_1 X(s)^{(\beta_1-\gamma)} ds+\int_0^t \mu_2 X(s)^{(\beta_2-\gamma)} ds[$]
[$] -\int_0^t .5 \gamma {\sigma}^2 X(s)^{(\gamma-1)} ds +\int_0^t \sigma dz(s)[$] 

To get a second order Ito-Taylor expansion, we will do a repeated Ito on each of the variable drift term of the above transformed SDE of [$]w(t)[$]
Taking the repeated Ito of the first term in the SDE, we get
[$]\int_0^t \mu_1 X(s)^{(\beta_1-\gamma)} ds= \mu_1 X(0)^{(\beta_1-\gamma)} \int_0^t ds[$]
[$] + \int_0^t \int_0^s \mu_1(\beta_1-\gamma) X(v)^{(\beta_1-\gamma-1)} (\mu_1 X(v)^{\beta_1} + \mu_2 X(v)^{\beta_2})dv ds[$]
[$]+\int_0^t \int_0^s \mu_1(\beta_1-\gamma) X(v)^{(\beta_1-\gamma-1)} \sigma X(v)^{\gamma} dz(v) ds[$]
[$]+\int_0^t \int_0^s .5 \mu_1  (\beta_1-\gamma)(\beta_1-\gamma-1) X(v)^{(\beta_1-\gamma-2)} {\sigma}^2 X(v)^{2 \gamma} dv ds[$]
notice that first term in the above expansion is calculated at initial time zero and rest of the terms are calculated along the evolution path of the stochastic integrals. We could continue higher order expansion inside the integrals but if the interval is small we can just stop the expansion and find a very good approximation by substituting the terms that have to be calculated along the evolution path by their values at time zero. Our above expansion would become
[$]\int_0^t \mu_1 X(s)^{(\beta_1-\gamma)} ds]= \mu_1 X(0)^{(\beta_1-\gamma)} \int_0^t ds[$]
[$] + \mu_1(\beta_1-\gamma) X(0)^{(\beta_1-\gamma-1)} (\mu_1 X(0)^{\beta_1} + \mu_2 X(0)^{\beta_2}) \int_0^t \int_0^s dv ds[$]
[$]+ \mu_1(\beta_1-\gamma) X(0)^{(\beta_1-\gamma-1)} \sigma X(0)^{\gamma} \int_0^t \int_0^s  dz(v) ds[$]
[$]+ .5 \mu_1(\beta_1-\gamma)(\beta_1-\gamma-1) X(0)^{(\beta_1-\gamma-2)} {\sigma}^2 X(0)^{2 \gamma} \int_0^t \int_0^s  dv ds[$]
We know from basic calculus and stochastic calculus that 
[$]\int_0^t \int_0^s dv ds = \frac{t^2}{2} [$] and
[$]\int_0^t \int_0^s dz(v) ds = (1-\frac{1}{\sqrt{3}}) t^{1.5} Z [$] where [$]Z[$] is a standard normal.

We have to find repeated Ito of all the three drift terms. There is no further expansion of the volatility terms as it has no variable coefficients. After the Ito Taylor expansion of the all three drift terms, the expansion of the transformed SDE becomes
[$]w(t)=\frac{X(t)^{(1-\gamma)}}{1-\gamma}=w(0)+ \mu_1 X(s)^{(\beta_1-\gamma)} {\Delta t}+ \mu_2 X(s)^{(\beta_2-\gamma)} {\Delta t}[$]
[$] - .5 \gamma {\sigma}^2 X(s)^{(\gamma-1)} {\Delta t} [$]
[$] + (\mu_1(\beta_1-\gamma) X(0)^{(\beta_1-\gamma-1)} +\mu_2(\beta_2-\gamma) X(0)^{(\beta_2-\gamma-1)} -.5\gamma (\gamma-1) X(0)^{(\gamma-2)})[$]
[$]* (\mu_1 X(0)^{\beta_1} + \mu_2 X(0)^{\beta_2})  \frac{{\Delta t}^2}{2}[$]
[$]+ .5 (\mu_1(\beta_1-\gamma)(\beta_1-\gamma-1) X(0)^{(\beta_1-\gamma-2)} +\mu_2(\beta_2-\gamma)(\beta_2-\gamma-1) X(0)^{(\beta_2-\gamma-2)} -.5\gamma (\gamma-1)(\gamma-2) X(0)^{(\gamma-3)} )[$]
[$]*{\sigma}^2 X(0)^{2 \gamma} \frac{{\Delta t}^2}{2}[$]
[$]+ \sigma \sqrt{\Delta t} Z[$] 
[$]+ ( \mu_1(\beta_1-\gamma) X(0)^{(\beta_1-\gamma-1)} +  \mu_2(\beta_2-\gamma) X(0)^{(\beta_2-\gamma-1)}-.5\gamma (\gamma-1) X(0)^{(\gamma-2)}) (\sigma X(0)^{\gamma} ) (1-\frac{1}{\sqrt{3}}) {\Delta t}^{1.5} Z [$]

We can write the one step differential evolution equation as 
[$]w(t+\Delta t)= w(t)+ \mu_{\Delta t} + \alpha_{\Delta t} Z [$]
where particularly 
[$]\alpha_{\Delta t}= \sigma \sqrt{\Delta t}+ ( \mu_1(\beta_1-\gamma) X(0)^{(\beta_1-\gamma-1)} +  \mu_2(\beta_2-\gamma) X(0)^{(\beta_2-\gamma-1)}-.5\gamma (\gamma-1) X(0)^{(\gamma-2)}) (\sigma X(0)^{\gamma} ) (1-\frac{1}{\sqrt{3}}) {\Delta t}^{1.5})[$]
After looking at Ito-Taylor expansion, we proceed towards the main density evolution algorithm that advances the density over each small time interval in our setting that is different from monte carlo simulation.
I have divided the density into N subdivisions with a value of normal random variable linearly associated with each subdivision. In my program, these subdivisions are linearly associated with normal RV. Let [$]Z_n[$] be value of normal RV associated with nth subdivision with roughly [$]Z_1=-4.0[$] and [$]Z_N=4.0[$] . In my current program I have taken [$]N=46[$] covering the gaussian from -4.5 to +4.5.
I took Lamperti transform of the SDE in original coordinates X and the transformed SDE variable w became very linear with respect to normal RV which greatly helped in analytics. In fact, in most of the SDEs only first order variance was of relevance. 
I define existing density as starting transformed density at any point [$]t[$] in evolution and it is written as [$]w(t)[$] , and then a differential term is added for an evolution period [$]\Delta t[$] and then we find density at next step [$]T=t+\Delta t[$] denoted as [$]w(T)[$]
In order to calculate the one step evolution of first order variance, we have to calculate the first order Hermite representation of the existing SDE along nth subdivision or nth point in the grid and also calculate the one step first order hermite representation of the innovation/differential term. It is very easy to calculate the hermite represntation of the differential term due to Ito-Taylor expansion of the SDE and there is a huge literature on that. With help of Ito-Taylor expansion we can easily calculate the evolution in terms of hermite polynomials and thereby decompose the variance into first order and higher orders based on association with each hermite polynomial of different order. After taking Lamperti transform, it was usually the first order variance that was significant. But despite taking Lamperti transform, on each step SDE still evolves slightly non linearly due to non-linearity in drift terms and it was not straightforward how to represent the existing SDE variable [$]w(t)[$] in terms of hermite polynomials. If we could represent the existing and differential parts of the SDE variable in terms of hermite polynomials, we could get a representation of the evolved variable [$]w_T(n)[$] in terms of hermite polynomial.
Hermite representation of existing SDE variable on each point of the evolving grid could be represented as [$]w_t(n)=\mu_w+\alpha_t(n)Z(n)[$] and hermite representation of the differential term could be written as [$]w_{\Delta t}(n)=\mu_{\Delta t}(n)+ \alpha_{\Delta t}(n) Z(n)[$], then representation of evolved SDE at next step [$]T=t+\Delta t[$]  i.e [$]w_T[$] could be written as [$]w_T(n)= \mu_{\Delta t}(n)+\mu_w+ \sqrt{{\alpha}^2_t(n)+{\alpha}^2_{\Delta t}(n)} Z(n)[$]
where [$]\mu_w[$] associated with existing density w(t) is defined as [$]\mu_w = \int_{-inf}^{+inf} w p(w) dw[$]
Here [$]\mu_{\Delta t}(n)[$] is the total deterministic contribution from all the drift terms in the SDE for the nth subdivision for time interval [$]\Delta t[$].
Again it was very easy to find a Hermite representation of the differential part of SDE due to Ito-Taylor but it was not simple at all to find a hermite representation of the existing SDE on each step. My contribution came in here as I noticed that [$]w(t)[$] was very close to scaled Z and it could be possible to omit the step of finding the hermite representation of the existing SDE and instead use the first order evolution equation for the SDE as
[$]w_T(n)=\mu_{\Delta t}(n)+\mu_{w} + \sqrt{(w_t(n)-\mu_{w})^2+(\alpha_{\Delta t}(n) Z(n))^2}[$]
This worked since [$]w_t[$] is reasonably linear with respect to Z due to Lamperti transform and [$]w_t(n)-\mu_w[$] was a valid proxy for [$]\alpha_t(n) Z(n)[$].
When I made this first step approximation, things worked perfectly with linear SDEs where [$]w(t)[$] remains linear with [$]Z[$]. But particularly for mean reverting type SDEs which have multiple terms in the drift and when volatility is high, due to non-linear dynamics in the drift, a huge quantity of second order derivative of [$]w(t)[$] with respect to [$]Z[$] is introduced and dealing with the first order is not enough to get good results for the density of mean reverting type SDEs.
Then a recent breakthrough came when I found from my experiments that second order hermite term in the expansion of the existing SDE [$]w(t)[$] was making a contribution to the SDE evolution. If we denote the hermite representation of [$]w(t)[$] up to second order as [$]w_t(n)=\alpha_t(n) Z(n) + \beta_t(n) (Z^2(n)-1)[$].
This contribution from the coefficients of second hermite polynomial to one step evolution of the SDE was empirically found up to a positive scaling constant [$]c[$]  as equal to [$] -.5 c \beta_t(n)(Z^2(n)-1) {\sigma}^2(t) \frac{\Delta t}{2}[$]
Here [$]\sigma(t)[$] is volatility in the noise term of the SDE. 
But again all of this was approximate since it is not very simple to find a hermite representation of the existing SDE and it was particularly difficult to find coefficients of second hermite polynomial accurately. So analogously to my approximation for coefficients of coefficients of first order hermite polynomials for the existing density, I devised a similar approximation to coefficients of second order hermite polynomial and used the formula
[$]\beta_t(n) (Z^2(n)-1) = (w_t(n)-\mu_w)^2- \sigma^2_w[$] which is very valid due to high relative linearity of the transformed variable.
So our contribution from the second order derivative from the surface of w(t) could now be written as
[$]-.5 c \beta_t(n)(Z^2(n)-1) {\sigma}^2(t) \frac{\Delta t}{2} = -.5 c ((w_t(n)-\mu_w)^2- \sigma^2_w)  {\sigma}^2(t) \frac{\Delta t}{2}[$]
Here [$]\sigma^2_w=\int_{-\inf}^{+\inf} (w-\mu_w)^2 p(w) dw[$]
if there is a contribution from the second order hermite polynomial of differential part of SDE evolution, it would be added to the above second order contribution from the existing surface w(t) in a squared fashion and resulting second order addition to one step evolution of the SDE would become as
[$]\sqrt{(.5 c ((w_t(n)-\mu_w)^2- \sigma^2_w)  {\sigma}^2(t) \frac{\Delta t}{2})^2 + (\beta_{\Delta t}(n) (Z(n)^2-1))^2}[$]
where [$] \beta_{\Delta t} [$] are the coefficients of second order hermite polynomial in the Ito-Taylor expansion of the differential term of the SDE. After the second order correction, the evolution equation for the SDE becomes
[$]w_T(n)=\mu_{\Delta t}(n)+\mu_{w} + \sqrt{(w_t(n)-\mu_{w})^2+(\alpha_{\Delta t}(n) Z(n))^2}[$]
[$]+ \sqrt{(.5 c ((w_t(n)-\mu_w)^2- \sigma^2_w)  {\sigma}^2(t) \frac{\Delta t}{2})^2 + (\beta_{\Delta t}(n) (Z(n)^2-1))^2}[$]

In our simulations we found that second order contribution from the differential term after taking Lamperti transform was completely miniscule as compared to the second order contribution from the existing density surface and this contribution from differential term to the density average never made a difference of more than one basis point (.0001) in our simulations with high volatility and mean reversion. We write this as
 [$]\beta_{\Delta t}(n) (Z(n)^2-1) << (.5 c ((w_t(n)-\mu_w)^2- \sigma^2_w)  {\sigma}^2(t) \frac{\Delta t}{2})[$]
So we could safely write the evolution equation when we are simulating density of a Lamperti transformed variable as
[$]w_T(n)=\mu_{\Delta t}(n)+\mu_{w} + \sqrt{(w_t(n)-\mu_{w})^2+(\alpha_{\Delta t}(n) Z(n))^2}[$]
[$]- .5 c ((w_t(n)-\mu_w)^2- \sigma^2_w)  {\sigma}^2(t) \frac{\Delta t}{2}[$]
Since the drift term is highly non-linear in cases of high volatility and mean-reversion of SDEs and it introduces higher derivatives despite the linearity of diffusion term when evolving a Lamperti transformed SDE, there is a possibility that there could be a slight amount of third order derivative contribution from the existing density [$]w(t)[$] surface when evolving the density over a small time interval. The form of this contribution can easily be guessed from the form of second order contribution and we could easily write third order hermite polynomial representation of the existing density analogously to what we have written for first order and second order as
[$] \gamma_t(n) (Z(n)^3-3Z(n))=(w_t(n)-\mu_w)^3- 3 \sigma^2_w  (w_t(n)-\mu_w)[$]  here [$]\gamma_t(n)[$] is the coefficient of third hermite polynomial in expansion of existing SDE surface w(t).
Then the total contribution to one step evolution from the existing surface [$]w(t)[$] becoming equivalent to the following term as
 [$]\frac{1}{6} c_2 ((w_t(n)-\mu_w)^3- 3 \sigma^2_w (w_t(n)-\mu_w))  {\sigma}^3(t) \frac{{\Delta t}^{1.5}}{6}[$]

Please notice that you have to account for signs in all the squared additions in the above post. Please look at the program in post 802 for the precise equations.
Last edited by Amin on September 2nd, 2019, 10:31 am, edited 26 times in total.
 
User avatar
Amin
Topic Author
Posts: 2122
Joined: July 14th, 2002, 3:00 am

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

September 1st, 2019, 2:46 pm

Here is a slightly modified program from previous post 800. It has corrected a mistake that was not making a difference in the evolution of transformed SDE but could make a difference if one tried to simulate SDE in original variables
.
function [] = ItoTaylorDensityMRRecombining14F4Wilmott()
%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)
%where mu1=kappa*theta and mu2=-kappa with
%beta1=0 and beta2=1;
%I have not directly simulated the SDE but simulated the transformed 
%Bessel 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.
%The present program is meant for mean-reversion SDEs that take two terms
%in non-linear drift. It will work for simpler SDEs but for that you will 
%have to turn SecondOC variable, second order correction variable to zero.
%
%For diffusions going into zero, you may see instabilities, I have not
%corrected for that yet in this initial version.
%Heston is still not good as substantial density goes into zero at 
%practically realistic vol.
%Heston when density does not hit zero will work perfectly.
%I try volatility power/gamma above .65 0r .70
%gamma=.95 is proxy for lognormal, the best and coveted SV model but
%not used due to lack of tractibility.
%HEston and even smaller gamma will work far better when
%x0=1 and theta=1 as in Interest rate LMM like models with local volatility.
%Very close to zero, heston does not work very well. Whenever substantial 
%density goes into zero, model performance deteriorates.
dt=.125/32;   % Simulation time interval.%For diffusions close to zero
             %decrease dt for accuracy.
Tt=32*80;     % 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/4;%Monte carlo time interval size dtM.
TtM=4*80;%Monte carlo number of simulation intervals.

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

x0=.250;   % starting value of SDE
beta1=0.0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.95;   % volatility power.
kappa=.950;   %mean reversion parameter.
theta=.10;   %mean reversion target

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;
y(1:Nn)=x0;                
w(1:Nn)=x0^(1-gamma)/(1-gamma);
Z(1:Nn)=(((1:Nn)-3.5)*dNn-NnMid);
H4(1:Nn)=Z(1:Nn).^4-6*Z(1:Nn)+3;

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;
SecondOC(1:Nn)=0;
SecondD(1:Nn)=0.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;
    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;
%    dw0(wnStart:Nn)=(sigma0.*Z(wnStart:Nn).^1.*dtsqrt+ ...
%        dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*Z(wnStart:Nn).^1.*dtonehalf + ... %;%+ ...
%        1*d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0.*(Z(wnStart:Nn).^2-1));
    dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn);% + ...
    d2w(wnStart:Nn)=c2(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
   dw2(wnStart:Nn)=c1(wnStart:Nn).^2.0.*Z(wnStart:Nn).^2;% + ...
   d2w2(wnStart:Nn)= c2(wnStart:Nn).^2.*(Z(wnStart:Nn).^4-2*Z(wnStart:Nn).^2+1);
   
   
   
   
    w(isnan(w)==1)=0;
    
    wMu0dt(isnan(wMu0dt)==1)=0;
    
    wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
    w2MeanPrev=sum(ZProb(wnStart:Nn).*(w(wnStart:Nn)-wMeanPrev).^2);
    SOT(wnStart:Nn)=-1/sqrt(2*pi)*.5*((w(wnStart:Nn)-wMeanPrev).^2-w2MeanPrev).*sigma0^2*dt/2;
    w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+  ...
        (sign(SOT(wnStart:Nn)+d2w(wnStart:Nn))).*...
        sqrt(abs(sign(SOT(wnStart:Nn)).*SOT(wnStart:Nn).^2+ ...
        sign(d2w(wnStart:Nn)).*d2w2(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)));  
   

%Calculation of dt-integral starts    
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+2)-w(NnMidl-2)-5*D1wl*dNn-5.0*D2wl*dNn.^2.0/1.0)/12.50/dNn;
    
    wll(NnMidl-1)=w(NnMidl-2)+(.5*(D1wl+D2wl*dNn/1)+.5*(D1wl+D2wl*dNn/1+DeltaD1w))*dNn;
    wll(NnMidl)=wll(NnMidl-1)+(.5*(D1wl+D2wl*dNn/1+1.0*DeltaD1w)+.5*(D1wl+D2wl*dNn/1+2.0*DeltaD1w))*dNn;
    wll(NnMidh)=wll(NnMidl)+(.5*(D1wl+D2wl*dNn/1+2.0*DeltaD1w)+.5*(D1wl+D2wl*dNn/1+3.0*DeltaD1w))*dNn;
    wll(NnMidh+1)=wll(NnMidh)+(.5*(D1wl+D2wl*dNn/1+3.0*DeltaD1w)+.5*(D1wl+D2wl*dNn/1+4.0*DeltaD1w))*dNn;

%         
    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+2)-w(NnMidl-2)-5*D1wh*dNn-5.0*D2wh*dNn.^2.0/1.0)/12.50/dNn;
% 

    whh(NnMidh+1)=w(NnMidh+2)-(.5*(D1wh+D2wh*dNn/1)+.5*(D1wh+D2wh*dNn/1+DeltaD1w))*dNn;
    whh(NnMidh)=whh(NnMidh+1)-(.5*(D1wh+D2wh*dNn/1+1*DeltaD1w)+.5*(D1wh+D2wh*dNn/1+2.0*DeltaD1w))*dNn;
    whh(NnMidl)=whh(NnMidh)-(.5*(D1wh+D2wh*dNn/1+2.0*DeltaD1w)+.5*(D1wh+D2wh*dNn/1+3.0*DeltaD1w))*dNn;
    whh(NnMidl-1)=whh(NnMidl)-(.5*(D1wh+D2wh*dNn/1+3.0*DeltaD1w)+.5*(D1wh+D2wh*dNn/1+4.0*DeltaD1w))*dNn;

    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
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')
 paths=200000;
 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;
 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;
     YY0=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);
                 %Above is original coordinates monte carlo equation.
                     end
                 end
             end
         end
     end
 end
 YY(YY<0)=0;
 X(X<0)=0;
 
 sum(YY(:))/paths %origianl coordinates monte carlo average.
 sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
 
 BinSize=.0075/4;%Please change this bin size accordingly. When data range is too small decrease it.
 %When density range is too large increase it.Roughly, When starting point of SDE is
 %1, set at .0075*2; when starting point is x0=.1, set at .0075/4.0;
 MaxCutOff=30;
 [YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
 plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 str=input('red line is density of SDE from orthogonal incrments method, green is monte carlo.');
 
 
end
 
User avatar
Amin
Topic Author
Posts: 2122
Joined: July 14th, 2002, 3:00 am

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

September 6th, 2019, 5:26 pm

For the friends, who have been following the research, there are two new developments that I want to present.

1. Density Moments Defect correcting Hermite Polynomials.
I worked with this idea several years ago for the mean of the distribution but generalized it to higher moments just a few days ago. I have still not worked with higher moments but this is a very good idea. When the density shape is very close to the true shape and moments are also quite close to the true value, moment defect correcting Hermite polynomials can be a wonderful approach. I will explain the idea with the mean first and then generalize to higher moments.
Suppose we know the true mean[$]\mu_1[$] of a density X, and we have a density grid on N points with each point represented by [$]X(n)[$] .We multiply the density with an expression of the form [$]a(X(n)-b)[$] and place the restrictions that equation for density sums to one and equation for mean sums to [$]\mu_1[$] .This can be written in equations. Since I used summations in my program, it can be written as
Equation 1.   [$]\sum_1^{N} a(X(n)-b) P(X(n))=1.0[$] where P(X(n)) is integrated probability mass in nth subdivision.
Equation 2.  [$]\sum_1^{N} a(X(n)-b) X(n) P(X(n))=\mu_1[$]
So we have two equations and two unknowns and we could very easily solve for the unknowns 'a' and 'b' . Even first order correction works wonders since mean is corrected under the constraint that density sums to one. In my implementation, I did not change the density but updated the density variable as
[$]X_{Corrected}(n)=a(X(n)-b)X(n)[$]. First Mean correcting polynomial is usually a line with miniscule slope increasing the density on one side of b and decreasing it on the other side of b with the constraint that density still sums to one. If the defect is small, density remains positive and b is chosen outside of the support of the density. But it would be better to add 1.0 to the correction polynomial as I have done in generalized example below. After I got the density scaling coefficients, I did not change the density and kept the density on the original Z grid and only scaled X as its scaling has exactly the same meanings as changing the amplitude of density.  
Generalizing to higher order moment defect correcting hermite polynomimals, we have after including a second order hermite polynomial that corrects for variance, we get the following three equations which can be algebraically solved for three unknowns as
Equation 1. [$]\sum_1^N (1.0+\sigma_c (X(n)-\mu_c) + \kappa_c( (X(n)-\mu_c)^2 - {\sigma}^2_c)) P(X_n)=1[$]
Equation 2.[$]\sum_1^N (1.0+\sigma_c (X(n)-\mu_c) + \kappa_c( (X(n)-\mu_c)^2 - {\sigma}^2_c)) X(n) P(X_n)=\mu_1[$]
Equation 2.[$]\sum_1^N (1.0+\sigma_c (X(n)-\mu_c) + \kappa_c( (X(n)-\mu_c)^2 - {\sigma}^2_c)) X^2(n) P(X_n)=\mu_2[$] 
Here [$]\mu_c[$], [$]\sigma_c[$] and [$]\kappa_c[$], the coefficients with subscript c are respectively mean, and coefficients of first and second mean correcting hermite polynomials. And corrected value of X(n) would be
[$]X_{Corrected}(n)=  (1.0+\sigma_c (X(n)-\mu_c) + \kappa_c( (X(n)-\mu_c)^2 - {\sigma}^2_c)) X(n)[$]
If we know the moments upto fourth order and we can easily solve the resulting five equations algebraically, we can write the equations for 
mean correction in terms of four hermite polynomijals and make sure first four moments are good. I have done this with only the mean yet. I will post a program with first order defect correction with graphs on the weekend.

2. Density Construction Algorithm Improvement.
I have been able to get analytic density shapes very close to the true density. Previously analytic density was usually more peaked in the middle which remained problematic but this problem has decreased to a high degree. I have been able to get far better density shapes by slightly changing the mean I have used in the one step density simulation algorithm. Usually mean is higher than median in the Lamperti transformed coordinates but the difference is very small. The difference is almost zero when second derivative in the density is not significant but the difference is slightly larger for mean reverting type densities which have a significant second derivative with respect to normal density. setting the mean at a slightly higher value helps make the shape of density more exact. This new mean can be written as 
new mean = true mean + .25*( true mean - median)
This addition to true mean is very small but helps bring the shape of the density closer to true shape especially where we saw a larger peak, the peak is closer to true value.
I will be putting more worked out code with these two above mentioned improvements over the weekend. I will also be showing a lot of graph comparisons with true density for a range of values of parameters for the SDEs.
Last edited by Amin on September 8th, 2019, 9:41 am, edited 1 time in total.
 
User avatar
Amin
Topic Author
Posts: 2122
Joined: July 14th, 2002, 3:00 am

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

September 8th, 2019, 7:38 am

I am writing this post to mention that I have added a 1.0 to moment defect correcting hermite polynomial equations. If you looked at the equations earlier, please note that I have made a slight change and edited the equations to add a one.
 
User avatar
Amin
Topic Author
Posts: 2122
Joined: July 14th, 2002, 3:00 am

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

Yesterday, 6:23 pm

Here are a few graphs showing a comparison between monte carlo density and Ito-Hermite method density. Here the SDE is 

X(t+1)=X(t)+kappa(theta-X(t))*dt +sigma * X(t)^gamma *dz(t);   with X(0)=x0; Graphs are self-explanatory. In the next post, I will give the code used to make these graphs.
The graphs for T>=2 years are slightly off at the peak of the density. This is because I chose some constants by trial and error and did not have enough time today to perfect the values to complete precision.
First moment is matched by moment matching technique as I had mentioned in a previous post.



Image



Image



Image
Image



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

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

Yesterday, 6:32 pm

I realized that second hermite polynomial related to existing density not only requires a volatility adjustment, it also requires a drift adjustment and I had far better results after including a drift adjustment for second hermite. I will post another updated program in a few days with updated coefficients as the current coefficients were chosen only with a rough trial and error and needed more work to be determined to perfect precision. I showed the graphs for x0=.25 and theta =.1. but for interest rate derivatives, the results were more impressive when x0=1 and theta =1 as usually taken in interest rate derivatives. I will be showing those graphs later. Here is the new code.
function [] = ItoHermiteMethodWilmott()
%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)
%where mu1=kappa*theta and mu2=-kappa with
%beta1=0 and beta2=1;
%I have not directly simulated the SDE but simulated the transformed 
%Bessel 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.
%The present program is meant for mean-reversion SDEs that take two terms
%in non-linear drift. It will work for simpler SDEs but for that you will 
%have to turn SecondOC variable, second order correction variable to zero.
%
%For diffusions going into zero, you may see instabilities, I have not
%corrected for that yet in this initial version.
%Heston is still not good as substantial density goes into zero at 
%practically realistic vol.
%Heston when density does not hit zero will work perfectly.
%I try volatility power/gamma above .65 0r .70
%gamma=.95 is proxy for lognormal, the best and coveted SV model but
%not used due to lack of tractibility.
%HEston and even smaller gamma will work far better when
%x0=1 and theta=1 as in Interest rate LMM like models with local volatility.
%Very close to zero, heston does not work very well. Whenever substantial 
%density goes into zero, model performance deteriorates.
dt=.125/8;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.
Tt=8*2;     % 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/4;%Monte carlo time interval size dtM.
TtM=4*2;%Monte carlo number of simulation intervals.



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

x0=.250;   % starting value of SDE
beta1=0.0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.95;   % volatility power.
kappa=1.0;   %mean reversion parameter.
theta=.10;   %mean reversion target
%MonteTp5dtp03125x0p1Thp1Ka1Gap75Sip8

mu1=1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
sigma0=.90;%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)-3.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;
%    dw0(wnStart:Nn)=(sigma0.*Z(wnStart:Nn).^1.*dtsqrt+ ...
%        dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*Z(wnStart:Nn).^1.*dtonehalf + ... %;%+ ...
%        1*d_dGG_VV_dyy(wnStart:Nn).*VV(wnStart:Nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0.*(Z(wnStart:Nn).^2-1));
    dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn);% + ...
   dw2(wnStart:Nn)=c1(wnStart:Nn).^2.0.*Z(wnStart:Nn).^2;        
%    d2w(wnStart:Nn)=c2(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
%   d2w2(wnStart:Nn)= c2(wnStart:Nn).^2.*(Z(wnStart:Nn).^4-2*Z(wnStart:Nn).^2+1);
   
    w(isnan(w)==1)=0;
  
    wMu0dt(isnan(wMu0dt)==1)=0;
    
    wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.
    w2MeanPrev=sum(ZProb(wnStart:Nn).*(w(wnStart:Nn)-wMeanPrev).^2);
    %SOT_dz2(wnStart:Nn)=-(1/sqrt(2*pi))*gamma*((w(wnStart:Nn)-wMeanPrev).^2-w2MeanPrev).*c1(wnStart:Nn).^2/2;
    %SOT_mudt(wnStart:Nn)=-(1/sqrt(2*pi))*.25*((w(wnStart:Nn)-wMeanPrev).^2-w2MeanPrev).*wMu0dt(wnStart:Nn);
    
    SOT_dz2(wnStart:Nn)=-.5*gamma^2*((w(wnStart:Nn)-wMeanPrev).^2-w2MeanPrev).*c1(wnStart:Nn).^2/2;
    SOT_mudt(wnStart:Nn)=-(1/sqrt(2*pi))*.33*((w(wnStart:Nn)-wMeanPrev).^2-w2MeanPrev).*wMu0dt(wnStart:Nn);
     
    %wPrev=w;
    w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+  ...
        SOT_dz2(wnStart:Nn)+ ...
        SOT_mudt(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)));  


%Calculation of dt-integral starts    
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+2)-w(NnMidl-2)-5*D1wl*dNn-5.0*D2wl*dNn.^2.0/1.0)/12.50/dNn;
    
    wll(NnMidl-1)=w(NnMidl-2)+(.5*(D1wl+D2wl*dNn/1)+.5*(D1wl+D2wl*dNn/1+DeltaD1w))*dNn;
    wll(NnMidl)=wll(NnMidl-1)+(.5*(D1wl+D2wl*dNn/1+1.0*DeltaD1w)+.5*(D1wl+D2wl*dNn/1+2.0*DeltaD1w))*dNn;
    wll(NnMidh)=wll(NnMidl)+(.5*(D1wl+D2wl*dNn/1+2.0*DeltaD1w)+.5*(D1wl+D2wl*dNn/1+3.0*DeltaD1w))*dNn;
    wll(NnMidh+1)=wll(NnMidh)+(.5*(D1wl+D2wl*dNn/1+3.0*DeltaD1w)+.5*(D1wl+D2wl*dNn/1+4.0*DeltaD1w))*dNn;

%         
    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+2)-w(NnMidl-2)-5*D1wh*dNn-5.0*D2wh*dNn.^2.0/1.0)/12.50/dNn;
%     DeltaD1w=(w(NnMidh+2)-w(NnMidl-1)-4*D1wh*dNn-16.0*D2wh*dNn.^2.0/1.0)/6.0/dNn;


    whh(NnMidh+1)=w(NnMidh+2)-(.5*(D1wh+D2wh*dNn/1)+.5*(D1wh+D2wh*dNn/1+DeltaD1w))*dNn;
    whh(NnMidh)=whh(NnMidh+1)-(.5*(D1wh+D2wh*dNn/1+1*DeltaD1w)+.5*(D1wh+D2wh*dNn/1+2.0*DeltaD1w))*dNn;
    whh(NnMidl)=whh(NnMidh)-(.5*(D1wh+D2wh*dNn/1+2.0*DeltaD1w)+.5*(D1wh+D2wh*dNn/1+3.0*DeltaD1w))*dNn;
    whh(NnMidl-1)=whh(NnMidl)-(.5*(D1wh+D2wh*dNn/1+3.0*DeltaD1w)+.5*(D1wh+D2wh*dNn/1+4.0*DeltaD1w))*dNn;



    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.
        
        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);
        
        
        
        
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')
 paths=500000;
 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;
 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;
  %   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);
                 %Above is original coordinates monte carlo equation.
                     end
                 end
             end
         end
     end
 end
 YY(YY<0)=0;
 
 sum(YY(:))/paths %origianl coordinates monte carlo average.
 sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates 
 theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average
 
 
 
 
 BinSize=.0075/2;%Please change this bin size accordingly. When data range is too small decrease it.
 %When density range is too large increase it.Roughly, When starting point of SDE is
 %x0=1, set at .0075*2; when starting point is x0=.1, set at .0075/2.0;
 MaxCutOff=30;
 [YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
 plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 title('x0=.25,theta=.1,kappa=1,sigma=.9,gamma=.95,T=.25')
 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
 
User avatar
Amin
Topic Author
Posts: 2122
Joined: July 14th, 2002, 3:00 am

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

Today, 6:58 pm

Here are some more graphs for interest to people in interest rate derivatives. These graphs are for diffusion with x0=theta=1. with volatility exponent gamma .55, .75 and .95 with various different maturities. There is one graph for gamma=.45<.5
Here are the graphs.

Image

Image
Image
Image
Image
Image
Image
Image
Image
Image
Image
Image
Image
Image
Image

and finally at gamma=.45, we have a graph for T=2

Image
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