Breakthrough in the theory of stochastic differential equations and their simulation

Amin · May 24th, 2022, 3:09 pm

Friends, every time I heuristically add derivatives from volatility of volatility term into 2D FPE solution, the series solution becomes a closer match with monte carlo density. I am re-doing some of the analytics so that effect of volatility of volatility term is very properly accounted for as including all such terms may possibly bring the results in line with monte carlo density. I could not work properly for past several days but was able to do some work today. I will let friends know how it goes after formally including all effects from volatility of volatility term in next few days.

Amin · May 27th, 2022, 12:36 pm

Friends, I continue to not succeed in proper solution of 2D FPE that would lead to a density match with monte carlo density. I also feel slightly burnt out and therefore I have tried to take a break for a week in which I would work on my market trading models. Meanwhile, I want to keep thinking about other possibilities to solve the 2D pair of SV SDEs.
Good thing is that, after my exploratory work, many friends have some basic framework of 2D series and they can try their own analytics to find a proper solution of 2D SDE. I wonder why people do not use 2D or 3D power series for solution of PDEs etc. I think they are a very powerful mathematical tool if the higher derivatives of the solution are decreasing fast enough.
Most of the 2D series functions are self-explanatory and mathematicians would not have any problems in understanding and generalizing them. However the function that converts 2D series solution to 1D asset solution probably needs some explanation and I would be doing that here on this forum in next few days.
I have many other ideas for interesting research in stochastics and I would continue to post my work here. Sorry for this break but I want to come back with fresh and better ideas and start posting my research again.

Amin · May 27th, 2022, 8:56 pm

Friends, I will be posting a new version of analytic solution to pair of SV SDEs. It is simpler since it depends only on addition of squared variances across hermite polynomials and yields somewhat better results than solution of Fokker-planck equation I posted earlier. I hope it would be interesting for friends to play around and explore it. It is late now but I would be posting it tomorrow.

Amin · May 29th, 2022, 4:41 am

Friends, I think I know the reason why our programs are not working well for 2D stochastic volatility SDE models. I think the integral of variance

[$] \int \, V \, dt [$]

in all of the solutions have to be taken in sense of a dt-integral and not merely V multiplied with dt. The relevant integral on every time interval would be addition to dt-integral on that time stage.
When I was playing with dt-integrals, I noticed that their volatility was reasonably lesser than simple integrals with multiplication with dt. And I believe that in order to get our programs to work, we need a smaller vol in a similar ratio that is between a properly calculated dt-integral and an integral obtained by simple multiplication of dt.
Good thing is that we did some good work on dt-integrals so that work will be very helpful.
I will be re-doing my programs (both derived from fokker-planck and with simple addition of hermite coefficients) to change them so that a proper dt-integral goes into the variance and then post those program in another 2-3 days.

Amin · May 29th, 2022, 11:02 am

What I mean in previous post is that when the equation below is integrated in time.
[$]=\int_0^{w_b(n)} \, \Big[\, -\mu_y \, P(Z_1) \, \frac{\partial Z_2}{\partial y} \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big] \, \Big] dZ_1 [$]
[$]+ \int_0^{w_b(n)} \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0 [$]

the integral
[$]\int_t^{t+dt} \, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V} dt [$]

has to be solved in the sense of proper stochastic dt-integral and not merely a multiplication with delta t.

Amin · May 30th, 2022, 3:33 pm

Friends, I am quite convinced that difference between a simple integral found by multiplying variance with delta t and a properly calculated stochastic dt-integral is the reason why we are not getting a good match with monte carlo in our analytic solution of stochastic vol SDEs. A proper stochastic dt-integral has to be calculated so that cumulative variance(over all past time steps) matches perfectly with variance of dt-integral on each time step. The integrlas when a stochastic term is simply multiplied with delta t has larger variance than a properly calculated stochastic integral but difference is not extremely huge and I would guess that it is of the same order as discrepancy we have(with monte carlo) in our case.
I have to be busy tomorrow to get my national ID card and I might have to stand in line for several hours. Later tomorrow night my sister is arriving from US after 3-4 years and I would have to be a bit busy but good thing is that I think I have a perfect handle on the analytics now and I will get our 2D program to work in next 2-3 days.

Just taking a very small two step toy problem. Let us suppose
[$] X(t_1) \,= a_{11} \, Z \,[$]
[$] X(t_2) \,= \sqrt{ {a_{11}}^2 + {a_{21}}^2 } \, Z \,[$]

A crude summation of both terms is given as
[$] X(t_1) \, \Delta t + X(t_2) \, \Delta t [$]
[$] = \big [a_{11} \, Z \, +\sqrt{ {a_{11}}^2 + {a_{21}}^2 } \, Z \big] \, \Delta t [$]
which has variance given as
[$] =\big [ {a_{11}}^2 \, +\, { {a_{11}}^2 + {a_{21}}^2 } \,+ 2 \, a_{11} \, \sqrt{ {a_{11}}^2 + {a_{21}}^2 } \big ] \, {\Delta t}^2 [$] Eq(1)

while true variance of dt-integral is given as
[$] =\big [ 4 {a_{11}}^2 \, + {a_{21}}^2 \,\big ] \, {\Delta t}^2[$] Eq(2)

We notice that variance of Eq(1) is slightly larger than that of Eq(2)
We did this calculation with just two terms and adding more terms has a similar effect of increasing the variance of delta t multiplied integral as compared to true variance of stochastic dt-integral. When calculating a proper stochastic dt-integral, the cumulative variance(over all time steps on time grid) must match perfectly with cumulative variance of dt-integral on the same time steps.

Amin · June 2nd, 2022, 5:44 am

Friends, sorry for all this delay but I have been able to verify that analytical 2D SV density remains in line with monte carlo when we treat the integration as a proper stochastic dt-integral on every step. There are still some problems like a huge tail for Asset density but that can be corrected (or hugely improved) by increasing the number of terms (of Z-series appropriately) in our analysis. I will be doing some more work and then post the results and the program later today or tomorrow. But I am very happy that we got close to completing of our project after I was almost losing hope.
I will be posting the new program with explanation of analytics very shortly.

Amin · June 3rd, 2022, 4:54 am

Before posting the graphs from my work, I would like to take this opportunity to thank those friends who with their help made my research possible. Several years ago, I had to spend a huge amount of time everyday by going into remote parts of Lahore city and try to find good food and water. Mind control agencies would drug food, water and beverages in the entire markets and it was extremely difficult to get good food and water. But as I tried to do good research, several American friends tried to help me by asking mind control agencies to become better and stop employing such tactics. Though food and water is drugged at a lot of places even today, it is easier these days than ever before to go out and be careful and easily be able to get good water and food (though I would still get hit occasionally sometimes). All other forms of torture also decreased with time especially when I would mention the torture here and friends would try to help me by asking secret services to become more civilized. I can never imagine that I could do any of the research I posted on this forum if circumstances were as bad as they were several years ago. I simply cannot have enough gratitude for American friends and others who have indirectly helped me in the past and thus made my research possible.

About the 2D SV SDEs, though I have been able to get a good match to monte carlo out to several years, we still have to do quite a bit of work to refine the model. As time increases, there are much larger tails especially when both volatility and asset are close to lognormal. To the extreme right volatility tail gets larger and asset tail has a compound effect on it and it increases by a huge amount but I believe both of these effects can be remarkably decreased by appropriately increasing the order of volatility Z-series and asset Z-series where required. Another reason behind the larger tail is that I used inversion of series to get the asset density and inversion of the series is very approximate in the tails and also contributes to larger tail. Increasing the order of the series would also help with inversion of the series and there could be other more exact ways to get the density without using inversion of the series.
I have used stochastic time integrals in the program and with mean-reversion, they are not so accurate after 4-5 years in the current program. I am sure this is not a major problem and we will easily solve this(I just could not find time yesterday to improve the previous work). There are also a lot of minor issues, there are some higher order effects (possibly second order) that I have yet not included in pricing and this also needs improvement. But the very good thing is that we have solved most of the difficult issues and work remains to be done is relatively very easy. But still the current work is only intermediate and quite some of the work(though relatively easy) still remains to be done.

Here are some graphs from the new program for 2D SV pricing. All the parameters of the two SDEs are shown in the title of the graphs. I will be posting the very intermediate matlab program used to make these graphs in a few hours.

Amin · June 3rd, 2022, 1:38 pm

Friends, the program below would work with mean-reverting volatility only when V0=theta.
The purpose behind writing this program was only to test whether calculating a proper dt-integral would fix our problem of discrepancy of analytical density with respect to monte carlo density. I hope to fix all the remaining shortcomings in next few days and will make it a very general program. It is just an intermediate version. I have posted all other necessary functions earlier.

.
.

function [] = FPESeriesCoeffsFromGaussiansWilmottNew02A_Lite_2DFPEA()

%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334 
%or skype ahsan.amin2999

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

dt=.125/2/2/2;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.
Tt=64*1;%16*2*4;%*4*4*1;%*4;%128*2*2*2;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year; 
T=Tt*dt;
% 
% %Below I have done calculations for smaller step size at start.
% 
% %ds(1:64)=dt/16;
% %ds(65:128)=dt/4;
% if(Tt<=4)
% Ts=(Tt*16);%+(64-16);
% ds(1:Ts)=dt/16;
% end
% if(Tt>4)&&(Tt<=20)
% Ts=(64)+((Tt-4)*4);
% ds(1:64)=dt/16;
% ds(65:Ts)=dt/4;
% end
% if(Tt>20)
% Ts=(64)+((20-4)*4)+(Tt-20);    
% ds(1:64)=dt/16;
% ds(65:128)=dt/4;
% ds(129:Ts)=dt;
% end

Ts=Tt;
ds(1:Tt)=dt;

T
sum(ds(1:Ts))
Ts

%Ts=4;
str=input('Check if time is right');

OrderA=4;  %
OrderM=4;  %

if(T>=1)
dtM=1/32;%.0625/1;
TtM=T/dtM;
else
dtM=T/32;%.0625/1;
TtM=T/dtM;
end

%dtM=dt;
%TtM=Tt;

dNn=.2/2;   % Normal density subdivisions width. would change with number of subdivisions
Nn=90;  % No of normal density subdivisions

NnMid=4.0;

v00=1.0;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.950;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=1.0;%.950;   %mean reversion parameter.
theta=1.0;%mean reversion target
sigma0=.9500;%Volatility value
NoOfCumulants=6;
SeriesOrder=NoOfCumulants-1;

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

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


x0=1.00;
gammaX=.75;
gammaV=.5;
sigmaX=.350;
thetaX=1.0;
kappaX=0.0;
mu1X=thetaX*kappaX;
mu2X=-1*kappaX;
beta1X=0.0;
beta2X=1.0;
q0=x0^(1-gammaX)/(1-gammaX);
alpha1X=1-gammaX;



Z(1:Nn)=(((1:Nn)-6.5/1)*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);
ZProb(2:Nn-1)=normcdf(.5*Z(2:Nn-1)+.5*Z(3:Nn),0,1)-normcdf(.5*Z(2:Nn-1)+.5*Z(1:Nn-2),0,1);


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

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


sigma11X(1:OrderA+1)=0;
mu11X(1:OrderA+1)=0;
mu22X(1:OrderA+1)=0;
sigma22X(1:OrderA+1)=0;

sigma11X(1)=1;
mu11X(1)=1;
mu22X(1)=1;
sigma22X(1)=1;

for k=1:(OrderA+1)
    if sigmaX~=0
        sigma11X(k)=sigmaX^(k-1);
    end
    if mu1X ~= 0
        mu11X(k)=mu1X^(k-1);
    end
    if mu2X ~= 0
        mu22X(k)=mu2X^(k-1);
    end
    if sigmaX~=0
        sigma22X(k)=sigmaX^(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.

Fp1X(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
%YCoeff0 and YCoeff are coefficents for original coordinates monte carlo.
%YqCoeff0 and YqCoeff are bessel/lamperti version monte carlo.

YCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
YqCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
YqCoeff0X(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;

%Pre-compute the time and power exponent values in small multi-dimensional arrays
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma); %expand y^alpha where alpha=1;
YqCoeff = ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);%expand y^alpha1 where alpha1=(1-gamma)
YqCoeff=YqCoeff/(1-gamma); %Transformed coordinates coefficients have to be 
%further divided by (1-gamma)
YqCoeffX = ItoTaylorCoeffsNew(alpha1X,beta1X,beta2X,gammaX);%expand y^alpha1 where alpha1=(1-gamma)
YqCoeffX=YqCoeffX/(1-gammaX); 
for k = 0 : (OrderA)
    for m = 0:k
        l4 = k - m + 1;
        for n = 0 : m
            l3 = m - n + 1;
            for j = 0:n
                l2 = n - j + 1;
                l1 = j + 1;
                %Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
                Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                
                Fp1X(l1,l2,l3,l4) = (alpha1X + (l1-1) * beta1X + (l2-1) * beta2X + (l3-1) * 2* gammaX + (l4-1) * gammaX ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                
                
                YCoeff0(l1,l2,l3,l4) =YCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
                YqCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
                YqCoeff0X(l1,l2,l3,l4) =YqCoeffX(l1,l2,l3,l4).*mu11X(l1).*mu22X(l2).*sigma22X(l3).*sigma11X(l4);
            end
        end
    end
end


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

ExpnOrder=5;
SeriesOrder=5;
c(1:SeriesOrder)=0;
c0=w0;
wnStart=1;
yydt(wnStart:Nn)=0.0;
tic

for tt=1:Ts
    tt

    if(tt<=1)
        [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
    cprev(1:5)=0;
    c0prev=w0;
        c0=wMu0dt+w0;
        c(1)=sigma0*sqrt(ds(tt));
        c(2:7)=0.0;
    
        w0=c0;
    elseif(tt<=3) 

    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),8);
     [b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
     
    c0prev=c0;
    cprev=c;
    c0=c0+b0;
     c(1)=b(1)+sqrt(c(1)^2+sigma0^2*ds(tt));
     c(2:SeriesOrder)=b(2:SeriesOrder);
     
     w0=c0;
    else
        
     

        
        
     [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
     [b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
     
     
     
    [wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
    
    [wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
    
    [g10,g1] = CalculateDriftbCoeffs08A(wVol0dt,dwVol0dtdw,c,SeriesOrder);
    g10=g10+sigma0*sqrt(ds(tt));
    [g20,g2] = CalculateDriftbCoeffs08A(wVol2dt,dwVol2dtdw,c,SeriesOrder);
    
cprev=c;
c0prev=c0;
  a0=c0+b0;
  a(1:SeriesOrder)=c(1:SeriesOrder)+b(1:SeriesOrder);
%     
      [a0,a] = AdvanceSeriesCumulantSolution02NewO6Params5A_Lite_0(a0,a,g10,g1,g20,g2,sigma0,gamma,ds(tt),SeriesOrder,NoOfCumulants);

 c0=a0;
 c(1:SeriesOrder)=a(1:SeriesOrder);
    
    w0=c0;
    end
    
     if(tt==1)
        cmid0=c0;
        cmid=c;
        
     else
       cmid0=(c0+c0prev)*.5;
        cmid=(c+cprev)*.5;  
        cmid0=c0;
        cmid=c;  
        
     end
    
 
     
      gammaV2=2*gammaV;
      [vw0,dvw0] = CalculatevgammaV2AndDerivatives(c0prev,gamma,gammaV2,SeriesOrder);
%      %[V2gamma0,V2gamma1] = CalculateDriftbCoeffs08A(vw0,dvw0,cmid,SeriesOrder);
      [V2gamma0,V2gamma1] = CalculateDriftbCoeffs08A(vw0,dvw0,cprev,SeriesOrder);
      
      V2gamma(1)=V2gamma0;
      V2gamma(2:6)=V2gamma1(1:5);
%      
     
     
    if(tt==1)

        cq(1:6,1:6)=0;    
        cq(1,1)=x0^(1-gammaX)/(1-gammaX);
        
        
    end
    %    cprev=c;
    
    q0=cq(1,1);

    
    if(tt==1)
        
        gammaV2=2*gammaV;
     [vw0,dvw0] = CalculatevgammaV2AndDerivatives(c0prev,gamma,gammaV2,SeriesOrder);
     %[V2gamma0,V2gamma1] = CalculateDriftbCoeffs08A(vw0,dvw0,cmid,SeriesOrder);
     [ch0,ch] = CalculateDriftbCoeffs08A(vw0,dvw0,cprev,SeriesOrder);
     [ch0,ch] = ConvertZCoeffsToHCoeffs(ch0,ch,SeriesOrder);    
    
     cH0(tt)=ch0;
     cH(tt,1:5)=ch(1:5);
     
%     tt
     
     
    else
    [cH0,cH,cVdt] = CalculateVPathStepIntegral(c0prev,cprev,c0,c,kappa,theta,gamma,dt,tt,cH0,cH,SeriesOrder);

    
    end
    
    
    

    DoubleExpnOrder=ExpnOrder*2;
    
    [qMu0dt,dqMu0dtdq,qMu1dt,dqMu1dtdq] = BesselDriftAndDerivatives08A0(q0,mu1X,mu2X,beta1X,beta2X,gammaX,DoubleExpnOrder);
    [q2Mu0dt,dq2Mu0dtdq,q2Mu1dt,dq2Mu1dtdq] = BesselDriftAndDerivatives08Aq(q0,sigmaX,gammaX,DoubleExpnOrder);
 
    [Muq0] = CalculateDriftbCoeffs08A2Dim02(qMu0dt,dqMu0dtdq,cq,SeriesOrder);
    [Muq2] = CalculateDriftbCoeffs08A2Dim02(q2Mu0dt,dq2Mu0dtdq,cq,SeriesOrder);
    [DMuq0] = CalculateDriftbCoeffs08A2Dim02(qMu1dt,dqMu1dtdq,cq,SeriesOrder);
    [DMuq2] = CalculateDriftbCoeffs08A2Dim02(q2Mu1dt,dq2Mu1dtdq,cq,SeriesOrder);
    
 
    
    
    [Muq2] =SeriesProduct2D1D(Muq2,V2gamma);
    [DMuq2] =SeriesProduct2D1D(DMuq2,V2gamma);
 
    
    
    Muq=Muq0+Muq2;%+Muq01*dt/2;%-tt*Muq02*dt/2+tt*Muq01*dt/2;
    DMuq=DMuq0+DMuq2;
    
    if(tt==1)
        cq=cq+Muq*dt+SeriesProduct2D(Muq,DMuq)*dt^2/2;
        cq(2,1)=cq(2,1)+sigmaX*v00^gammaV*sqrt(dt);
    elseif(tt<=3)
%        cq
        [D1cq] = Series2DZ1Derivative(cq);
%        D1cq
        [D1cqInv] = Series2DReciprocal02(D1cq);
%        D1cqInv
        [D1cqInvZ] = MultiplySeries2DWithZ1(D1cqInv);
%        D1cqInvZ
        fq=Muq+.5*sigmaX^2*SeriesProduct2D1D(D1cqInvZ,V2gamma);
%        fq*dt
        Dfq=DMuq;
        cq=cq+fq*dt+SeriesProduct2D(fq,Dfq)*dt^2/2;
        
        
        SeriesProduct2D1D(D1cqInvZ,V2gamma)
%        str=input('Look at numbers-1111');
%        D1cqInvZ
%        str=input('Look at numbers-2222');
    else
        [D1cq] = Series2DZ1Derivative(cq);
    
        [D1cqInv] = Series2DReciprocal02(D1cq);
        [D1cqInvZ] = MultiplySeries2DWithZ1(D1cqInv);
        [D2cq] = Series2DZ1Derivative(D1cq);
        [D1cqInvSqr]=SeriesProduct2D(D1cqInv,D1cqInv);
        [D1cqInvCub]=SeriesProduct2D(D1cqInvSqr,D1cqInv);
        [D2cqD1cqInvSqr]=SeriesProduct2D(D2cq,D1cqInvSqr);
        [D2cqD1cqInvCub]=SeriesProduct2D(D2cq,D1cqInvCub);
        D2cqT2=D1cqInvSqr-MultiplySeries2DWithZ1(D2cqD1cqInvCub);
        Dfq=DMuq*dt+.5*sigmaX^2*SeriesProduct2D1D(D2cqT2,cVdt);
        D2cqT1=D1cqInvZ+D2cqD1cqInvSqr;
        

        fq=Muq*dt+.5*sigmaX^2*SeriesProduct2D1D(D2cqT1,cVdt);    
 
       cq=cq+(fq+SeriesProduct2D(fq,Dfq)*dt^0/2);%.*cprev(1)/c(1);
        
    end
 
 


end

Z1=Z;

if(tt*dt>1)
    cq(:,5)=cq(:,5)/4;
    cq(5,:)=cq(5,:)/4;
end

[b] = EvaluateSeriesForZ2(cq,Z1,Nn);
% Mm=301;
% MmMid=151;
% Q(MmMid+1:Mm)=cq(1,1)+cq(2,1).*.1.*(1:MmMid-1);
% Q(MmMid)=cq(1,1);
% %Q(1:70)=cq(1,1)-cq(2,1).*.25.*(70:-1:1);
% for mm=1:MmMid-1
% 
%     Q(MmMid-mm)=cq(1,1)-cq(2,1).*.1.*mm;
% end
% 
% 
% Qa=Q-cq(2,1).*.1/2;
% Qb=Q+cq(2,1).*.1/2;



Mm=601;
MmMid=301;
dMm=.05;
Q(MmMid+1:Mm)=cq(1,1)+cq(2,1).*dMm.*(1:MmMid-1);
Q(MmMid)=cq(1,1);
%Q(1:70)=cq(1,1)-cq(2,1).*.25.*(70:-1:1);
for mm=1:MmMid-1

    Q(MmMid-mm)=cq(1,1)-cq(2,1).*dMm.*mm;
end


Qa=Q-cq(2,1).*dMm/2;
Qb=Q+cq(2,1).*dMm/2;






Z1Prob(1:Mm)=0;
for nn=1:Nn
    bq(1:6)=b(1:6,nn);

    
    [Z0Prob] = CalculateProbabilitySeriesInv(Qa,Qb,bq);
    
%    Z0Prob
    
%str=input('Look at numbers-tttt')
    
    
    Z1Prob(1:Mm)=Z1Prob(1:Mm)+Z0Prob(1:Mm).*ZProb(nn);
    
end

pQ=Z1Prob/(cq(2,1).*dMm);
    
xx(1:Mm)=((1-gammaX).*Q(1:Mm)).^(1/(1-gammaX));

pxx(1:Mm)=pQ(1:Mm).*xx(1:Mm).^(-gammaX);
    
    
Mean=sum(xx(1:Mm).*Z1Prob(1:Mm));
Mean
  w(1:Nn)=c0;
  for nn=1:SeriesOrder
      w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^nn;
  end
  
  
 % q(1:Nn)=cq0;
 % for nn=1:SeriesOrder
 %     q(1:Nn)=q(1:Nn)+cq(nn)*Z(1:Nn).^nn;
 % end
  
  
  
%  %Flag=0;
%  %for nn=ceil(Nn/2)-1:-1:1
%  %    if((w(nn)<0)&&(Flag==0))
%  %        wnStart=nn+1;
%  %        Flag=1;
%  %    end
%  %end
%  
%  
%  %c0
%  %c
%  
%  
%  %str=input('Look at numbers')
  yy0(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
  
 % xx0(wnStart:Nn)=((1-gammaX).*q(wnStart:Nn)).^(1/(1-gammaX));
%  
%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%   yw0=((1-gamma).*c0).^(1/(1-gamma));
%   dyw0(1)=((1-gamma).*c0).^(1/(1-gamma)-1);
%   dyw0(2)=(1/(1-gamma)-1).*((1-gamma).*c0).^(1/(1-gamma)-2)*(1-gamma);
%   dyw0(3)=(1/(1-gamma)-1).*(1/(1-gamma)-2).*((1-gamma).*c0).^(1/(1-gamma)-3)*(1-gamma)^2;
%   dyw0(4)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*((1-gamma).*c0).^(1/(1-gamma)-4)*(1-gamma)^3;
%   dyw0(5)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*((1-gamma).*c0).^(1/(1-gamma)-5)*(1-gamma)^4;
%   dyw0(6)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*((1-gamma).*c0).^(1/(1-gamma)-6)*(1-gamma)^5;
%   dyw0(7)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*(1/(1-gamma)-6)*((1-gamma).*c0).^(1/(1-gamma)-7)*(1-gamma)^6;
%   dyw0(8)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*(1/(1-gamma)-6)*(1/(1-gamma)-7)*((1-gamma).*c0).^(1/(1-gamma)-8)*(1-gamma)^7;
% %  
% % c(6)=0;
% % c(7)=0;
% % c(8)=0;
% % 
%      [y10,y1] = CalculateDriftbCoeffs08A(yw0,dyw0,c,SeriesOrder);
%      [y20,y2] = CalculateDriftbCoeffs08A(yw0,dyw0,c,5);
% %  
%  
%  yy1(1:Nn)=y10;
%  %y(1:SeriesOrder)=c(1:SeriesOrder);
%  for nn=1:SeriesOrder
%      yy1(1:Nn)=yy1(1:Nn)+y1(nn)*Z(1:Nn).^nn;
%  end
%  
%  yy2(1:Nn)=y20;
%  %y(1:SeriesOrder)=c(1:SeriesOrder);
%  for nn=1:5
%      yy2(1:Nn)=yy2(1:Nn)+y2(nn)*Z(1:Nn).^nn;
%  end
%  plot(Z(1:Nn),yy0(1:Nn),'r',Z(1:Nn),yy1(1:Nn),'b',Z(1:Nn),yy2(1:Nn),'k');
% title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f', x0,theta,kappa,gamma,sigma0,T,dt));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
%  
% legend({'SDE Variable Converted on Grid','SDE Variable From Eight Cumulants','SDE Variable From Six Cumulants'},'Location','northeast')
%  
% %str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
% 
%  
%  str=input('red line is SDE from Ito-Hermite method.');
%  
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
 yy=yy0;
% xx=xx0;
Dfyy(wnStart:Nn)=0;
%Dfyy1(wnStart:Nn)=0;
%Dfyy2(wnStart:Nn)=0;
%Dfxx(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    
    Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%    Dfyy1(nn) = (yy1(nn + 1) - yy1(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%    Dfyy2(nn) = (yy2(nn + 1) - yy2(nn - 1))/(Z(nn + 1) - Z(nn - 1));
 %   Dfxx(nn) = (xx(nn + 1) - xx(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    
    %Change of variable derivative for densities
end
%Dfyy(Nn)=Dfyy(Nn-1);
%Dfyy(1)=Dfyy(2);
%Dfyydt(Nn)=Dfyydt(Nn-1);
%Dfyydt(1)=Dfyydt(2);


pyy(1:Nn)=0;
%pyy1(1:Nn)=0;
%pyy2(1:Nn)=0;
%pxx(1:Nn)=0;
for nn = wnStart+1:Nn-1
    
    pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
%    pyy1(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy1(nn));
%    pyy2(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy2(nn));
 %   pxx(nn) = (normpdf(Z(nn),0, 1))/abs(Dfxx(nn));
end


toc

ItoHermiteMeanVar=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
ItoHermiteMeanAsset=Mean
disp('true Mean only applicable to standard SV mean reverting type models otherwise disregard');
%TrueMean=theta+(v0-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average
%yy0


theta1=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%rng(29079137, 'twister')
rng(15898837, 'twister')
paths=100000;
V(1:paths)=v00;  %Original process monte carlo.
X=0.0;
X(1:paths)=x0;
alpha1=0;
alpha2=1;
a=mu1X;
b=mu2X;
rho=0;
sigma1=sigmaX;
gammaV=.5;
Random1(1:paths)=0;
Random2(1:paths)=0;

for ttM=1:TtM
Random1=randn(size(Random1));
Random2=randn(size(Random2));



X(1:paths)=X(1:paths)+ ...
    (a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dtM + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) * sqrt(dtM) + ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*sqrt(dtM) + ...
    (a*alpha1* X(1:paths).^(alpha1-1)+b*alpha2* X(1:paths).^(alpha2-1)).* ...
    (((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dtM^2/2)+ ...
    (sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) *(1-1/sqrt(3)).*dtM^1.5+ ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*(1-1/sqrt(3)).*dtM^1.5))+ ...
    .5*(a*alpha1*(alpha1-1)* X(1:paths).^(alpha1-2)+b*alpha2*(alpha2-1).* X(1:paths).^(alpha2-2)).* ...
    ( sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX)) *dtM^2/2 + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
    ((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random1(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*(Random1(1:paths).^2-1) * dtM/2 + ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths)*dtM/2)+ ...
    .5*sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
    (sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random1(1:paths).*1/sqrt(3).*dtM^1.5)+ ...
    sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
    ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random1(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
    sigma0*V(1:paths).^gamma.*Random1(1:paths).*Random2(1:paths)*dtM/2)+ ...
    .5*sqrt(1-rho^2)* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
    (sigma0^2*V(1:paths).^(2*gamma).*Random1(1:paths)*1/sqrt(3)*dtM^1.5)+ ...
    sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
    rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random1(1:paths)*1/sqrt(3)*dtM^1.5+ ...
    rho* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
    ((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random2(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths) * dtM/2 + ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*(Random2(1:paths).^2-1)*dtM/2)+ ...
    .5*rho* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
    (sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random2(1:paths).*1/sqrt(3).*dtM^1.5)+ ...
    rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
    ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
    sigma0*V(1:paths).^gamma.*(Random2(1:paths).^2-1)*dtM/2)+ ...
    .5*rho* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
    sigma0^2.*V(1:paths).^(2*gamma).*Random2(1:paths) * 1/sqrt(3).* dtM^1.5+ ...
    rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
    rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random2(1:paths)*1/sqrt(3)*dtM^1.5;
    
VBefore=V;

    V(1:paths)=V(1:paths)+ ...
        (mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dtM + ...
        sigma0*V(1:paths).^gamma .*Random2(1:paths)*sqrt(dtM) + ...
        (mu1.*beta1*V(1:paths).^(beta1-1) + mu2.*beta2.*V(1:paths).^(beta2-1)).* ...
        ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dtM^2/2 + ...
        sigma0*V(1:paths).^gamma .*Random2(1:paths)*(1-1/sqrt(3))*dtM^1.5) + ...
        .5*(mu1.*beta1.*(beta1-1).*V(1:paths).^(beta1-2) + mu2.*beta2.*(beta2-1).*V(1:paths).^(beta2-2)).* ...    
        sigma0^2.*V(1:paths).^(2*gamma).*dtM^2/2 + ...
        sigma0*gamma*V(1:paths).^(gamma-1) .* ...
        ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths).*1/sqrt(3)*dtM^1.5 + ...
        sigma0.*V(1:paths).^gamma .*(Random2(1:paths).^2-1)*dtM/2) + ...
        .5*sigma0*gamma*(gamma-1).*V(1:paths).^(gamma-2) .* ...
        sigma0^2.*V(1:paths).^(2*gamma) .*Random2(1:paths).*1/sqrt(3)*dtM^1.5;

%Vm=.5*VBefore+.5*V;
Vm=V;
    
%    X(1:paths)=X(1:paths)+ ...
%    sigma1* Vm(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) * sqrt(dtM) + ...
%    sigma1* Vm(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
%    (sigma1* Vm(1:paths).^gammaV.* X(1:paths).^gammaX .*(Random1(1:paths).^2-1) * dtM/2)+ ...
%    .5*sigma1* Vm(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
%    (sigma1^2* Vm(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random1(1:paths).*1/sqrt(3).*dtM^1.5);%+ ...

    
    
    
end

%SVolMeanAnalytic=thetaV+(V0-thetaV)*exp(-kappaV*dt*Tt)
SVolMeanAnalytic=theta+(v00-theta)*exp(-kappa*dt*Tt)
SVolMeanMC=sum(V(1:paths))/paths
AssetMeanAnalytic=x0
AssetMeanMC=sum(X(1:paths))/paths

% theta
% v00
% kappa
% dt
% Tt
% T


%MeanX

disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(v00-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average


MaxCutOff=30;
NoOfBins=2*300;%round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too 
[XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti_NEW(X,paths,NoOfBins,MaxCutOff );
plot(xx(wnStart+1:Mm-1),pxx(wnStart+1:Mm-1),'r',IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g');
 %plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');
 
title(sprintf('x0 = %.2f,thetaX=%.2f,kappaX=%.2f,gammaX=%.2f,sigmaX=%.2f,v0 =%.2f,kappa=%.2f,theta=%.2f,gamma=%.2f,sigma0=%.2f,T=%.2f,dt=%.3f,M=%.4f',x0,thetaX,kappaX,gammaX,sigmaX,v00,kappa,theta,gamma,sigma0,T,dt,ItoHermiteMeanAsset));
%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
 
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');


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

NoOfBins=round(2*500*gamma^2*4*sigma0/sqrt(SVolMeanMC)/(1+kappa))/10;
[VDensity,IndexOutV,IndexMaxV] = MakeDensityFromSimulation_Infiniti_NEW(V,paths,NoOfBins,MaxCutOff );

%[YdtDensity,IndexOutYdt,IndexMaxYdt] = MakeDensityFromSimulation_Infiniti_NEW(YYdt,paths,NoOfBins,MaxCutOff );
plot(yy(wnStart+1:Nn-1),pyy(wnStart+1:Nn-1),'r',IndexOutV(1:IndexMaxV),VDensity(1:IndexMaxV),'g');
title(sprintf('v0 =%.2f,kappa=%.2f,theta=%.2f,gamma=%.2f,sigma0=%.2f,T=%.2f,dt=%.3f,M=%.4f',v00,kappa,theta,gamma,sigma0,T,dt,ItoHermiteMeanVar));



end

.
.
.

function [cH0,cH,cVdt] = CalculateVPathStepIntegral(c0prev,cprev,c0,c,kappa,theta,gamma,dt,Ts,cH0,cH,SeriesOrder)


    gammaV2New=1;
    [vw0,dvw0] = CalculatevgammaV2AndDerivatives(c0,gamma,gammaV2New,SeriesOrder);
    %[V2gamma0,V2gamma1] = CalculateDriftbCoeffs08A(vw0,dvw0,cmid,SeriesOrder);
    [cv0,cv] = CalculateDriftbCoeffs08A(vw0,dvw0,c,SeriesOrder);

    gammaV2New=1;
    [vw0,dvw0] = CalculatevgammaV2AndDerivatives(c0prev,gamma,gammaV2New,SeriesOrder);
    %[V2gamma0,V2gamma1] = CalculateDriftbCoeffs08A(vw0,dvw0,cmid,SeriesOrder);
    [cv0prev,cvprev] = CalculateDriftbCoeffs08A(vw0,dvw0,cprev,SeriesOrder);
    
   % b0=cv0prev*(exp(-kappa*dt)-1)-theta*(exp(-kappa*dt)-1);
   % b=cvprev*(exp(-kappa*dt)-1);

    [fH0,fH] = ConvertZCoeffsToHCoeffs(cv0,cv,SeriesOrder);    
    
    [fH0prev,fHprev] = ConvertZCoeffsToHCoeffs(cv0prev,cvprev,SeriesOrder);    
 
 
   % [bH0,bH] = ConvertZCoeffsToHCoeffs(b0,b,SeriesOrder);   

    dH0=fH0-fH0prev;
 

    %dH(1:SeriesOrder)=sign(fH(1:SeriesOrder)-bH(1:SeriesOrdery)-exp(-1*kappa*dt).*(fHprev(1:SeriesOrder))).*sqrt(abs(sign(fH(1:SeriesOrdery)-bH(1:SeriesOrder)).*(fH(1:SeriesOrdery)-bH(1:SeriesOrdery)).^2-sign(exp(-kappa*dt).*(fHprev(1:SeriesOrdery))).*(exp(-2*1*kappa*dt).*(fHprev(1:SeriesOrdery)).^2))); 
    dH(1:SeriesOrder)=sign(fH(1:SeriesOrder)-exp(-1*kappa*dt).*(fHprev(1:SeriesOrder))).*sqrt(abs(sign(fH(1:SeriesOrder)).*(fH(1:SeriesOrder)).^2-sign(exp(-kappa*dt).*(fHprev(1:SeriesOrder))).*(exp(-2*1*kappa*dt).*(fHprev(1:SeriesOrder)).^2)));    
 
    cH0(Ts)=dH0;
    cH(Ts,1:SeriesOrder)=dH(1:SeriesOrder);
 
 
     Coeff1(1:Ts)=0;
    Coeff2(1:Ts)=0;

     
      for tt=1:Ts-1
      for ss=tt:Ts-1
          
          Coeff1(tt)=Coeff1(tt)+(exp(0*kappa*dt*(tt))).*(exp(-2*kappa*dt*(ss-tt))).*dt^2;    
          for rr=tt:Ts-1
             if(rr==tt)
             
             else
             Coeff1(tt)=Coeff1(tt)+1*(exp(-1*kappa*dt*(ss)))*(exp(-0*kappa*dt*(tt))).*dt^2;
             end
          end
      end
      end
 
     for tt=1:Ts
     for ss=tt:Ts
         
         Coeff2(tt)=Coeff2(tt)+(exp(0*kappa*dt*(tt))).*(exp(-2*kappa*dt*(ss-tt))).*dt^2;    
         for rr=tt:Ts
            if(rr==tt)
            
            else
            Coeff2(tt)=Coeff2(tt)+1*(exp(-1*kappa*dt*(ss)))*(exp(-0*kappa*dt*(tt))).*dt^2;
            end
         end
     end

     end

 %   Coeff0(1)=Coeff2(1);

     for tt=2:Ts
%         Coeff0(tt)=(Coeff2(tt)-Coeff2(tt-1));
     end
    
    
%     for tt1=1:tt-1
%         for ss=tt1:tt-1
%             for rr=ss:tt-1
%                 if(rr==ss)
%             
%                     Coeff1(tt1)=Coeff1(tt1)+(exp(0*kappa*dt*(tt1))).*(exp(-2*kappa*dt*(rr-ss))).*dt^2;    
%                 else
%            
%                     Coeff1(tt1)=Coeff1(tt1)+2*(exp(-2*kappa*dt*(rr-ss)))*(exp(-1*kappa*dt*(tt1))).*dt^2;
%                 end
%             end
%          
%         end
%      
%     end
%     for tt1=1:tt
%         for ss=tt1:tt
%             for rr=ss:tt
%                 if(rr==ss)
%             
%                     Coeff2(tt1)=Coeff2(tt1)+(exp(0*kappa*dt*(tt1))).*(exp(-2*kappa*dt*(rr-ss))).*dt^2;    
%                 else
%            
%                     Coeff2(tt1)=Coeff2(tt1)+2*(exp(-2*kappa*dt*(rr-ss)))*(exp(-1*kappa*dt*(tt1))).*dt^2;
%                 end
%             end
%          
%         end
%      
%     end
    
    eH10=0;
    eH20=0;
    eH1(1:SeriesOrder)=0;
    eH2(1:SeriesOrder)=0;
    for tt=2:Ts
    ts=Ts-1-tt+1;
    eH10=eH10+cH0(tt)*ts*dt;
    
   % if(tt*dt*kappa<4)
    eH1(1:SeriesOrder)=eH1(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*Coeff1(tt).*cH(tt,1:SeriesOrder).^2;   
%    eH1(1:SeriesOrder)=eH1(1:SeriesOrder)+sign(cH(tt1,1:SeriesOrder)).*Coeff0(tt1).*cH(tt1,1:SeriesOrder);   
   % else
   % %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery).*(1+exp(2*kappa*dt).*dt^2);
   % eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery).*(1+sqrt(2)/kappa*dt^2);
    %end
    end 
eH0=0;
eH2(1:SeriesOrder)=0;
for tt=1:Ts
    ts1=Ts-1-tt+1;
    ts2=Ts-tt+1;
    eH0=eH0+cH0(tt)*(ts2-ts1)*dt;
    
   % if(tt*dt*kappa<4)
    eH2(1:SeriesOrder)=eH2(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*Coeff2(tt).*cH(tt,1:SeriesOrder).^2;   
   % else
   % %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery).*(1+exp(2*kappa*dt).*dt^2);
   % eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery).*(1+sqrt(2)/kappa*dt^2);
    %end
end


%eH1
%eH2

%Ts
%str=input('Look at two numbers ---000');



eH(1:SeriesOrder)=0.0;
%eH0=0;

%eH20=0;
%eH10=0;


% for tt=1:Ts
%     ts=(Ts-tt+1)*dt;
%     eH2(1:SeriesOrder)=eH2(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*ts^2.*cH(tt,1:SeriesOrder).^2;   
% %    eH0=eH0+exp(-2.*kappa*dt).*cH0(tt)*dt;%*ts*dt;
%     eH20=eH20+cH0(tt)*ts;
%  
% %     if(tt>=1)
% %     eH2(1:SeriesOrder)=eH2(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*Coeff2(tt).*cH(tt,1:SeriesOrder).^2;   
% %     else
% %     eH2(1:SeriesOrder)=eH2(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*Coeff2(tt).*cH(tt,1:SeriesOrder).^2;       
% %     end
%  
% end
% for tt=1:Ts-1
%  %   ts=Ts-tt+1;
%     
%  %   eH0=eH0+exp(-2.*kappa*dt).*cH0(tt)*dt;%*ts*dt;
%     ts=(Ts-tt)*dt;
%     eH1(1:SeriesOrder)=eH1(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*ts^2.*cH(tt,1:SeriesOrder).^2;   
% %    eH0=eH0+exp(-2.*kappa*dt).*cH0(tt)*dt;%*ts*dt;
%     eH10=eH10+cH0(tt)*ts;
%  
% %     if(tt>=1)
% %     eH1(1:SeriesOrder)=eH1(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*Coeff1(tt).*cH(tt,1:SeriesOrder).^2;   
% %     else
% %     eH1(1:SeriesOrder)=eH1(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*Coeff1(tt).*cH(tt,1:SeriesOrder).^2;       
% %     end
%  
% end
% 

%eH0=eH20-eH10;
%eH0=dH0;
%eH0=fH0*dt;

eH(1:SeriesOrder)=sign(eH2(1:SeriesOrder)).*sqrt(abs(sign(eH2(1:SeriesOrder)).*eH2(1:SeriesOrder)))-sign(eH1(1:SeriesOrder)).*sqrt(abs(sign(eH1(1:SeriesOrder)).*eH1(1:SeriesOrder)));
%eH(1:SeriesOrder)=sign(eH2(1:SeriesOrder)).*sqrt(abs(sign(eH2(1:SeriesOrder)).*eH2(1:SeriesOrder).^2))-sign(eH1(1:SeriesOrder)).*sqrt(abs(sign(eH1(1:SeriesOrder)).*eH1(1:SeriesOrder).^2));
%eH(1:SeriesOrder)=sign(eH2(1:SeriesOrder)-eH1(1:SeriesOrder)).*sqrt(abs(sign(eH2(1:SeriesOrder)).*eH2(1:SeriesOrder).^2-sign(eH1(1:SeriesOrder)).*eH1(1:SeriesOrder).^2));

%eH(1:SeriesOrder)=sign(eH2(1:SeriesOrder)).*sqrt(abs(eH2(1:SeriesOrder)));

%eH(1:SeriesOrder)=sign(eH2(1:SeriesOrder)-eH1(1:SeriesOrder)).*sqrt(abs(eH2(1:SeriesOrder)-eH1(1:SeriesOrder)));

%eH0=eH20-eH10;
%eH0=eH20;
%eH0=cH0(tt);
 
 %We convert hermite represenation given by bH0 and bH(1:5) into 
 %dt-integral Z-series Coeffs, e0 and e(1:5)
 %eH(SeriesOrdery)=eH(SeriesOrdery)/2;
[c0vdt,cvdt] = ConvertHCoeffsToZCoeffs(eH0,eH,SeriesOrder);%  
%[c0vdt,cvdt] = ConvertHCoeffsToZCoeffs(cH0,cH,SeriesOrder);%  

 cVdt(1)=c0vdt;
 cVdt(2:6)=cvdt(1:5);

 
 %cVdt
 
 
% cVdt(1)=cv0*dt;
% cVdt(2:6)=cvprev(1:5)*dt;

% cVdt
 
 
% cH
% Ts
 
% str=input('Look at numbers');

end

.
.
.

function [vw0,dvw0] = CalculatevgammaV2AndDerivatives(c0,gamma,gammaV2,SeriesOrder)


    
    vw0=((1-gamma).*c0).^(gammaV2/(1-gamma));
    dvw0(1)=gammaV2*((1-gamma).*c0).^(gammaV2/(1-gamma)-1);
    dvw0(2)=gammaV2*(gammaV2/(1-gamma)-1).*((1-gamma).*c0).^(gammaV2/(1-gamma)-2)*(1-gamma);
    dvw0(3)=gammaV2*(gammaV2/(1-gamma)-1).*(gammaV2/(1-gamma)-2).*((1-gamma).*c0).^(gammaV2/(1-gamma)-3)*(1-gamma)^2;
    dvw0(4)=gammaV2*(gammaV2/(1-gamma)-1).*(gammaV2/(1-gamma)-2)*(gammaV2/(1-gamma)-3)*((1-gamma).*c0).^(gammaV2/(1-gamma)-4)*(1-gamma)^3;
    dvw0(5)=gammaV2*(gammaV2/(1-gamma)-1).*(gammaV2/(1-gamma)-2)*(gammaV2/(1-gamma)-3)*(gammaV2/(1-gamma)-4)*((1-gamma).*c0).^(gammaV2/(1-gamma)-5)*(1-gamma)^4;
    dvw0(6)=gammaV2*(gammaV2/(1-gamma)-1).*(gammaV2/(1-gamma)-2)*(gammaV2/(1-gamma)-3)*(gammaV2/(1-gamma)-4)*(gammaV2/(1-gamma)-5)*((1-gamma).*c0).^(gammaV2/(1-gamma)-6)*(1-gamma)^5;
    

end

.
.
.
Here is graph of the density that will follow when you run this program.

Amin · June 4th, 2022, 6:09 am

Friends, I will be explaining things later but in the meanwhile want to ask friends to replace the old function I posted yesterday with this new function for calculation of time step value of dt-integral.
.
.

function [cH0,cH,cVdt] = CalculateVPathStepIntegral02(c0prev,cprev,c0,c,kappa,theta,gamma,dt,Ts,cH0,cH,SeriesOrder)


    gammaV2New=1;
    [vw0,dvw0] = CalculatevgammaV2AndDerivatives(c0,gamma,gammaV2New,SeriesOrder);
    %[V2gamma0,V2gamma1] = CalculateDriftbCoeffs08A(vw0,dvw0,cmid,SeriesOrder);
    [cv0,cv] = CalculateDriftbCoeffs08A(vw0,dvw0,c,SeriesOrder);

    gammaV2New=1;
    [vw0,dvw0] = CalculatevgammaV2AndDerivatives(c0prev,gamma,gammaV2New,SeriesOrder);
    %[V2gamma0,V2gamma1] = CalculateDriftbCoeffs08A(vw0,dvw0,cmid,SeriesOrder);
    [cv0prev,cvprev] = CalculateDriftbCoeffs08A(vw0,dvw0,cprev,SeriesOrder);
    
   % b0=cv0prev*(exp(-kappa*dt)-1)-theta*(exp(-kappa*dt)-1);
   % b=cvprev*(exp(-kappa*dt)-1);

    [fH0,fH] = ConvertZCoeffsToHCoeffs(cv0,cv,SeriesOrder);    
    
    [fH0prev,fHprev] = ConvertZCoeffsToHCoeffs(cv0prev,cvprev,SeriesOrder);    
 
 
   % [bH0,bH] = ConvertZCoeffsToHCoeffs(b0,b,SeriesOrder);   

    dH0=fH0-fH0prev;
 

    %dH(1:SeriesOrder)=sign(fH(1:SeriesOrder)-bH(1:SeriesOrdery)-exp(-1*kappa*dt).*(fHprev(1:SeriesOrder))).*sqrt(abs(sign(fH(1:SeriesOrdery)-bH(1:SeriesOrder)).*(fH(1:SeriesOrdery)-bH(1:SeriesOrdery)).^2-sign(exp(-kappa*dt).*(fHprev(1:SeriesOrdery))).*(exp(-2*1*kappa*dt).*(fHprev(1:SeriesOrdery)).^2))); 
    dH(1:SeriesOrder)=sign(fH(1:SeriesOrder)-exp(-1*kappa*dt).*(fHprev(1:SeriesOrder))).*sqrt(abs(sign(fH(1:SeriesOrder)).*(fH(1:SeriesOrder)).^2-sign(exp(-kappa*dt).*(fHprev(1:SeriesOrder))).*(exp(-2*1*kappa*dt).*(fHprev(1:SeriesOrder)).^2)));    
 
    cH0(Ts)=dH0;
    cH(Ts,1:SeriesOrder)=dH(1:SeriesOrder);
 
 
     Coeff1(1:Ts)=0;
    Coeff2(1:Ts)=0;

for tt=1:Ts-1
     for ss=tt:Ts-1
         
         Coeff1(tt)=Coeff1(tt)+(exp(.5*kappa*dt*(tt))).*(exp(-2*kappa*dt*(ss-tt-1))).*dt^2;    
         for rr=ss:Ts-1
            if(rr==ss)
            
            else
            Coeff1(tt)=Coeff1(tt)+2*(exp(-1*kappa*dt*(ss-tt-1)))*(exp(-1*kappa*dt*(rr-tt-1))).*dt^2.*exp(0*kappa*dt*(tt));
            end
         end
     end
end
 
for tt=1:Ts
     for ss=tt:Ts
         
         Coeff2(tt)=Coeff2(tt)+(exp(.5*kappa*dt*(tt))).*(exp(-2*kappa*dt*(ss-tt-1))).*dt^2;    
         for rr=ss:Ts
            if(rr==ss)
            
            else
            Coeff2(tt)=Coeff2(tt)+2*(exp(-1*kappa*dt*(ss-tt-1)))*(exp(-1*kappa*dt*(rr-tt-1))).*dt^2.*exp(0*kappa*dt*(tt));
            end
         end
     end
 end
%       for tt=1:Ts-1
%       for ss=tt:Ts-1
%           
%           Coeff1(tt)=Coeff1(tt)+(exp(0*kappa*dt*(tt))).*(exp(-2*kappa*dt*(ss-tt))).*dt^2;    
%           for rr=tt:Ts-1
%              if(rr==tt)
%              
%              else
%              Coeff1(tt)=Coeff1(tt)+1*(exp(-1*kappa*dt*(ss)))*(exp(-0*kappa*dt*(tt))).*dt^2;
%              end
%           end
%       end
%       end
%  
%      for tt=1:Ts
%      for ss=tt:Ts
%          
%          Coeff2(tt)=Coeff2(tt)+(exp(0*kappa*dt*(tt))).*(exp(-2*kappa*dt*(ss-tt))).*dt^2;    
%          for rr=tt:Ts
%             if(rr==tt)
%             
%             else
%             Coeff2(tt)=Coeff2(tt)+1*(exp(-1*kappa*dt*(ss)))*(exp(-0*kappa*dt*(tt))).*dt^2;
%             end
%          end
%      end
% 
%      end

 %   Coeff0(1)=Coeff2(1);

     for tt=2:Ts
%         Coeff0(tt)=(Coeff2(tt)-Coeff2(tt-1));
     end
    
    
%     for tt1=1:tt-1
%         for ss=tt1:tt-1
%             for rr=ss:tt-1
%                 if(rr==ss)
%             
%                     Coeff1(tt1)=Coeff1(tt1)+(exp(0*kappa*dt*(tt1))).*(exp(-2*kappa*dt*(rr-ss))).*dt^2;    
%                 else
%            
%                     Coeff1(tt1)=Coeff1(tt1)+2*(exp(-2*kappa*dt*(rr-ss)))*(exp(-1*kappa*dt*(tt1))).*dt^2;
%                 end
%             end
%          
%         end
%      
%     end
%     for tt1=1:tt
%         for ss=tt1:tt
%             for rr=ss:tt
%                 if(rr==ss)
%             
%                     Coeff2(tt1)=Coeff2(tt1)+(exp(0*kappa*dt*(tt1))).*(exp(-2*kappa*dt*(rr-ss))).*dt^2;    
%                 else
%            
%                     Coeff2(tt1)=Coeff2(tt1)+2*(exp(-2*kappa*dt*(rr-ss)))*(exp(-1*kappa*dt*(tt1))).*dt^2;
%                 end
%             end
%          
%         end
%      
%     end
    
    eH10=0;
    eH20=0;
    eH1(1:SeriesOrder)=0;
    eH2(1:SeriesOrder)=0;
    for tt=2:Ts
    ts=Ts-1-tt+1;
    eH10=eH10+cH0(tt)*ts*dt;
    
   % if(tt*dt*kappa<4)
    eH1(1:SeriesOrder)=eH1(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*Coeff1(tt).*cH(tt,1:SeriesOrder).^2;   
%    eH1(1:SeriesOrder)=eH1(1:SeriesOrder)+sign(cH(tt1,1:SeriesOrder)).*Coeff0(tt1).*cH(tt1,1:SeriesOrder);   
   % else
   % %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery).*(1+exp(2*kappa*dt).*dt^2);
   % eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery).*(1+sqrt(2)/kappa*dt^2);
    %end
    end 
eH0=0;
eH2(1:SeriesOrder)=0;
for tt=1:Ts
    ts1=Ts-1-tt+1;
    ts2=Ts-tt+1;
    eH0=eH0+cH0(tt)*(ts2-ts1)*dt;
    
   % if(tt*dt*kappa<4)
    eH2(1:SeriesOrder)=eH2(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*Coeff2(tt).*cH(tt,1:SeriesOrder).^2;   
   % else
   % %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery).*(1+exp(2*kappa*dt).*dt^2);
   % eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery).*(1+sqrt(2)/kappa*dt^2);
    %end
end


%eH1
%eH2

%Ts
%str=input('Look at two numbers ---000');



eH(1:SeriesOrder)=0.0;
%eH0=0;

%eH20=0;
%eH10=0;


% for tt=1:Ts
%     ts=(Ts-tt+1)*dt;
%     eH2(1:SeriesOrder)=eH2(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*ts^2.*cH(tt,1:SeriesOrder).^2;   
% %    eH0=eH0+exp(-2.*kappa*dt).*cH0(tt)*dt;%*ts*dt;
%     eH20=eH20+cH0(tt)*ts;
%  
% %     if(tt>=1)
% %     eH2(1:SeriesOrder)=eH2(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*Coeff2(tt).*cH(tt,1:SeriesOrder).^2;   
% %     else
% %     eH2(1:SeriesOrder)=eH2(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*Coeff2(tt).*cH(tt,1:SeriesOrder).^2;       
% %     end
%  
% end
% for tt=1:Ts-1
%  %   ts=Ts-tt+1;
%     
%  %   eH0=eH0+exp(-2.*kappa*dt).*cH0(tt)*dt;%*ts*dt;
%     ts=(Ts-tt)*dt;
%     eH1(1:SeriesOrder)=eH1(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*ts^2.*cH(tt,1:SeriesOrder).^2;   
% %    eH0=eH0+exp(-2.*kappa*dt).*cH0(tt)*dt;%*ts*dt;
%     eH10=eH10+cH0(tt)*ts;
%  
% %     if(tt>=1)
% %     eH1(1:SeriesOrder)=eH1(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*Coeff1(tt).*cH(tt,1:SeriesOrder).^2;   
% %     else
% %     eH1(1:SeriesOrder)=eH1(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*Coeff1(tt).*cH(tt,1:SeriesOrder).^2;       
% %     end
%  
% end
% 

%eH0=eH20-eH10;
%eH0=dH0;
%eH0=fH0*dt;

eH(1:SeriesOrder)=sign(eH2(1:SeriesOrder)).*sqrt(abs(sign(eH2(1:SeriesOrder)).*eH2(1:SeriesOrder)))-sign(eH1(1:SeriesOrder)).*sqrt(abs(sign(eH1(1:SeriesOrder)).*eH1(1:SeriesOrder)));
%eH(1:SeriesOrder)=sign(eH2(1:SeriesOrder)).*sqrt(abs(sign(eH2(1:SeriesOrder)).*eH2(1:SeriesOrder).^2))-sign(eH1(1:SeriesOrder)).*sqrt(abs(sign(eH1(1:SeriesOrder)).*eH1(1:SeriesOrder).^2));
%eH(1:SeriesOrder)=sign(eH2(1:SeriesOrder)-eH1(1:SeriesOrder)).*sqrt(abs(sign(eH2(1:SeriesOrder)).*eH2(1:SeriesOrder).^2-sign(eH1(1:SeriesOrder)).*eH1(1:SeriesOrder).^2));

%eH(1:SeriesOrder)=sign(eH2(1:SeriesOrder)).*sqrt(abs(eH2(1:SeriesOrder)));

%eH(1:SeriesOrder)=sign(eH2(1:SeriesOrder)-eH1(1:SeriesOrder)).*sqrt(abs(eH2(1:SeriesOrder)-eH1(1:SeriesOrder)));

%eH0=eH20-eH10;
%eH0=eH20;
%eH0=cH0(tt);
 
 %We convert hermite represenation given by bH0 and bH(1:5) into 
 %dt-integral Z-series Coeffs, e0 and e(1:5)
 %eH(SeriesOrdery)=eH(SeriesOrdery)/2;
[c0vdt,cvdt] = ConvertHCoeffsToZCoeffs(eH0,eH,SeriesOrder);%  
%[c0vdt,cvdt] = ConvertHCoeffsToZCoeffs(cH0,cH,SeriesOrder);%  

 cVdt(1)=c0vdt;
 cVdt(2:6)=cvdt(1:5);

 
 %cVdt
 
 
% cVdt(1)=cv0*dt;
% cVdt(2:6)=cvprev(1:5)*dt;

% cVdt
 
 
% cH
% Ts
 
% str=input('Look at numbers');

end

Amin · June 4th, 2022, 7:51 am

Friends, please use this new improved program. Bessel drift term multiplies with variance and that multiplication also has to be taken in proper dt-integral sense. I was not doing that in previous program. This program is much improved especially for terminal time beyond two years.
.
.
.

function [] = FPESeriesCoeffsFromGaussiansWilmottNew02A_Lite_2DFPEA2()

%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334 
%or skype ahsan.amin2999

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

dt=.125/2/2/2;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.
Tt=64*5;%16*2*4;%*4*4*1;%*4;%128*2*2*2;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year; 
T=Tt*dt;
% 
% %Below I have done calculations for smaller step size at start.
% 
% %ds(1:64)=dt/16;
% %ds(65:128)=dt/4;
% if(Tt<=4)
% Ts=(Tt*16);%+(64-16);
% ds(1:Ts)=dt/16;
% end
% if(Tt>4)&&(Tt<=20)
% Ts=(64)+((Tt-4)*4);
% ds(1:64)=dt/16;
% ds(65:Ts)=dt/4;
% end
% if(Tt>20)
% Ts=(64)+((20-4)*4)+(Tt-20);    
% ds(1:64)=dt/16;
% ds(65:128)=dt/4;
% ds(129:Ts)=dt;
% end

Ts=Tt;
ds(1:Tt)=dt;

T
sum(ds(1:Ts))
Ts

%Ts=4;
str=input('Check if time is right');

OrderA=4;  %
OrderM=4;  %

if(T>=1)
dtM=1/32;%.0625/1;
TtM=T/dtM;
else
dtM=T/32;%.0625/1;
TtM=T/dtM;
end

%dtM=dt;
%TtM=Tt;

dNn=.2/2;   % Normal density subdivisions width. would change with number of subdivisions
Nn=90;  % No of normal density subdivisions

NnMid=4.0;

v00=.1250;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.950;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=1.0;%.950;   %mean reversion parameter.
theta=.1250;%mean reversion target
sigma0=.9500;%Volatility value
NoOfCumulants=6;
SeriesOrder=NoOfCumulants-1;

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

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


x0=1.00;
gammaX=.75;
gammaV=.5;
sigmaX=1.0;
thetaX=.50;
kappaX=.0;
mu1X=thetaX*kappaX;
mu2X=-1*kappaX;
beta1X=0.0;
beta2X=1.0;
q0=x0^(1-gammaX)/(1-gammaX);
alpha1X=1-gammaX;



Z(1:Nn)=(((1:Nn)-6.5/1)*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);
ZProb(2:Nn-1)=normcdf(.5*Z(2:Nn-1)+.5*Z(3:Nn),0,1)-normcdf(.5*Z(2:Nn-1)+.5*Z(1:Nn-2),0,1);


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

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


sigma11X(1:OrderA+1)=0;
mu11X(1:OrderA+1)=0;
mu22X(1:OrderA+1)=0;
sigma22X(1:OrderA+1)=0;

sigma11X(1)=1;
mu11X(1)=1;
mu22X(1)=1;
sigma22X(1)=1;

for k=1:(OrderA+1)
    if sigmaX~=0
        sigma11X(k)=sigmaX^(k-1);
    end
    if mu1X ~= 0
        mu11X(k)=mu1X^(k-1);
    end
    if mu2X ~= 0
        mu22X(k)=mu2X^(k-1);
    end
    if sigmaX~=0
        sigma22X(k)=sigmaX^(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.

Fp1X(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
%YCoeff0 and YCoeff are coefficents for original coordinates monte carlo.
%YqCoeff0 and YqCoeff are bessel/lamperti version monte carlo.

YCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
YqCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
YqCoeff0X(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;

%Pre-compute the time and power exponent values in small multi-dimensional arrays
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma); %expand y^alpha where alpha=1;
YqCoeff = ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);%expand y^alpha1 where alpha1=(1-gamma)
YqCoeff=YqCoeff/(1-gamma); %Transformed coordinates coefficients have to be 
%further divided by (1-gamma)
YqCoeffX = ItoTaylorCoeffsNew(alpha1X,beta1X,beta2X,gammaX);%expand y^alpha1 where alpha1=(1-gamma)
YqCoeffX=YqCoeffX/(1-gammaX); 
for k = 0 : (OrderA)
    for m = 0:k
        l4 = k - m + 1;
        for n = 0 : m
            l3 = m - n + 1;
            for j = 0:n
                l2 = n - j + 1;
                l1 = j + 1;
                %Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
                Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                
                Fp1X(l1,l2,l3,l4) = (alpha1X + (l1-1) * beta1X + (l2-1) * beta2X + (l3-1) * 2* gammaX + (l4-1) * gammaX ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                
                
                YCoeff0(l1,l2,l3,l4) =YCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
                YqCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
                YqCoeff0X(l1,l2,l3,l4) =YqCoeffX(l1,l2,l3,l4).*mu11X(l1).*mu22X(l2).*sigma22X(l3).*sigma11X(l4);
            end
        end
    end
end


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

ExpnOrder=5;
SeriesOrder=5;
c(1:SeriesOrder)=0;
c0=w0;
wnStart=1;
yydt(wnStart:Nn)=0.0;
tic

for tt=1:Ts
    tt

    if(tt<=1)
        [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
    cprev(1:5)=0;
    c0prev=w0;
        c0=wMu0dt+w0;
        c(1)=sigma0*sqrt(ds(tt));
        c(2:7)=0.0;
    
        w0=c0;
    elseif(tt<=3) 

    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),8);
     [b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
     
    c0prev=c0;
    cprev=c;
    c0=c0+b0;
     c(1)=b(1)+sqrt(c(1)^2+sigma0^2*ds(tt));
     c(2:SeriesOrder)=b(2:SeriesOrder);
     
     w0=c0;
    else
        
     

        
        
     [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
     [b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
     
     
     
    [wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
    
    [wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
    
    [g10,g1] = CalculateDriftbCoeffs08A(wVol0dt,dwVol0dtdw,c,SeriesOrder);
    g10=g10+sigma0*sqrt(ds(tt));
    [g20,g2] = CalculateDriftbCoeffs08A(wVol2dt,dwVol2dtdw,c,SeriesOrder);
    
cprev=c;
c0prev=c0;
  a0=c0+b0;
  a(1:SeriesOrder)=c(1:SeriesOrder)+b(1:SeriesOrder);
%     
      [a0,a] = AdvanceSeriesCumulantSolution02NewO6Params5A_Lite_0(a0,a,g10,g1,g20,g2,sigma0,gamma,ds(tt),SeriesOrder,NoOfCumulants);

 c0=a0;
 c(1:SeriesOrder)=a(1:SeriesOrder);
    
    w0=c0;
    end
    
     if(tt==1)
        cmid0=c0;
        cmid=c;
        
     else
       cmid0=(c0+c0prev)*.5;
        cmid=(c+cprev)*.5;  
        cmid0=c0;
        cmid=c;  
        
     end
    
 
     
      gammaV2=2*gammaV;
      [vw0,dvw0] = CalculatevgammaV2AndDerivatives(c0prev,gamma,gammaV2,SeriesOrder);
%      %[V2gamma0,V2gamma1] = CalculateDriftbCoeffs08A(vw0,dvw0,cmid,SeriesOrder);
      [V2gamma0,V2gamma1] = CalculateDriftbCoeffs08A(vw0,dvw0,cprev,SeriesOrder);
      
      V2gamma(1)=V2gamma0;
      V2gamma(2:6)=V2gamma1(1:5);
%      
     
     
    if(tt==1)

        cq(1:6,1:6)=0;    
        cq(1,1)=x0^(1-gammaX)/(1-gammaX);
        
        
    end
    %    cprev=c;
    
    q0=cq(1,1);

    
    if(tt==1)
        
        gammaV2=2*gammaV;
     [vw0,dvw0] = CalculatevgammaV2AndDerivatives(c0prev,gamma,gammaV2,SeriesOrder);
     %[V2gamma0,V2gamma1] = CalculateDriftbCoeffs08A(vw0,dvw0,cmid,SeriesOrder);
     [ch0,ch] = CalculateDriftbCoeffs08A(vw0,dvw0,cprev,SeriesOrder);
     [ch0,ch] = ConvertZCoeffsToHCoeffs(ch0,ch,SeriesOrder);    
    
     cH0(tt)=ch0;
     cH(tt,1:5)=ch(1:5);
     
%     tt
     
     
    else
    [cH0,cH,cVdt] = CalculateVPathStepIntegral02(c0prev,cprev,c0,c,kappa,theta,gamma,dt,tt,cH0,cH,SeriesOrder);

    
    end
    
    
    

    DoubleExpnOrder=ExpnOrder*2;
    
    [qMu0dt,dqMu0dtdq,qMu1dt,dqMu1dtdq] = BesselDriftAndDerivatives08A0(q0,mu1X,mu2X,beta1X,beta2X,gammaX,DoubleExpnOrder);
    [q2Mu0dt,dq2Mu0dtdq,q2Mu1dt,dq2Mu1dtdq] = BesselDriftAndDerivatives08Aq(q0,sigmaX,gammaX,DoubleExpnOrder);
 
    [Muq0] = CalculateDriftbCoeffs08A2Dim02(qMu0dt,dqMu0dtdq,cq,SeriesOrder);
    [Muq2] = CalculateDriftbCoeffs08A2Dim02(q2Mu0dt,dq2Mu0dtdq,cq,SeriesOrder)/dt;
    [DMuq0] = CalculateDriftbCoeffs08A2Dim02(qMu1dt,dqMu1dtdq,cq,SeriesOrder);
    [DMuq2] = CalculateDriftbCoeffs08A2Dim02(q2Mu1dt,dq2Mu1dtdq,cq,SeriesOrder)/dt;
    
 
    
    if(tt<=4)
    [Muq2] =SeriesProduct2D1D(Muq2,V2gamma)*dt;
    [DMuq2] =SeriesProduct2D1D(DMuq2,V2gamma)*dt;
    else
    [Muq2] =SeriesProduct2D1D(Muq2,cVdt);
    [DMuq2] =SeriesProduct2D1D(DMuq2,cVdt);
    end
    
    
    Muq=Muq0+Muq2;%+Muq01*dt/2;%-tt*Muq02*dt/2+tt*Muq01*dt/2;
    DMuq=DMuq0+DMuq2;
    
    if(tt==1)
        cq=cq+Muq*dt+SeriesProduct2D(Muq,DMuq)*dt^2/2;
        cq(2,1)=cq(2,1)+sigmaX*v00^gammaV*sqrt(dt);
    elseif(tt<=3)
%        cq
        [D1cq] = Series2DZ1Derivative(cq);
%        D1cq
        [D1cqInv] = Series2DReciprocal02(D1cq);
%        D1cqInv
        [D1cqInvZ] = MultiplySeries2DWithZ1(D1cqInv);
%        D1cqInvZ
        fq=Muq+.5*sigmaX^2*SeriesProduct2D1D(D1cqInvZ,V2gamma);
%        fq*dt
        Dfq=DMuq;
        cq=cq+fq*dt+SeriesProduct2D(fq,Dfq)*dt^2/2;
        
        
        SeriesProduct2D1D(D1cqInvZ,V2gamma)
%        str=input('Look at numbers-1111');
%        D1cqInvZ
%        str=input('Look at numbers-2222');
    else
        [D1cq] = Series2DZ1Derivative(cq);
    
        [D1cqInv] = Series2DReciprocal02(D1cq);
        [D1cqInvZ] = MultiplySeries2DWithZ1(D1cqInv);
        [D2cq] = Series2DZ1Derivative(D1cq);
        [D1cqInvSqr]=SeriesProduct2D(D1cqInv,D1cqInv);
        [D1cqInvCub]=SeriesProduct2D(D1cqInvSqr,D1cqInv);
        [D2cqD1cqInvSqr]=SeriesProduct2D(D2cq,D1cqInvSqr);
        [D2cqD1cqInvCub]=SeriesProduct2D(D2cq,D1cqInvCub);
        D2cqT2=D1cqInvSqr-MultiplySeries2DWithZ1(D2cqD1cqInvCub);
        Dfq=DMuq*dt+.5*sigmaX^2*SeriesProduct2D1D(D2cqT2,cVdt);
        D2cqT1=D1cqInvZ+D2cqD1cqInvSqr;
        

        fq=Muq*dt+.5*sigmaX^2*SeriesProduct2D1D(D2cqT1,cVdt);    
 
       cq=cq+(fq+SeriesProduct2D(fq,Dfq)*dt^0/2);%.*cprev(1)/c(1);
        
    end
 
 


end

Z1=Z;

if(tt*dt>1)
    cq(:,5)=cq(:,5)/4;
    cq(5,:)=cq(5,:)/4;
end

[b] = EvaluateSeriesForZ2(cq,Z1,Nn);
% Mm=301;
% MmMid=151;
% Q(MmMid+1:Mm)=cq(1,1)+cq(2,1).*.1.*(1:MmMid-1);
% Q(MmMid)=cq(1,1);
% %Q(1:70)=cq(1,1)-cq(2,1).*.25.*(70:-1:1);
% for mm=1:MmMid-1
% 
%     Q(MmMid-mm)=cq(1,1)-cq(2,1).*.1.*mm;
% end
% 
% 
% Qa=Q-cq(2,1).*.1/2;
% Qb=Q+cq(2,1).*.1/2;



Mm=601;
MmMid=301;
dMm=.05;
Q(MmMid+1:Mm)=cq(1,1)+cq(2,1).*dMm.*(1:MmMid-1);
Q(MmMid)=cq(1,1);
%Q(1:70)=cq(1,1)-cq(2,1).*.25.*(70:-1:1);
for mm=1:MmMid-1

    Q(MmMid-mm)=cq(1,1)-cq(2,1).*dMm.*mm;
end


Qa=Q-cq(2,1).*dMm/2;
Qb=Q+cq(2,1).*dMm/2;






Z1Prob(1:Mm)=0;
for nn=1:Nn
    bq(1:6)=b(1:6,nn);

    
    [Z0Prob] = CalculateProbabilitySeriesInv(Qa,Qb,bq);
    
%    Z0Prob
    
%str=input('Look at numbers-tttt')
    
    
    Z1Prob(1:Mm)=Z1Prob(1:Mm)+Z0Prob(1:Mm).*ZProb(nn);
    
end

pQ=Z1Prob/(cq(2,1).*dMm);
    
xx(1:Mm)=((1-gammaX).*Q(1:Mm)).^(1/(1-gammaX));

pxx(1:Mm)=pQ(1:Mm).*xx(1:Mm).^(-gammaX);
    
    
Mean=sum(xx(1:Mm).*Z1Prob(1:Mm));
Mean
  w(1:Nn)=c0;
  for nn=1:SeriesOrder
      w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^nn;
  end
  
  
 % q(1:Nn)=cq0;
 % for nn=1:SeriesOrder
 %     q(1:Nn)=q(1:Nn)+cq(nn)*Z(1:Nn).^nn;
 % end
  
  
  
%  %Flag=0;
%  %for nn=ceil(Nn/2)-1:-1:1
%  %    if((w(nn)<0)&&(Flag==0))
%  %        wnStart=nn+1;
%  %        Flag=1;
%  %    end
%  %end
%  
%  
%  %c0
%  %c
%  
%  
%  %str=input('Look at numbers')
  yy0(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
  
 % xx0(wnStart:Nn)=((1-gammaX).*q(wnStart:Nn)).^(1/(1-gammaX));
%  
%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%   yw0=((1-gamma).*c0).^(1/(1-gamma));
%   dyw0(1)=((1-gamma).*c0).^(1/(1-gamma)-1);
%   dyw0(2)=(1/(1-gamma)-1).*((1-gamma).*c0).^(1/(1-gamma)-2)*(1-gamma);
%   dyw0(3)=(1/(1-gamma)-1).*(1/(1-gamma)-2).*((1-gamma).*c0).^(1/(1-gamma)-3)*(1-gamma)^2;
%   dyw0(4)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*((1-gamma).*c0).^(1/(1-gamma)-4)*(1-gamma)^3;
%   dyw0(5)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*((1-gamma).*c0).^(1/(1-gamma)-5)*(1-gamma)^4;
%   dyw0(6)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*((1-gamma).*c0).^(1/(1-gamma)-6)*(1-gamma)^5;
%   dyw0(7)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*(1/(1-gamma)-6)*((1-gamma).*c0).^(1/(1-gamma)-7)*(1-gamma)^6;
%   dyw0(8)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*(1/(1-gamma)-6)*(1/(1-gamma)-7)*((1-gamma).*c0).^(1/(1-gamma)-8)*(1-gamma)^7;
% %  
% % c(6)=0;
% % c(7)=0;
% % c(8)=0;
% % 
%      [y10,y1] = CalculateDriftbCoeffs08A(yw0,dyw0,c,SeriesOrder);
%      [y20,y2] = CalculateDriftbCoeffs08A(yw0,dyw0,c,5);
% %  
%  
%  yy1(1:Nn)=y10;
%  %y(1:SeriesOrder)=c(1:SeriesOrder);
%  for nn=1:SeriesOrder
%      yy1(1:Nn)=yy1(1:Nn)+y1(nn)*Z(1:Nn).^nn;
%  end
%  
%  yy2(1:Nn)=y20;
%  %y(1:SeriesOrder)=c(1:SeriesOrder);
%  for nn=1:5
%      yy2(1:Nn)=yy2(1:Nn)+y2(nn)*Z(1:Nn).^nn;
%  end
%  plot(Z(1:Nn),yy0(1:Nn),'r',Z(1:Nn),yy1(1:Nn),'b',Z(1:Nn),yy2(1:Nn),'k');
% title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f', x0,theta,kappa,gamma,sigma0,T,dt));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
%  
% legend({'SDE Variable Converted on Grid','SDE Variable From Eight Cumulants','SDE Variable From Six Cumulants'},'Location','northeast')
%  
% %str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
% 
%  
%  str=input('red line is SDE from Ito-Hermite method.');
%  
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
 yy=yy0;
% xx=xx0;
Dfyy(wnStart:Nn)=0;
%Dfyy1(wnStart:Nn)=0;
%Dfyy2(wnStart:Nn)=0;
%Dfxx(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    
    Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%    Dfyy1(nn) = (yy1(nn + 1) - yy1(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%    Dfyy2(nn) = (yy2(nn + 1) - yy2(nn - 1))/(Z(nn + 1) - Z(nn - 1));
 %   Dfxx(nn) = (xx(nn + 1) - xx(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    
    %Change of variable derivative for densities
end
%Dfyy(Nn)=Dfyy(Nn-1);
%Dfyy(1)=Dfyy(2);
%Dfyydt(Nn)=Dfyydt(Nn-1);
%Dfyydt(1)=Dfyydt(2);


pyy(1:Nn)=0;
%pyy1(1:Nn)=0;
%pyy2(1:Nn)=0;
%pxx(1:Nn)=0;
for nn = wnStart+1:Nn-1
    
    pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
%    pyy1(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy1(nn));
%    pyy2(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy2(nn));
 %   pxx(nn) = (normpdf(Z(nn),0, 1))/abs(Dfxx(nn));
end


toc

ItoHermiteMeanVar=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
ItoHermiteMeanAsset=Mean
disp('true Mean only applicable to standard SV mean reverting type models otherwise disregard');
%TrueMean=theta+(v0-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average
%yy0


theta1=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%rng(29079137, 'twister')
rng(15898837, 'twister')
paths=100000;
V(1:paths)=v00;  %Original process monte carlo.
X=0.0;
X(1:paths)=x0;
alpha1=0;
alpha2=1;
a=mu1X;
b=mu2X;
rho=0;
sigma1=sigmaX;
gammaV=.5;
Random1(1:paths)=0;
Random2(1:paths)=0;

for ttM=1:TtM
Random1=randn(size(Random1));
Random2=randn(size(Random2));



X(1:paths)=X(1:paths)+ ...
    (a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dtM + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) * sqrt(dtM) + ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*sqrt(dtM) + ...
    (a*alpha1* X(1:paths).^(alpha1-1)+b*alpha2* X(1:paths).^(alpha2-1)).* ...
    (((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dtM^2/2)+ ...
    (sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) *(1-1/sqrt(3)).*dtM^1.5+ ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*(1-1/sqrt(3)).*dtM^1.5))+ ...
    .5*(a*alpha1*(alpha1-1)* X(1:paths).^(alpha1-2)+b*alpha2*(alpha2-1).* X(1:paths).^(alpha2-2)).* ...
    ( sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX)) *dtM^2/2 + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
    ((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random1(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*(Random1(1:paths).^2-1) * dtM/2 + ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths)*dtM/2)+ ...
    .5*sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
    (sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random1(1:paths).*1/sqrt(3).*dtM^1.5)+ ...
    sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
    ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random1(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
    sigma0*V(1:paths).^gamma.*Random1(1:paths).*Random2(1:paths)*dtM/2)+ ...
    .5*sqrt(1-rho^2)* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
    (sigma0^2*V(1:paths).^(2*gamma).*Random1(1:paths)*1/sqrt(3)*dtM^1.5)+ ...
    sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
    rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random1(1:paths)*1/sqrt(3)*dtM^1.5+ ...
    rho* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
    ((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random2(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths) * dtM/2 + ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*(Random2(1:paths).^2-1)*dtM/2)+ ...
    .5*rho* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
    (sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random2(1:paths).*1/sqrt(3).*dtM^1.5)+ ...
    rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
    ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
    sigma0*V(1:paths).^gamma.*(Random2(1:paths).^2-1)*dtM/2)+ ...
    .5*rho* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
    sigma0^2.*V(1:paths).^(2*gamma).*Random2(1:paths) * 1/sqrt(3).* dtM^1.5+ ...
    rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
    rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random2(1:paths)*1/sqrt(3)*dtM^1.5;
    
VBefore=V;

    V(1:paths)=V(1:paths)+ ...
        (mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dtM + ...
        sigma0*V(1:paths).^gamma .*Random2(1:paths)*sqrt(dtM) + ...
        (mu1.*beta1*V(1:paths).^(beta1-1) + mu2.*beta2.*V(1:paths).^(beta2-1)).* ...
        ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dtM^2/2 + ...
        sigma0*V(1:paths).^gamma .*Random2(1:paths)*(1-1/sqrt(3))*dtM^1.5) + ...
        .5*(mu1.*beta1.*(beta1-1).*V(1:paths).^(beta1-2) + mu2.*beta2.*(beta2-1).*V(1:paths).^(beta2-2)).* ...    
        sigma0^2.*V(1:paths).^(2*gamma).*dtM^2/2 + ...
        sigma0*gamma*V(1:paths).^(gamma-1) .* ...
        ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths).*1/sqrt(3)*dtM^1.5 + ...
        sigma0.*V(1:paths).^gamma .*(Random2(1:paths).^2-1)*dtM/2) + ...
        .5*sigma0*gamma*(gamma-1).*V(1:paths).^(gamma-2) .* ...
        sigma0^2.*V(1:paths).^(2*gamma) .*Random2(1:paths).*1/sqrt(3)*dtM^1.5;

%Vm=.5*VBefore+.5*V;
Vm=V;
    
%    X(1:paths)=X(1:paths)+ ...
%    sigma1* Vm(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) * sqrt(dtM) + ...
%    sigma1* Vm(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
%    (sigma1* Vm(1:paths).^gammaV.* X(1:paths).^gammaX .*(Random1(1:paths).^2-1) * dtM/2)+ ...
%    .5*sigma1* Vm(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
%    (sigma1^2* Vm(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random1(1:paths).*1/sqrt(3).*dtM^1.5);%+ ...

    
    
    
end

%SVolMeanAnalytic=thetaV+(V0-thetaV)*exp(-kappaV*dt*Tt)
SVolMeanAnalytic=theta+(v00-theta)*exp(-kappa*dt*Tt)
SVolMeanMC=sum(V(1:paths))/paths
AssetMeanAnalytic=x0
AssetMeanMC=sum(X(1:paths))/paths

% theta
% v00
% kappa
% dt
% Tt
% T


%MeanX

disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(v00-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average


MaxCutOff=30;
NoOfBins=2*300;%round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too 
[XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti_NEW(X,paths,NoOfBins,MaxCutOff );
plot(xx(wnStart+1:Mm-1),pxx(wnStart+1:Mm-1),'r',IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g');
 %plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');
 
title(sprintf('x0 = %.2f,thetaX=%.2f,kappaX=%.2f,gammaX=%.2f,sigmaX=%.2f,v0 =%.2f,kappa=%.2f,theta=%.2f,gamma=%.2f,sigma0=%.2f,T=%.2f,dt=%.3f,M=%.4f',x0,thetaX,kappaX,gammaX,sigmaX,v00,kappa,theta,gamma,sigma0,T,dt,ItoHermiteMeanAsset));
%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
 
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');


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

NoOfBins=round(2*500*gamma^2*4*sigma0/sqrt(SVolMeanMC)/(1+kappa))/10;
[VDensity,IndexOutV,IndexMaxV] = MakeDensityFromSimulation_Infiniti_NEW(V,paths,NoOfBins,MaxCutOff );

%[YdtDensity,IndexOutYdt,IndexMaxYdt] = MakeDensityFromSimulation_Infiniti_NEW(YYdt,paths,NoOfBins,MaxCutOff );
plot(yy(wnStart+1:Nn-1),pyy(wnStart+1:Nn-1),'r',IndexOutV(1:IndexMaxV),VDensity(1:IndexMaxV),'g');
title(sprintf('v0 =%.2f,kappa=%.2f,theta=%.2f,gamma=%.2f,sigma0=%.2f,T=%.2f,dt=%.3f,M=%.4f',v00,kappa,theta,gamma,sigma0,T,dt,ItoHermiteMeanVar));



end

.
.
.
Here is graph of a 5-yr asset from pair of SV SDEs, parameters are on the graph title.
.
.
.

Amin · June 5th, 2022, 8:49 am

Friends, you can freely use V0!=theta (Starting vol V0 not equal to theta) for SV SDE by using the last set of functions I posted yesterday.
But dt-integral step value function needs improvement.
After overcoming these problems, I want to write a more general program where friends could specify expansion order of 2D Z-series in both volatility and asset direction. I also want to improve the part where I have used series inversion to calculate asset density from 2D density. There will also be a lot of small improvements that I want to continue to do during next week.

Amin · June 6th, 2022, 1:57 pm

Friends, today I continued to work on improvement of dt-integral. I think my analytics are sound but still my results for dt-integral for mean-reverting SDEs start to slightly diverge from monte carlo results of dt-integral after 3-4 years. with kappa=1, the difference is small but continues to increase with time and after six years, they are not close to each other anymore. Tomorrow, I will outline my analytics and related solution for dt-integral for mean-reverting SDEs.
Meanwhile mind control torture continues. I have been sleeping in the car garage (porch) since the summers started. Yesterday when I went to sleep in the garage, it was very hot and I decided to sleep inside and turned on air-conditioner in my room. There was no smell in the air but I started to feel very miserable. There were some gases in the air-conditioner. I would not be able to sleep but I was not in my proper senses either. I knew it was a bad decision to sleep inside and wanted to go to porch again but I had inhaled enough gas and therefore I had no energy left to get up and go outside and I continued to simply lie devoid of any energy. After a very miserable one and a half hour, I somehow mustered enough energy to turn off the air-conditioner and went to sleep in the porch again. Mind control guy was saying that they would not go soft on me as every time they become soft with me and let me work with ease, it results in a nightmare for them. I want to request friends to ask the mind control agencies to become better and let me use air-conditioners now as it is already very very hot at this time of the year.

Amin · June 7th, 2022, 11:35 am

Friends for a context to post below and in order to get oriented with the problem, please read the following posts first
https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1200#p867718
https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1215#p867940
https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1215#p867941
https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1290#p868483
https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1290#p868484
https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1290#p868487
Posts related to above posts and also recent posts

Amin · June 7th, 2022, 12:19 pm

We have a general two SDE problem given as

[$] dV(t)= \mu_1 \, V(t)^{\beta_1} \, dt +\, \mu_2 \, V(t)^{\beta_2} \, dt \, +\, \sigma \, V(t)^{\gamma} \, dZ_1(t)[$]

our asset SDE is given as

[$] dX(t)= \mu_{x1} \, X(t)^{\beta_{x1}} \, dt +\, \mu_{x2} \, X(t)^{\beta_{x2}} \, dt \, +\, \sigma_x \, V(t)^{\gamma_v} \, X(t)^{\gamma_x} \, dZ_2(t)[$]

We do not deal with above SDEs in original coordinates and rather convert them to Bessel coordinates.
Here I will take a break and tell friends that whenever I tried solving even 1D SDEs analytically in original coordinates, I failed however I could solve 1D SDEs in Bessel coordinates very easily. Only know when I was working on this 2D SV project that I realized that original coordinates in 1D require solution of proper stochastic dt-integral in variance as opposed to constant squared volatility in Bessel coordinates. At that time I did not have a proper understanding of the problem and I would try by calculating dt-integral merely by a multiplication with delta t and it never worked. Working on this project gave me understanding that 1D SDEs can also be solved in original coordinates in various ways but we have to be careful with stochastic dt-integral.

For the current 2D SV project, we remain in Bessel framework for both asset and stochastic volatility SDEs.

For the stochastic volatility SDE we have the new transformed SDE in bessel coordinates given as

[$] dw(t)= \mu_w(t) dt +\, \sigma \, dZ_1(t)[$]

where [$]\mu_w \, = \mu_1 \, {((1-\gamma)w(t))}^{\beta_1-\gamma} \, dt +\, \mu_2 \, {((1-\gamma)w(t))}^{\beta_2-\gamma} \, dt \, - \, \frac{ .5 \,{\sigma}^2 \, \gamma}{((1-\gamma)w(t))} \, dt [$]

while our asset SDE in Bessel coordinates becomes

[$] dq(t)= \mu_q(t) dt +\, \sigma_X \, V(t)^{\gamma_v} \, dZ_2(t)[$]

where
[$]\mu_q \, = \mu_{x1} \, {((1-\gamma_x)q(t))}^{\beta_{x1}-\gamma_x} \, dt +\, \mu_{x2} \, {((1-\gamma_x)q(t))}^{\beta_{x2}-\gamma_x} \, dt \, - \, \frac{ .5 \,{\sigma}^2 \,\gamma_x \,{V(t)}^{2 \, \gamma_v} }{((1-\gamma_x)q(t))} \, dt [$]

Breakthrough in the theory of stochastic differential equations and their simulation

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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