function [] = FPESeriesCoeffsFromGaussiansWilmottNew02A_Lite_PIntegral01()

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

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

dt=.125/2/2/2;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.
Tt=32*2*2;%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;
Ts=Tt;
ds(1:Ts)=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*8);%+(64-16);
% ds(1:Ts)=dt/8;
% end
% if(Tt>4)&&(Tt<=20)
% Ts=(32)+((Tt-4)*2);
% ds(1:32)=dt/8;
% ds(33:Ts)=dt/2;
% end
% if(Tt>20)
% Ts=(32)+((20-4)*2)+(Tt-20);    
% ds(1:32)=dt/8;
% ds(33:64)=dt/2;
% ds(65:Ts)=dt;
% end
% T
% sum(ds(1:Ts))
% Ts
%str=input('Check if time is right');

OrderA=4;  %
OrderM=4;  %

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

dtM=dt;
TtM=Tt;



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

NnMid=4.0;

x0=1.0;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.75;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=0.0;%.950;   %mean reversion parameter.
theta=1.00;%mean reversion target
sigma0=.650;%Volatility value


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

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

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
%Ft(1:TtM+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0; %General time powers on hermite polynomials
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers on coefficients of hermite polynomials.
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers for bessel transformed coordinates.

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

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

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

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

ExpnOrder=5;
SeriesOrder=5;

wnStart=1;
tic

for tt=1:Ts
    tt

    if(tt==1)
        %dt=ds(1);
        %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder);
 %       [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),ExpnOrder);
        
        [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
    %    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
        
      %  dwMu0dtdw
      %  dwMu0dtdwb
      %  wMu0dt
      %  wMu0dtb
      %  str=input('Look at drift');
        c0=wMu0dt+w0;
        c(1)=sigma0*sqrt(ds(tt));
        c(2:10)=0.0;
        w0=c0;
    
    elseif(tt<=6) 

    %dt=ds(tt);    
        
    %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
     %[b0,b] = CalculateDriftbCoeffs08O4(wMu0dt,dwMu0dtdw,c);
     [b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
     c0=c0+b0;
     c(1)=b(1)+sqrt(c(1)^2+sigma0^2*ds(tt));
     c(2)=b(2);
     c(3)=b(3);
     c(4)=b(4);
     c(5)=b(5);
     
     w0=c0;
            
    else
        
      
        
        
        
        
    %below wMu0dt is drift and dwMu0dtdw are its derivatives.
   %  [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
     [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
     [b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
     
     %size(b)
     %ba(1:6)=b(1:6);
     
     
    [wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
    
    [wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
    
    %[wVol3dt,dwVol3dtdw] = BesselVolAndDerivativesH3(w0,YqCoeff0,Fp1,gamma,ds(tt),8);

    [g10,g1] = CalculateDriftbCoeffs08A(wVol0dt,dwVol0dtdw,c,SeriesOrder);
    
    [g20,g2] = CalculateDriftbCoeffs08A(wVol2dt,dwVol2dtdw,c,SeriesOrder);
 

 a0=c0+b0;
 
 a(1:5)=c(1:5)+b(1:5);
    
     [a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A_Lite(a0,a,g10,g1,g20,g2,sigma0,ds(tt));
     %[a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A(a0,a,b0,ba,g0,g,h0,h,sigma0,ds(tt));
     
     
 c0=a0;
 c(1)=a1;
 c(2)=a2;
 c(3)=a3;
 c(4)=a4;
 c(5)=a5;
    
    w0=c0;
    
    end
    
    
  %Below Calculate the Z-series expansion for SDE variable in original coordinates 
  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;

 if(tt>1)
    fH0prev=fH0;
    fHprev=fH;
 end

 %Below f0 and f define the Z_series in original coordinates
 [f0,f] = CalculateDriftbCoeffs08A(yw0,dyw0,c,SeriesOrder);
 %fH0 and fH(1:5) define the coefficients of hermite polynomial 
 %representation of SDE variable in original coordinates.  
 [fH0,fH] = ConvertZCoeffsToHCoeffs(f0,f,SeriesOrder);    

%Below calculate the Hermite representation of orthogonal increments in the
%original coordinates variable at each time step by squared subtraction of
%coefficients of hermite polynomials.
%Here cH0 and cH define the hermite polynomial coordinates of orthogonal
 %increments at each simulation step. The take double array index in 
 %simulation level, tt and hermite polynomial count (1:5).  
if(tt==1) 
 cH0=fH0;
 cH(tt,1:SeriesOrder)=fH(1:SeriesOrder);

end
if(tt>1)       
 dH0=fH0-fH0prev;
    
 dH(1:SeriesOrder)=sign(fH(1:SeriesOrder)-fHprev(1:SeriesOrder)).*sqrt(abs(sign(fH(1:SeriesOrder)).*fH(1:SeriesOrder).^2-sign(fHprev(1:SeriesOrder)).*fHprev(1:SeriesOrder).^2));    
 
 cH0(tt)=dH0;
 cH(tt,1:SeriesOrder)=dH(1:SeriesOrder);
end



end
wnStart=1;


yy0=((1-gamma).*w0).^(1/(1-gamma));

 w(1:Nn)=c0;
 for nn=1:ExpnOrder
     w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^nn;
 end
 
  %Below calculate the hermite representation of dt-integral in original
 %coordinates
 %Here cH0 and cH define the hermite polynomial coordinates of orthogonal
 %increments at each simulation step and bH0 and bH define the hermite
 %representation of the dt-integral after continued squared addition of 
 %scaled/weighted values of hermite coefficients of orthogonal increments.
 
 
 bH0=cH0(1)*dt*Ts;
 bH(1:SeriesOrder)=cH(1,1:SeriesOrder)*dt*Ts;
 
 for tt=2:Ts
    ts=Ts-tt+1;
    bH0=bH0+cH0(tt)*dt*ts;
    bH(1:SeriesOrder)=sign(bH(1:SeriesOrder)+dt*(ts).*cH(tt,1:SeriesOrder)).*sqrt(abs(sign(bH(1:SeriesOrder)).*bH(1:SeriesOrder).^2+sign(dt*(ts).*cH(tt,1:SeriesOrder)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrder).^2));   
 end
 
 %We convert hermite represenation given by bH0 and bH(1:5) into 
 %dt-integral Z-series Coeffs, e0 and e(1:5)
 
 [e0,e] = ConvertHCoeffsToZCoeffs(bH0,bH,SeriesOrder);
  yydt(1:Nn)=e0;
 for nn=1:ExpnOrder
     yydt(1:Nn)=yydt(1:Nn)+e(nn)*Z(1:Nn).^nn;
 end

 
 %c0
 %c
 
 
 %str=input('Look at numbers')
 yy0(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));

 yy=yy0;
Dfyy(wnStart:Nn)=0;
Dfyydt(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    
    Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    Dfyydt(nn) = (yydt(nn + 1) - yydt(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;
pyydt(1:Nn)=0;
for nn = wnStart+1:Nn-1
    
    pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
    pyydt(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyydt(nn));
end

toc

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


theta1=1;
rng(29079137, 'twister')
paths=50000;
YY(1:paths)=x0;  %Original process monte carlo.
YYdt(1:paths)=0;   

Random1(1:paths)=0;
for tt=1:TtM
    
    
    Random1=randn(size(Random1));
    HermiteP1(1,1:paths)=1;
    HermiteP1(2,1:paths)=Random1(1:paths);
    HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
    HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
    HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;


    YYdt(1:paths)=YYdt(1:paths)+.5*YY(1:paths)*dtM;   

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

    
end




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


MaxCutOff=30;
NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too 
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
plot(yy(wnStart+1:Nn-1),pyy(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 %plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');
 
title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
 
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');

NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa))/2;
[YdtDensity,IndexOutYdt,IndexMaxYdt] = MakeDensityFromSimulation_Infiniti_NEW(YYdt,paths,NoOfBins,MaxCutOff );
plot(yydt(wnStart+1:Nn-1),pyydt(wnStart+1:Nn-1),'r',IndexOutYdt(1:IndexMaxYdt),YdtDensity(1:IndexMaxYdt),'g');

title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
legend({'Ito-Hermite Path Integral Density','Monte Carlo Path Integral Density'},'Location','northeast')
 
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');



end

function [b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,a,SeriesOrder)



b0=wMu0dt;

b(1)=a(1) *dwMu0dtdw(1);

b(2)=1/2*(2* a(2) *dwMu0dtdw(1)+a(1)^2 *dwMu0dtdw(2));

b(3)=1/6*(6 *a(3) *dwMu0dtdw(1)+6 *a(1)* a(2) *dwMu0dtdw(2)+a(1)^3 *dwMu0dtdw(3));

b(4)=1/24*(24* a(4) *dwMu0dtdw(1)+12* a(2)^2 *dwMu0dtdw(2)+24 *a(1)* a(3) *dwMu0dtdw(2)+ ...
    12* a(1)^2* a(2) *dwMu0dtdw(3)+a(1)^4 *dwMu0dtdw(4));

if(SeriesOrder>=5)
b(5)=1/120*(120* a(5) *dwMu0dtdw(1)+120* a(2) *a(3) *dwMu0dtdw(2)+120* a(1)* a(4) *dwMu0dtdw(2)+ ...
    60 *a(1)* a(2)^2 *dwMu0dtdw(3)+60* a(1)^2* a(3) *dwMu0dtdw(3)+ ...
    20 *a(1)^3* a(2) *dwMu0dtdw(4)+a(1)^5 *dwMu0dtdw(5));
end
if(SeriesOrder>=6)
b(6)=1/720*(720* a(6) *dwMu0dtdw(1)+360* a(3)^2 *dwMu0dtdw(2)+720* a(2)* a(4) *dwMu0dtdw(2)+720* a(1)* a(5) *dwMu0dtdw(2)+ ...
    120* a(2)^3 *dwMu0dtdw(3)+720* a(1)* a(2)* a(3) *dwMu0dtdw(3)+360* a(1)^2* a(4) *dwMu0dtdw(3)+ ...
    180* a(1)^2* a(2)^2 *dwMu0dtdw(4)+120* a(1)^3* a(3) *dwMu0dtdw(4)+ ...
    30* a(1)^4 *a(2) *dwMu0dtdw(5)+a(1)^6 *dwMu0dtdw(6));
end


if(SeriesOrder>=7)
  b(7)=1/(720*7)*(5040* a(7)*dwMu0dtdw(1)+5040* a(3)* a(4)*dwMu0dtdw(2)+5040* a(2)* a(5) *dwMu0dtdw(2)+5040 *a(1)* a(6) *dwMu0dtdw(2)+ ...
      2520* a(2)^2* a(3) *dwMu0dtdw(3)+2520* a(1)* a(3)^2 *dwMu0dtdw(3)+5040 *a(1)* a(2)* a(4) *dwMu0dtdw(3)+2520* a(1)^2 *a(5) *dwMu0dtdw(3)+ ...
      840* a(1)* a(2)^3 *dwMu0dtdw(4)+2520* a(1)^2* a(2)* a(3) *dwMu0dtdw(4)+840* a(1)^3* a(4) *dwMu0dtdw(4)+ ...
      420* a(1)^3* a(2)^2 *dwMu0dtdw(5)+210* a(1)^4* a(3) *dwMu0dtdw(5)+ ...
      42* a(1)^5* a(2) *dwMu0dtdw(6)+a(1)^7 *dwMu0dtdw(7));
  
end
  
if(SeriesOrder>=8)
  b(8)=1/(720*7*8)*(40320* a(8) *dwMu0dtdw(1)+20160* a(4)^2 *dwMu0dtdw(2)+40320* a(3)* a(5) *dwMu0dtdw(2)+40320* a(2) *a(6) *dwMu0dtdw(2)+ ...
      40320* a(1) *a(7) *dwMu0dtdw(2)+ ...
      20160* a(2)* a(3)^2 *dwMu0dtdw(3)+20160* a(2)^2* a(4) *dwMu0dtdw(3)+40320* a(1)* a(3)* a(4) *dwMu0dtdw(3)+40320* a(1)* a(2)* a(5) *dwMu0dtdw(3)+ ...
      20160* a(1)^2* a(6) *dwMu0dtdw(3)+ ...
      1680* a(2)^4 *dwMu0dtdw(4)+20160* a(1)* a(2)^2* a(3) *dwMu0dtdw(4)+10080* a(1)^2* a(3)^2 *dwMu0dtdw(4)+20160* a(1)^2* a(2) *a(4) *dwMu0dtdw(4)+ ...
      6720* a(1)^3* a(5) *dwMu0dtdw(4)+ ...
      3360* a(1)^2* a(2)^3 *dwMu0dtdw(5)+6720* a(1)^3 *a(2)* a(3) *dwMu0dtdw(5)+1680* a(1)^4* a(4) *dwMu0dtdw(5)+ ...
      840* a(1)^4* a(2)^2 *dwMu0dtdw(6)+336* a(1)^5* a(3) *dwMu0dtdw(6)+ ...
      56* a(1)^6* a(2) *dwMu0dtdw(7)+a(1)^8 *dwMu0dtdw(8));
end
%b(9:10)=0.0; 

%Bound=abs(b(1)/b0)*gamma/(1-gamma);
%Bound=1/abs(b0/b(1))*gamma;%/(1-gamma);

%for nn=1:9
%   
%    if(abs(b(nn+1))>abs(b(nn)*Bound))
%  %      b(nn+1)=sign(b(nn+1))*abs(b(nn))*Bound;    
%    end
%end

%[b0,b] = ApplyBoundsOnSeriesCoeffs(b0,b);


 
 % 
% 
% 
% b(9)=1/(720*7*8*9)*(362880* a(9) *dwMu0dtdw(1)+ ...
%     362880* a(4)* a(5) *dwMu0dtdw(2)+ 362880* a(3)* a(6) *dwMu0dtdw(2) +362880* a(2)* a(7) *dwMu0dtdw(2) +362880* a(1) *a(8) *dwMu0dtdw(2) + ...
%     60480* a(3)^3 *dwMu0dtdw(3)+362880* a(2)* a(3)* a(4) *dwMu0dtdw(3)+181440* a(1)* a(4)^2 *dwMu0dtdw(3)+181440* a(2)^2* a(5) *dwMu0dtdw(3) + ...
%     362880* a(1)* a(3)* a(5) *dwMu0dtdw(3)+362880 *a(1)* a(2)* a(6) *dwMu0dtdw(3) +181440 *a(1)^2* a(7) *dwMu0dtdw(3)+ ...
%     60480* a(2)^3* a(3) *dwMu0dtdw(4)+181440* a(1)* a(2)* a(3)^2 *dwMu0dtdw(4)+181440* a(1)* a(2)^2 *a(4) *dwMu0dtdw(4)+ ...
%     181440* a(1)^2* a(3)* a(4) *dwMu0dtdw(4)+181440* a(1)^2* a(2)* a(5) *dwMu0dtdw(4)+60480* a(1)^3* a(6) *dwMu0dtdw(4)+ ...
%     15120* a(1)* a(2)^4 *dwMu0dtdw(5)+90720* a(1)^2* a(2)^2* a(3) *dwMu0dtdw(5)+30240* a(1)^3* a(3)^2 *dwMu0dtdw(5)+ ...
%     60480 *a(1)^3* a(2)* a(4) *dwMu0dtdw(5)+15120* a(1)^4* a(5) *dwMu0dtdw(5)+ ...
%     10080* a(1)^3* a(2)^3 *dwMu0dtdw(6)+15120* a(1)^4 *a(2)* a(3) *dwMu0dtdw(6)+3024* a(1)^5* a(4) *dwMu0dtdw(6)+ ...
%     1512* a(1)^5* a(2)^2 *dwMu0dtdw(7)+504* a(1)^6* a(3) *dwMu0dtdw(7)+ ...
%     72* a(1)^7 *a(2) *dwMu0dtdw(8)+a(1)^9 *dwMu0dtdw(9));
% 
% 
% 
% 
% b(10)=1/(720*7*8*9*10)*(3628800* a(10)* dwMu0dtdw(1)+ ...
%     1814400 * a(5)^2 *dwMu0dtdw(2)+3628800 *a(4)* a(6)* dwMu0dtdw(2)+3628800* a(3)* a(7)* dwMu0dtdw(2)+3628800 *a(2)* a(8)* dwMu0dtdw(2)+ ...
%     3628800* a(1)* a(9)* dwMu0dtdw(2)+ ...
%     1814400* a(3)^2* a(4)* dwMu0dtdw(3)+1814400* a(2)* a(4)^2* dwMu0dtdw(3)+3628800* a(2)* a(3)* a(5)* dwMu0dtdw(3)+3628800* a(1)* a(4)* a(5) *dwMu0dtdw(3)+ ...
%     1814400* a(2)^2* a(6)* dwMu0dtdw(3)+3628800* a(1)* a(3)* a(6)* dwMu0dtdw(3)+3628800 *a(1)* a(2)* a(7)* dwMu0dtdw(3)+1814400* a(1)^2* a(8)* dwMu0dtdw(3)+ ...
%     907200* a(2)^2 *a(3)^2* dwMu0dtdw(4)+604800* a(1)* a(3)^3* dwMu0dtdw(4)+604800* a(2)^3* a(4)* dwMu0dtdw(4)+3628800* a(1)* a(2)* a(3)* a(4)* dwMu0dtdw(4)+ ...
%     907200* a(1)^2* a(4)^2* dwMu0dtdw(4)+1814400* a(1)* a(2)^2* a(5)* dwMu0dtdw(4)+1814400* a(1)^2* a(3)* a(5)* dwMu0dtdw(4)+ ...
%     1814400* a(1)^2* a(2)* a(6)* dwMu0dtdw(4)+ 604800* a(1)^3* a(7) *dwMu0dtdw(4)+ ...
%     30240 *a(2)^5* dwMu0dtdw(5)+604800* a(1)* a(2)^3* a(3)* dwMu0dtdw(5)+907200 *a(1)^2* a(2)* a(3)^2* dwMu0dtdw(5)+907200* a(1)^2* a(2)^2* a(4) *dwMu0dtdw(5)+ ...
%     604800* a(1)^3* a(3)* a(4)* dwMu0dtdw(5)+604800* a(1)^3* a(2)* a(5)* dwMu0dtdw(5)+151200* a(1)^4* a(6)* dwMu0dtdw(5)+ ...
%     75600* a(1)^2* a(2)^4* dwMu0dtdw(6)+302400 *a(1)^3* a(2)^2* a(3)* dwMu0dtdw(6)+75600* a(1)^4 *a(3)^2* dwMu0dtdw(6)+151200* a(1)^4* a(2)* a(4) *dwMu0dtdw(6)+ ...
%     30240* a(1)^5* a(5)* dwMu0dtdw(6)+ ...
%     25200* a(1)^4* a(2)^3* dwMu0dtdw(7)+30240* a(1)^5* a(2)* a(3)* dwMu0dtdw(7)+5040* a(1)^6* a(4)* dwMu0dtdw(7)+ ...
%     2520* a(1)^6* a(2)^2* dwMu0dtdw(8)+720* a(1)^7* a(3)* dwMu0dtdw(8)+ ...
%     90 *a(1)^8* a(2)* dwMu0dtdw(9)+ ...
%     a(1)^10* dwMu0dtdw(10));



end

function [aH0,aH] = ConvertZCoeffsToHCoeffs(a0,a,SeriesOrder)

if(SeriesOrder==4)
aH0=a0+a(2)+3*a(4);
aH(1)=a(1)+3*a(3);
aH(2)=a(2)+6*a(4);
aH(3)=a(3);
aH(4)=a(4);
end

if(SeriesOrder==5)
aH0=a0+a(2)+3*a(4);
aH(1)=a(1)+3*a(3)+15*a(5);
aH(2)=a(2)+6*a(4);
aH(3)=a(3)+10*a(5);
aH(4)=a(4);
aH(5)=a(5);
end



if(SeriesOrder==6)
aH0=a0+a(2)+3*a(4)+15*a(6);
aH(1)=a(1)+3*a(3)+15*a(5);
aH(2)=a(2)+6*a(4)+45*a(6);
aH(3)=a(3)+10*a(5);
aH(4)=a(4)+15*a(6);
aH(5)=a(5);
aH(6)=a(6);
end


end

function [a0,a] = ConvertHCoeffsToZCoeffs(aH0,aH,SeriesOrder)

if(SeriesOrder==4)
a0=aH0-aH(2)+3*aH(4);
a(1)=aH(1)-3*aH(3);
a(2)=aH(2)-6*aH(4);
a(3)=aH(3);
a(4)=aH(4);
end

if(SeriesOrder==5)
a0=aH0-aH(2)+3*aH(4);
a(1)=aH(1)-3*aH(3)+15*aH(5);
a(2)=aH(2)-6*aH(4);
a(3)=aH(3)-10*aH(5);
a(4)=aH(4);
a(5)=aH(5);
end

if(SeriesOrder==6)
a0=aH0-aH(2)+3*aH(4)-15*aH(6);
a(1)=aH(1)-3*aH(3)+15*aH(5);
a(2)=aH(2)-6*aH(4)+45*aH(6);
a(3)=aH(3)-10*aH(5);
a(4)=aH(4)-15*aH(6);
a(5)=aH(5);
a(6)=aH(6);
end


end

function [] = FPESeriesCoeffsFromGaussiansWilmottNew02A_Lite_PIntegral02()

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

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

dt=.125/2/2/2;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.
Tt=32*2*2;%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;
Ts=Tt;
ds(1:Ts)=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*8);%+(64-16);
% ds(1:Ts)=dt/8;
% end
% if(Tt>4)&&(Tt<=20)
% Ts=(32)+((Tt-4)*2);
% ds(1:32)=dt/8;
% ds(33:Ts)=dt/2;
% end
% if(Tt>20)
% Ts=(32)+((20-4)*2)+(Tt-20);    
% ds(1:32)=dt/8;
% ds(33:64)=dt/2;
% ds(65:Ts)=dt;
% end
% T
% sum(ds(1:Ts))
% Ts
%str=input('Check if time is right');

OrderA=4;  %
OrderM=4;  %

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

dtM=dt;
TtM=Tt;



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

NnMid=4.0;

x0=1.0;   % starting value of SDE
beta1=.650;
beta2=0.0;   % Second drift term power.
gamma=.75;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=0.250;%.950;   %mean reversion parameter.
theta=1.00;%mean reversion target
sigma0=.450;%Volatility value


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

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

Z(1:Nn)=(((1:Nn)-7.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
%Ft(1:TtM+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0; %General time powers on hermite polynomials
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers on coefficients of hermite polynomials.
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers for bessel transformed coordinates.

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

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

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

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

ExpnOrder=5;
SeriesOrder=5;

wnStart=1;
tic

for tt=1:Ts
    tt

    if(tt==1)
        %dt=ds(1);
        %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder);
 %       [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),ExpnOrder);
        
        [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
    %    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
        
      %  dwMu0dtdw
      %  dwMu0dtdwb
      %  wMu0dt
      %  wMu0dtb
      %  str=input('Look at drift');
        c0=wMu0dt+w0;
        c(1)=sigma0*sqrt(ds(tt));
        c(2:10)=0.0;
        w0=c0;
        b0=wMu0dt;
        b(1:5)=0;
    
    elseif(tt<=3) 

    %dt=ds(tt);    
        
    %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
     %[b0,b] = CalculateDriftbCoeffs08O4(wMu0dt,dwMu0dtdw,c);
     [b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
     c0=c0+b0;
     c(1)=b(1)+sqrt(c(1)^2+sigma0^2*ds(tt));
     c(2)=b(2);
     c(3)=b(3);
     c(4)=b(4);
     c(5)=b(5);
     
     w0=c0;
            
    else
        
      
        
        
        
        
    %below wMu0dt is drift and dwMu0dtdw are its derivatives.
   %  [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
     [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
     [b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
     
     %size(b)
     %ba(1:6)=b(1:6);
     
     
    [wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
    
    [wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
    
    %[wVol3dt,dwVol3dtdw] = BesselVolAndDerivativesH3(w0,YqCoeff0,Fp1,gamma,ds(tt),8);

    [g10,g1] = CalculateDriftbCoeffs08A(wVol0dt,dwVol0dtdw,c,SeriesOrder);
    
    [g20,g2] = CalculateDriftbCoeffs08A(wVol2dt,dwVol2dtdw,c,SeriesOrder);
 

 a0=c0+b0;
 
 a(1:5)=c(1:5)+b(1:5);
    
     [a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A_Lite(a0,a,g10,g1,g20,g2,sigma0,ds(tt));
     %[a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A(a0,a,b0,ba,g0,g,h0,h,sigma0,ds(tt));
     
     
 c0=a0;
 c(1)=a1;
 c(2)=a2;
 c(3)=a3;
 c(4)=a4;
 c(5)=a5;
    
    w0=c0;
    
    end
   
   
    
  %Below Calculate the Z-series expansion for SDE variable in original coordinates 
  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;

  
    
 if(tt>1)
    fH0prev=fH0;%+bH0;
    fHprev=fH;%+bH;
 end

 %Below f0 and f define the Z_series in original coordinates
 [f0,f] = CalculateDriftbCoeffs08A(yw0,dyw0,c,SeriesOrder);
 %fH0 and fH(1:5) define the coefficients of hermite polynomial 
 %representation of SDE variable in original coordinates.  
 [fH0,fH] = ConvertZCoeffsToHCoeffs(f0,f,SeriesOrder);    
 
 [yMu0dt,dyMu0dtdy] = OriginalDriftAndDerivativesH0(yw0,YCoeff0,Fp,ds(tt),ExpnOrder)
  [by0,by] = CalculateDriftbCoeffs08A(yMu0dt,dyMu0dtdy,f,SeriesOrder);
 [byH0,byH] = ConvertZCoeffsToHCoeffs(by0,by,SeriesOrder);   

 by
 by0
 byH0
 byH
 
 %str=input('Look at original drift');
 
 if(tt>1)
    fH0prev=fH0prev+byH0;
    fHprev=fHprev+byH;
 end
 
 
%Below calculate the Hermite representation of orthogonal increments in the
%original coordinates variable at each time step by squared subtraction of
%coefficients of hermite polynomials.
%Here cH0 and cH define the hermite polynomial coordinates of orthogonal
 %increments at each simulation step. The take double array index in 
 %simulation level, tt and hermite polynomial count (1:5).  
if(tt==1) 
 cH0=fH0;
 cH(tt,1:SeriesOrder)=fH(1:SeriesOrder);

 bt0=by0;
 bt(tt,1:SeriesOrder)=by(1:SeriesOrder);
end
if(tt>1)       
 dH0=fH0-fH0prev;
    
 dH(1:SeriesOrder)=sign(fH(1:SeriesOrder)-fHprev(1:SeriesOrder)).*sqrt(abs(sign(fH(1:SeriesOrder)).*fH(1:SeriesOrder).^2-sign(fHprev(1:SeriesOrder)).*fHprev(1:SeriesOrder).^2));    
 
 cH0(tt)=dH0;
 cH(tt,1:SeriesOrder)=dH(1:SeriesOrder);
 
 bt0(tt)=by0;
 bt(tt,1:SeriesOrder)=by(1:SeriesOrder);
end



end
wnStart=1;


yy0=((1-gamma).*w0).^(1/(1-gamma));

 w(1:Nn)=c0;
 for nn=1:ExpnOrder
     w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^nn;
 end
 
  %Below calculate the hermite representation of dt-integral in original
 %coordinates
 %Here cH0 and cH define the hermite polynomial coordinates of orthogonal
 %increments at each simulation step and bH0 and bH define the hermite
 %representation of the dt-integral after continued squared addition of 
 %scaled/weighted values of hermite coefficients of orthogonal increments.
 
 
 eH0=cH0(1)*dt*Ts;
 eH(1:SeriesOrder)=cH(1,1:SeriesOrder)*dt*Ts;
 
 for tt=2:Ts
    ts=Ts-tt+1;
    eH0=eH0+cH0(tt)*dt*ts;
    eH(1:SeriesOrder)=sign(eH(1:SeriesOrder)+dt*(ts).*cH(tt,1:SeriesOrder)).*sqrt(abs(sign(eH(1:SeriesOrder)).*eH(1:SeriesOrder).^2+sign(dt*(ts).*cH(tt,1:SeriesOrder)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrder).^2));   
 end
 
 %We convert hermite represenation given by bH0 and bH(1:5) into 
 %dt-integral Z-series Coeffs, e0 and e(1:5)
 
 [e0,e] = ConvertHCoeffsToZCoeffs(eH0,eH,SeriesOrder);%  

 
   e0=e0+bt0(1)*dt*Ts;
  e(1:SeriesOrder)=e(1:SeriesOrder)+bt(1:SeriesOrder)*dt*Ts;
  for tt=2:Ts
     ts=Ts-tt+1;
     e0=e0+bt0(tt)*dt*ts;
     e(1:SeriesOrder)=e(1:SeriesOrder)+bt(tt,1:SeriesOrder)*dt*ts;
     %eH(1:SeriesOrder)=sign(eH(1:SeriesOrder)+dt*(ts).*cH(tt,1:SeriesOrder)).*sqrt(abs(sign(eH(1:SeriesOrder)).*eH(1:SeriesOrder).^2+sign(dt*(ts).*cH(tt,1:SeriesOrder)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrder).^2));   
  end
  
 
 
  yydt(1:Nn)=e0;
 for nn=1:ExpnOrder
     yydt(1:Nn)=yydt(1:Nn)+e(nn)*Z(1:Nn).^nn;
 end

 
 %c0
 %c
 
 
 %str=input('Look at numbers')
 yy0(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));

 yy=yy0;
Dfyy(wnStart:Nn)=0;
Dfyydt(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    
    Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    Dfyydt(nn) = (yydt(nn + 1) - yydt(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;
pyydt(1:Nn)=0;
for nn = wnStart+1:Nn-1
    
    pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
    pyydt(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyydt(nn));
end

toc

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


theta1=1;
rng(29079137, 'twister')
paths=50000;
YY(1:paths)=x0;  %Original process monte carlo.
YYdt(1:paths)=0;   

Random1(1:paths)=0;
for tt=1:TtM
    
    
    Random1=randn(size(Random1));
    HermiteP1(1,1:paths)=1;
    HermiteP1(2,1:paths)=Random1(1:paths);
    HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
    HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
    HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;


    YYdt(1:paths)=YYdt(1:paths)+.5*YY(1:paths)*dtM;   

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

    
end




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


MaxCutOff=30;
NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too 
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
plot(yy(wnStart+1:Nn-1),pyy(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 %plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');
 
%title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
title(sprintf('x0 = %.4f,mu1=%.3f,beta1=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,mu1,beta1,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,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.');

NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa))/2;
[YdtDensity,IndexOutYdt,IndexMaxYdt] = MakeDensityFromSimulation_Infiniti_NEW(YYdt,paths,NoOfBins,MaxCutOff );
plot(yydt(wnStart+1:Nn-1),pyydt(wnStart+1:Nn-1),'r',IndexOutYdt(1:IndexMaxYdt),YdtDensity(1:IndexMaxYdt),'g');

%title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
title(sprintf('x0 = %.4f,mu1=%.3f,beta1=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,mu1,beta1,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));


legend({'Ito-Hermite Path Integral Density','Monte Carlo Path Integral Density'},'Location','northeast')
 
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');



end

function [] = FPESeriesCoeffs_Lite_PIntegral08AB_Drifting()

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

%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=32*2*8;%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;
Ts=Tt;
ds(1:Ts)=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*8);%+(64-16);
% ds(1:Ts)=dt/8;
% end
% if(Tt>4)&&(Tt<=20)
% Ts=(32)+((Tt-4)*2);
% ds(1:32)=dt/8;
% ds(33:Ts)=dt/2;
% end
% if(Tt>20)
% Ts=(32)+((20-4)*2)+(Tt-20);    
% ds(1:32)=dt/8;
% ds(33:64)=dt/2;
% ds(65:Ts)=dt;
% end
% T
% sum(ds(1:Ts))
% Ts
%str=input('Check if time is right');

OrderA=4;  %
OrderM=4;  %

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

dtM=dt;
TtM=Tt;



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

NnMid=4.0;

x0=1;   % starting value of SDE

%beta1=0.0;
beta2=0.0;   % Second drift term power.
beta1=0.9950;
%beta2=0.0;   % Second drift term power.
gamma=.6;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=0.0;%.950;   %mean reversion parameter.
theta=1.0;%mean reversion target
sigma0=.250;%Volatility value


%you can specify any general mu1 and mu2 and beta1 and beta2.
mu1=-.150;   %first drift coefficient.
mu2=0;%-1*kappa;    % Second drift coefficient.
%mu1=theta*kappa;
%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=x0^(1-gamma)/(1-gamma);

Z(1:Nn)=(((1:Nn)-7.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
%Ft(1:TtM+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0; %General time powers on hermite polynomials
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers on coefficients of hermite polynomials.
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers for bessel transformed coordinates.

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

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

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

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

ExpnOrder=5;
SeriesOrder=5;

wnStart=1;
tic

for tt=1:Ts
    tt

    if(tt==1)
        %dt=ds(1);
        %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder);
 %       [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),ExpnOrder);
        
        [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
    %    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
        
      %  dwMu0dtdw
      %  dwMu0dtdwb
      %  wMu0dt
      %  wMu0dtb
      %  str=input('Look at drift');
        c0=wMu0dt+w0;
        c(1)=sigma0*sqrt(ds(tt));
        c(2:10)=0.0;
        w0=c0;
        b0=wMu0dt;
        b(1:5)=0;
    
    elseif(tt<=3) 

    %dt=ds(tt);    
        
    %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
     %[b0,b] = CalculateDriftbCoeffs08O4(wMu0dt,dwMu0dtdw,c);
     [b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
     c0=c0+b0;
     c(1)=b(1)+sqrt(c(1)^2+sigma0^2*ds(tt));
     c(2)=b(2);
     c(3)=b(3);
     c(4)=b(4);
     c(5)=b(5);
     
     w0=c0;
            
    else
        
      
        
        
        
        
    %below wMu0dt is drift and dwMu0dtdw are its derivatives.
   %  [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
     [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
     [b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
     
     %size(b)
     %ba(1:6)=b(1:6);
     
     
    [wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
    
    [wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
    
    %[wVol3dt,dwVol3dtdw] = BesselVolAndDerivativesH3(w0,YqCoeff0,Fp1,gamma,ds(tt),8);

    [g10,g1] = CalculateDriftbCoeffs08A(wVol0dt,dwVol0dtdw,c,SeriesOrder);
    
    [g20,g2] = CalculateDriftbCoeffs08A(wVol2dt,dwVol2dtdw,c,SeriesOrder);
 

 a0=c0+b0;
 
 a(1:5)=c(1:5)+b(1:5);
    
     [a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A_Lite(a0,a,g10,g1,g20,g2,sigma0,ds(tt));
     %[a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A(a0,a,b0,ba,g0,g,h0,h,sigma0,ds(tt));
     
     
 c0=a0;
 c(1)=a1;
 c(2)=a2;
 c(3)=a3;
 c(4)=a4;
 c(5)=a5;
    
    w0=c0;
    
    end
   
   
    
  %Below Calculate the Z-series expansion for SDE variable in original coordinates 
  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;
 
  
  if(tt==1)
    bSumH0=0;%+bH0;
    bSumH(1:SeriesOrder+2)=0;%+bH;
    f0prev=x0;
    fprev(1:SeriesOrder)=0;
 end
    
 if(tt>1)
    fH0prev=fH0;%+bH0;
    fHprev=fH;%+bH;
    f0prev=f0;
    fprev=f;
 end

 %Below f0 and f define the Z_series in original coordinates
 [f0,f] = CalculateDriftbCoeffs08A(yw0,dyw0,c,SeriesOrder);
 f(6:8)=0; 
 
 SeriesOrdery=SeriesOrder;
 
 
 %fH0 and fH(1:5) define the coefficients of hermite polynomial 
 %representation of SDE variable in original coordinates.  
 [fH0,fH] = ConvertZCoeffsToHCoeffs(f0,f,SeriesOrdery);    
 
 [yMu0dt,dyMu0dtdy] = OriginalDriftAndDerivativesH0(f0prev,YCoeff0,Fp,ds(tt),7)
  [by0,by] = CalculateDriftbCoeffs08A(yMu0dt,dyMu0dtdy,fprev,SeriesOrdery);
 [byH0,byH] = ConvertZCoeffsToHCoeffs(by0,by,SeriesOrdery);   

 %by
 %by0
 %byH0
 %byH
 
 fH(SeriesOrdery+1:7)=0;
 byH(SeriesOrdery+1:7)=0;
 
 
 %str=input('Look at original drift');
 
% if(tt>1)
%    fH0prev=fH0prev+byH0;
%    fHprev=fHprev+byH;
% end
 
 if(tt>=1)
    
     bSumH0prev=bSumH0;
     bSumHprev=bSumH; 
     %fH0prev=fH0prev+byH0;
     %fHprev=fHprev;%-exp(-kappa*dt).*bSumH;
     bSumH0=bSumH0+byH0;
     bSumH=bSumH+byH; 
    
    
 end
 
%Below calculate the Hermite representation of orthogonal increments in the
%original coordinates variable at each time step by squared subtraction of
%coefficients of hermite polynomials.
%Here cH0 and cH define the hermite polynomial coordinates of orthogonal
 %increments at each simulation step. The take double array index in 
 %simulation level, tt and hermite polynomial count (1:5).  
if(tt==1) 
 cH0(tt)=fH0;
 cH(tt,1:SeriesOrdery)=fH(1:SeriesOrdery)-byH(1:SeriesOrdery);

 bt0=by0;
 bt(tt,1:SeriesOrdery)=by(1:SeriesOrdery);
end
if(tt>1)       
 %dH0=fH0-fH0prev;
 dH0=fH0-bSumH0-fH0prev+bSumH0prev;
    
 %dH(1:SeriesOrdery)=sign(fH(1:SeriesOrdery)-fHprev(1:SeriesOrdery)).*sqrt(abs(sign(fH(1:SeriesOrdery)).*fH(1:SeriesOrdery).^2-sign(fHprev(1:SeriesOrdery)).*fHprev(1:SeriesOrdery).^2));    
 dH(1:SeriesOrdery)=sign(fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)-(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery))).*sqrt(abs(sign(fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)).*(fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)).^2-sign(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery)).*(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery)).^2));    
 
 
 cH0(tt)=dH0;
 cH(tt,1:SeriesOrdery)=dH(1:SeriesOrdery);
 
 bt0(tt)=by0;
 bt(tt,1:SeriesOrdery)=by(1:SeriesOrdery);
end



end
wnStart=1;


yy0=((1-gamma).*w0).^(1/(1-gamma));

 w(1:Nn)=c0;
 for nn=1:ExpnOrder
     w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^nn;
 end
 
  %Below calculate the hermite representation of dt-integral in original
 %coordinates
 %Here cH0 and cH define the hermite polynomial coordinates of orthogonal
 %increments at each simulation step and bH0 and bH define the hermite
 %representation of the dt-integral after continued squared addition of 
 %scaled/weighted values of hermite coefficients of orthogonal increments.
 
 
 eH0=cH0(1)*dt*Ts;
 %eH(1:SeriesOrdery)=cH(1,1:SeriesOrdery)*dt*Ts;
 
 for tt=2:Ts
 %   ts=Ts-tt+1;
 %   eH0=eH0+cH0(tt)*dt*ts;
 %   eH(1:SeriesOrdery)=sign(eH(1:SeriesOrdery)+dt*(ts).*cH(tt,1:SeriesOrdery)).*sqrt(abs(sign(eH(1:SeriesOrdery)).*eH(1:SeriesOrdery).^2+sign(dt*(ts).*cH(tt,1:SeriesOrdery)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrdery).^2));   
 end
 
 
 %eH2(1:SeriesOrdery)=sign(cH(1,1:SeriesOrdery)*dt*Ts).*(cH(1,1:SeriesOrdery)*dt*Ts).^2;
 
 
 eH2(1:SeriesOrdery)=0;
 eH0=0;
 
 
 
 for tt=2:Ts
 %   ts=Ts-tt+1;
 %   eH0=eH0+cH0(tt)*dt*ts;
 %   eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+sign(dt*(ts).*cH(tt,1:SeriesOrdery)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrdery).^2;   
    %  for ss=1:tt-1
    %      tss=tt-ss;
    %      %tss=1;
    %      %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(tt+ss)).*sign(dt*(tss).*cH(ss,1:SeriesOrdery)).*(dt).^2.*tss.^2.*cH(tt,1:SeriesOrdery).^2;   
    %      eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(tt*.5)).*sign(dt.*cH(ss,1:SeriesOrdery)).*dt^2.*2.*tss.*cH(ss,1:SeriesOrdery).^2;   
    %  end
 end
 Coeff0(1:Ts)=0;
 Coeff(1:Ts)=0;
 Coeff1(1:Ts)=0;
 Coeff2(1:Ts)=0;
%  for tt=1:Ts
%      Coeff2(tt)=Coeff2(tt)+0.*(Ts-tt+1).*dt^2;
%  end
 for tt=1:Ts
     for ss=tt:Ts
         for rr=ss:Ts
            if(rr==ss)
            
            Coeff2(tt)=Coeff2(tt)+dt^2;    
            else
            
            Coeff2(tt)=Coeff2(tt)+2.*dt^2;
            end

         end
     end
 end
 
 for ss=1:Ts
     for tt=ss:Ts
     
        Coeff1(ss)=Coeff1(ss)+dt;
     end
 end
 
 
 for tt=1:Ts
    ts=Ts-tt+1;
    eH0=eH0+cH0(tt)*Coeff1(tt);
    eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+sign(cH(tt,1:SeriesOrdery)).*Coeff2(tt).*cH(tt,1:SeriesOrdery).^2;   
 end
 
 
 

 eH(1:SeriesOrdery)=sign(eH2(1:SeriesOrdery)).*sqrt(abs(eH2(1:SeriesOrdery)));
 
 
 %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;
 [e0,e] = ConvertHCoeffsToZCoeffs(eH0,eH,SeriesOrdery);%  


for ss=1:Ts
     %for ss=1:tt
     for tt=ss:Ts    
     
        e(1:SeriesOrdery)=e(1:SeriesOrdery)+bt(ss,1:SeriesOrdery)*dt;   
        e0=e0+bt0(ss)*dt; 
        
     end
end

  yydt(1:Nn)=e0;
 for nn=1:SeriesOrdery
     yydt(1:Nn)=yydt(1:Nn)+e(nn)*Z(1:Nn).^nn;
 end

 
 %c0
 %c
 
 
 %str=input('Look at numbers')
 yy0(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));

 yy=yy0;
Dfyy(wnStart:Nn)=0;
Dfyydt(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    
    Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    Dfyydt(nn) = (yydt(nn + 1) - yydt(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;
pyydt(1:Nn)=0;
for nn = wnStart+1:Nn-1
    
    pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
    pyydt(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyydt(nn));
end

toc

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


theta1=1;
rng(29079137, 'twister')
paths=50000;
YY(1:paths)=x0;  %Original process monte carlo.
YYdt(1:paths)=0;   

Random1(1:paths)=0;
for tt=1:TtM
    
    
    Random1=randn(size(Random1));
    HermiteP1(1,1:paths)=1;
    HermiteP1(2,1:paths)=Random1(1:paths);
    HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
    HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
    HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;


    %YYdt(1:paths)=YYdt(1:paths)+.5*YY(1:paths)*dtM;   

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

    
end




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


MaxCutOff=100;
NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too 
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
plot(yy(wnStart+1:Nn-1),pyy(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 %plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');
 
%title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
title(sprintf('x0 = %.4f,mu1=%.3f,beta1=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,mu1,beta1,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,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.');

NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa))/2;
[YdtDensity,IndexOutYdt,IndexMaxYdt] = MakeDensityFromSimulation_Infiniti_NEW(YYdt,paths,NoOfBins,MaxCutOff );
plot(yydt(wnStart+1:Nn-1),pyydt(wnStart+1:Nn-1),'r',IndexOutYdt(1:IndexMaxYdt),YdtDensity(1:IndexMaxYdt),'g');

%title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
title(sprintf('x0 = %.4f,mu1=%.3f,beta1=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,mu1,beta1,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));


legend({'Ito-Hermite Path Integral Density','Monte Carlo Path Integral Density'},'Location','northeast')
 
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');



end

function [] = FPESeriesCoeffsFromGaussiansWilmottNew02A_Lite_PIntegral08CC()

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

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

dt=.125/2/2/2;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.
Tt=32*2*3;%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;
Ts=Tt;
ds(1:Ts)=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*8);%+(64-16);
% ds(1:Ts)=dt/8;
% end
% if(Tt>4)&&(Tt<=20)
% Ts=(32)+((Tt-4)*2);
% ds(1:32)=dt/8;
% ds(33:Ts)=dt/2;
% end
% if(Tt>20)
% Ts=(32)+((20-4)*2)+(Tt-20);    
% ds(1:32)=dt/8;
% ds(33:64)=dt/2;
% ds(65:Ts)=dt;
% end
% T
% sum(ds(1:Ts))
% Ts
%str=input('Check if time is right');

OrderA=4;  %
OrderM=4;  %

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

dtM=dt;
TtM=Tt;



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

NnMid=4.0;

x0=.10;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
%beta1=0.950;
%beta2=0.0;   % Second drift term power.
gamma=.85;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=2.00;%.950;   %mean reversion parameter.
theta=.10;%mean reversion target
sigma0=.750;%Volatility value


%you can specify any general mu1 and mu2 and beta1 and beta2.
mu1=.50;   %first drift coefficient.
mu2=0;%-1*kappa;    % Second drift coefficient.
mu1=theta*kappa;
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=x0^(1-gamma)/(1-gamma);

Z(1:Nn)=(((1:Nn)-7.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
%Ft(1:TtM+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0; %General time powers on hermite polynomials
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers on coefficients of hermite polynomials.
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers for bessel transformed coordinates.

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

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

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

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

ExpnOrder=5;
SeriesOrder=5;

wnStart=1;
tic

for tt=1:Ts
    tt

    if(tt==1)
        %dt=ds(1);
        %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder);
 %       [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),ExpnOrder);
        
        [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
    %    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
        
      %  dwMu0dtdw
      %  dwMu0dtdwb
      %  wMu0dt
      %  wMu0dtb
      %  str=input('Look at drift');
        c0=wMu0dt+w0;
        c(1)=sigma0*sqrt(ds(tt));
        c(2:10)=0.0;
        w0=c0;
        b0=wMu0dt;
        b(1:5)=0;
    
    elseif(tt<=3) 

    %dt=ds(tt);    
        
    %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
     %[b0,b] = CalculateDriftbCoeffs08O4(wMu0dt,dwMu0dtdw,c);
     [b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
     c0=c0+b0;
     c(1)=b(1)+sqrt(c(1)^2+sigma0^2*ds(tt));
     c(2)=b(2);
     c(3)=b(3);
     c(4)=b(4);
     c(5)=b(5);
     
     w0=c0;
            
    else
        
      
        
        
        
        
    %below wMu0dt is drift and dwMu0dtdw are its derivatives.
   %  [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
     [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
     [b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
     
     %size(b)
     %ba(1:6)=b(1:6);
     
     
    [wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
    
    [wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
    
    %[wVol3dt,dwVol3dtdw] = BesselVolAndDerivativesH3(w0,YqCoeff0,Fp1,gamma,ds(tt),8);

    [g10,g1] = CalculateDriftbCoeffs08A(wVol0dt,dwVol0dtdw,c,SeriesOrder);
    
    [g20,g2] = CalculateDriftbCoeffs08A(wVol2dt,dwVol2dtdw,c,SeriesOrder);
 

 a0=c0+b0;
 
 a(1:5)=c(1:5)+b(1:5);
    
     [a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A_Lite(a0,a,g10,g1,g20,g2,sigma0,ds(tt));
     %[a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A(a0,a,b0,ba,g0,g,h0,h,sigma0,ds(tt));
     
     
 c0=a0;
 c(1)=a1;
 c(2)=a2;
 c(3)=a3;
 c(4)=a4;
 c(5)=a5;
    
    w0=c0;
    
    end
   
   
    
  %Below Calculate the Z-series expansion for SDE variable in original coordinates 
  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;
 
  
  if(tt==1)
    bSumH0=0;%+bH0;
    bSumH(1:SeriesOrder+2)=0;%+bH;
    fH0prev=x0;%+bH0;
    fHprev(1:SeriesOrder)=0;
    f0prev=x0;
    fprev(1:SeriesOrder)=0;
 end
    
 if(tt>1)
    fH0prev=fH0;%+bH0;
    fHprev=fH;%+bH;
    f0prev=f0;%+bH0;
    fprev=f;%+bH;
 end

 %Below f0 and f define the Z_series in original coordinates
 [f0,f] = CalculateDriftbCoeffs08A(yw0,dyw0,c,SeriesOrder);
 f(6:8)=0; 
 
 SeriesOrdery=SeriesOrder;
 
 
 %fH0 and fH(1:5) define the coefficients of hermite polynomial 
 %representation of SDE variable in original coordinates.  
 [fH0,fH] = ConvertZCoeffsToHCoeffs(f0,f,SeriesOrdery);    
 
% [yMu0dt,dyMu0dtdy] = OriginalDriftAndDerivativesH0(f0prev,YCoeff0,Fp,ds(tt),7)
%  [by0,by] = CalculateDriftbCoeffs08A(yMu0dt,dyMu0dtdy,fprev,SeriesOrdery);
% [byH0,byH] = ConvertZCoeffsToHCoeffs(by0,by,SeriesOrdery);   

% fH0
% byH0
% fH
% 
% byH
% fH./byH
% 
% 1/(exp(kappa*dt)-1)
% 1/(1-exp(-kappa*dt))
% byH./fH
% (1-exp(-kappa*dt))
%  (fHprev(1:SeriesOrder)+byH(1:SeriesOrder))./fHprev(1:SeriesOrder)
%  exp(-kappa*dt)
%  %fH+byH
%  %(fH+byH)./fH
%  
% by0
% byH0
% byH
 
 fH(SeriesOrdery+1:7)=0;
% byH(SeriesOrdery+1:7)=0;
 
 
% str=input('Look at original drift');
 
% if(tt>1)
%    fH0prev=fH0prev+byH0;
%    fHprev=fHprev+byH;
% end
 
 if(tt>=1)
    
     bSumH0prev=0;
     bSumHprev(1:SeriesOrdery)=0; 
     %fH0prev=fH0prev+byH0;
     %fHprev=fHprev;%-exp(-kappa*dt).*bSumH;
     bSumH0=0;
     bSumH(1:SeriesOrdery)=0; 
    
    
 end
 
%Below calculate the Hermite representation of orthogonal increments in the
%original coordinates variable at each time step by squared subtraction of
%coefficients of hermite polynomials.
%Here cH0 and cH define the hermite polynomial coordinates of orthogonal
 %increments at each simulation step. The take double array index in 
 %simulation level, tt and hermite polynomial count (1:5).  
if(tt==1) 
 cH0(tt)=fH0;
 cH(tt,1:SeriesOrdery)=fH(1:SeriesOrdery);%-byH(1:SeriesOrdery);

% bt0=by0;
% bt(tt,1:SeriesOrdery)=by(1:SeriesOrdery);
end
if(tt>1)       
 %dH0=fH0-fH0prev;
 dH0=fH0-bSumH0-fH0prev+bSumH0prev;
    
 %dH(1:SeriesOrdery)=sign(fH(1:SeriesOrdery)-fHprev(1:SeriesOrdery)).*sqrt(abs(sign(fH(1:SeriesOrdery)).*fH(1:SeriesOrdery).^2-sign(fHprev(1:SeriesOrdery)).*fHprev(1:SeriesOrdery).^2));    
 dH(1:SeriesOrdery)=sign(fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)-exp(-1*kappa*dt).*(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery))).*sqrt(abs(sign(fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)).*(exp(0*kappa*dt).*fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)).^2-sign(exp(-kappa*dt).*fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery)).*exp(-2*1*kappa*dt).*(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery)).^2));    
 %dH(1:SeriesOrdery)=sign(fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)-exp(-kappa*dt).*(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery))).*sqrt(abs(sign(fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)).*(exp(1*kappa*tt*dt).*fH(1:SeriesOrdery)-bSumH(1:SeriesOrdery)).^2-sign(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery)).*exp(2*(tt-1)*kappa*dt).*(fHprev(1:SeriesOrdery)-bSumHprev(1:SeriesOrdery)).^2));    
 
 
 cH0(tt)=dH0;
 cH(tt,1:SeriesOrdery)=dH(1:SeriesOrdery);
 
% bt0(tt)=by0;
% bt(tt,1:SeriesOrdery)=by(1:SeriesOrdery);
end



end
wnStart=1;


yy0=((1-gamma).*w0).^(1/(1-gamma));

 w(1:Nn)=c0;
 for nn=1:ExpnOrder
     w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^nn;
 end
 
  %Below calculate the hermite representation of dt-integral in original
 %coordinates
 %Here cH0 and cH define the hermite polynomial coordinates of orthogonal
 %increments at each simulation step and bH0 and bH define the hermite
 %representation of the dt-integral after continued squared addition of 
 %scaled/weighted values of hermite coefficients of orthogonal increments.
 
 
 eH0=cH0(1)*dt*Ts;
 %eH(1:SeriesOrdery)=cH(1,1:SeriesOrdery)*dt*Ts;
 
 for tt=2:Ts
 %   ts=Ts-tt+1;
 %   eH0=eH0+cH0(tt)*dt*ts;
 %   eH(1:SeriesOrdery)=sign(eH(1:SeriesOrdery)+dt*(ts).*cH(tt,1:SeriesOrdery)).*sqrt(abs(sign(eH(1:SeriesOrdery)).*eH(1:SeriesOrdery).^2+sign(dt*(ts).*cH(tt,1:SeriesOrdery)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrdery).^2));   
 end
 
 
 %eH2(1:SeriesOrdery)=sign(cH(1,1:SeriesOrdery)*dt*Ts).*(cH(1,1:SeriesOrdery)*dt*Ts).^2;
 
 
 eH2(1:SeriesOrdery)=0;
 eH0=0;
 
 
 
 for tt=2:Ts
 %   ts=Ts-tt+1;
 %   eH0=eH0+cH0(tt)*dt*ts;
 %   eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+sign(dt*(ts).*cH(tt,1:SeriesOrdery)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrdery).^2;   
    %  for ss=1:tt-1
    %      tss=tt-ss;
    %      %tss=1;
    %      %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(tt+ss)).*sign(dt*(tss).*cH(ss,1:SeriesOrdery)).*(dt).^2.*tss.^2.*cH(tt,1:SeriesOrdery).^2;   
    %      eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(tt*.5)).*sign(dt.*cH(ss,1:SeriesOrdery)).*dt^2.*2.*tss.*cH(ss,1:SeriesOrdery).^2;   
    %  end
 end
 Coeff0(1:Ts)=0;
 Coeff(1:Ts)=0;
 Coeff1(1:Ts)=0;
 Coeff2(1:Ts)=0;
 for tt=1:Ts
     Coeff2(tt)=Coeff2(tt)+0.*(Ts-tt+1).*dt^2;
 end
 for tt=1:Ts
     for ss=tt:Ts
         for rr=ss:Ts
            %Coeff2(tt)=Coeff2(tt)+0*.75* exp(-.5*kappa*(dt*((rr)-ss))).*dt^2;
            %Coeff2(tt)=Coeff2(tt)+2/.77* exp(-1*kappa*(dt*((rr)-ss))).*dt^2;
            %coef=max(.8-.8/sqrt(1)/kappa,0);
            %coef=max(.75-.8/sqrt(2)/kappa,0);
            %coef=max(.75-.8/sqrt(4)/sqrt(kappa)+2/3*sqrt(kappa)*sqrt(Ts*dt)-sqrt(kappa)*sqrt(.5),0);
            %Coeff2(tt)=Coeff2(tt)+(2+coef)* exp(-1*kappa*(dt*((rr)-ss))).*dt^2;
            %Coeff2(tt)=Coeff2(tt)+(2+coef)* exp(-1*kappa*(dt*((rr)-ss))).*dt^2;
            if(rr==ss)
            
            Coeff2(tt)=Coeff2(tt)+(exp(-0*kappa*dt*(tt))).*(exp(-2*kappa*dt*(rr-ss))).*dt^2;    
            else
            %Coeff2(tt)=Coeff2(tt)+sqrt(3)/sqrt(kappa)*(exp(-2*kappa*dt*(rr-ss))).*dt^2;
            Coeff2(tt)=Coeff2(tt)+2*(exp(-2*kappa*dt*(rr-ss)))*(exp(-1*kappa*dt*(tt))).*dt^2;
            end
          %  rr
          %  ss
          %  (exp(-1*kappa*dt*(rr-ss)))
          %  (exp(-2*kappa*dt*(rr-ss)))
          %  str=input('Look at numbers');
            %Coeff2(tt)=Coeff2(tt)+(2)* (exp(-kappa*dt*(rr-ss)))*(1+.125*kappa*dt*(ss)).*dt^2;
            %Coeff2(tt)=Coeff2(tt)+(2)* (exp(-kappa*dt*(rr-ss)))*(1/.8-.125/4*kappa*dt*(tt)/2).*dt^2;
            %Coeff2(tt)=Coeff2(tt)+1.*(exp(-kappa*dt*(tt/2))).*dt^2;
            
            %Coeff2(tt)=Coeff2(tt)+1/(2).* dt^2;
         end
     end
 end
 
 
 
 for tt=1:Ts
    ts=Ts-tt+1;
    eH0=eH0+cH0(tt)*ts*dt;
    eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+sign(cH(tt,1:SeriesOrdery)).*Coeff2(tt).*cH(tt,1:SeriesOrdery).^2;   
 end
 
 
 
%   for ss=1:Ts-1
%         %tss=tt-ss;
%         tss=Ts-ss+1;
%         %tss=1;
%         %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(tt+ss)).*sign(dt*(tss).*cH(ss,1:SeriesOrdery)).*(dt).^2.*tss.^2.*cH(tt,1:SeriesOrdery).^2;   
%         eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(Ts+ss)).*sign(dt.*cH(ss,1:SeriesOrdery)).*dt^2.*tss^2.*cH(ss,1:SeriesOrdery).^2;   
%  end

 eH(1:SeriesOrdery)=sign(eH2(1:SeriesOrdery)).*sqrt(abs(eH2(1:SeriesOrdery)));
 
 
 %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;
 [e0,e] = ConvertHCoeffsToZCoeffs(eH0,eH,SeriesOrdery);%  


for ss=1:Ts
     %for ss=1:tt
     for tt=ss:Ts    
     
        %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(tt+ss)).*sign(dt*(tss).*cH(ss,1:SeriesOrdery)).*(dt).^2.*tss.^2.*cH(tt,1:SeriesOrdery).^2;   
        %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(Ts*.5)).*sign(dt.*cH(ss,1:SeriesOrdery)).*dt^2.*tss^2.*cH(ss,1:SeriesOrdery).^2;   
        %e(1:SeriesOrdery)=e(1:SeriesOrdery)+1/sqrt(2^0)*exp(-kappa*dt*(tt-ss)).*bt(ss,1:SeriesOrdery)*dt;   
%        e(1:SeriesOrdery)=e(1:SeriesOrdery)+exp(-1*kappa*dt*(tt-ss)).*bt(ss,1:SeriesOrdery)*dt;   
%        e0=e0+exp(-kappa*dt*(tt-ss)).*bt0(ss)*dt; 
%        e(1:SeriesOrdery)=0;   
%        e0=0; 

        
     end
end

% for tt=1:Ts
%      for ss=tt:Ts
%      
%         %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(tt+ss)).*sign(dt*(tss).*cH(ss,1:SeriesOrdery)).*(dt).^2.*tss.^2.*cH(tt,1:SeriesOrdery).^2;   
%         %eH2(1:SeriesOrdery)=eH2(1:SeriesOrdery)+exp(-kappa*dt*(Ts*.5)).*sign(dt.*cH(ss,1:SeriesOrdery)).*dt^2.*tss^2.*cH(ss,1:SeriesOrdery).^2;   
%         %e(1:SeriesOrdery)=e(1:SeriesOrdery)+exp(-kappa*dt*(tt-ss)).*bt(ss,1:SeriesOrdery)*dt;   
%         e(1:SeriesOrdery)=e(1:SeriesOrdery)+exp(-kappa*dt*(ss-tt)).*bt(tt,1:SeriesOrdery)*dt;   
%         %e0=e0+exp(-kappa*dt*(tt-ss)).*bt0(ss)*dt; 
%         e0=e0+exp(-kappa*dt*(ss-tt)).*bt0(tt)*dt; 
%      end
% end

%   e0=e0+bt0(1)*dt*Ts;
%   e(1:SeriesOrdery)=e(1:SeriesOrdery)+bt(1:SeriesOrdery)*dt*Ts;
%   for tt=2:Ts-1
%      ts=Ts-tt;
%      e0=e0+bt0(tt)*dt*ts;
%      e(1:SeriesOrdery)=e(1:SeriesOrdery)+bt(tt,1:SeriesOrdery)*dt*ts.*exp(-kappa*ts*dt);
%      %eH(1:SeriesOrder)=sign(eH(1:SeriesOrder)+dt*(ts).*cH(tt,1:SeriesOrder)).*sqrt(abs(sign(eH(1:SeriesOrder)).*eH(1:SeriesOrder).^2+sign(dt*(ts).*cH(tt,1:SeriesOrder)).*(dt).^2.*ts.^2.*cH(tt,1:SeriesOrder).^2));   
%   end
  
 %e(5)=1/3*e(5);
 
  yydt(1:Nn)=e0;
 for nn=1:SeriesOrdery
     yydt(1:Nn)=yydt(1:Nn)+e(nn)*Z(1:Nn).^nn;
 end

 
 %c0
 %c
 
 
 %str=input('Look at numbers')
 yy0(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));

 yy=yy0;
Dfyy(wnStart:Nn)=0;
Dfyydt(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    
    Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    Dfyydt(nn) = (yydt(nn + 1) - yydt(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;
pyydt(1:Nn)=0;
for nn = wnStart+1:Nn-1
    
    pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
    pyydt(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyydt(nn));
end

toc

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


theta1=1;
rng(29079137, 'twister')
paths=50000;
YY(1:paths)=x0;  %Original process monte carlo.
YYdt(1:paths)=0;   

Random1(1:paths)=0;
for tt=1:TtM
    
    
    Random1=randn(size(Random1));
    HermiteP1(1,1:paths)=1;
    HermiteP1(2,1:paths)=Random1(1:paths);
    HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
    HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
    HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;


    YYdt(1:paths)=YYdt(1:paths)+.5*YY(1:paths)*dtM;   

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

    
end




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


MaxCutOff=30;
NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too 
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
plot(yy(wnStart+1:Nn-1),pyy(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 %plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');
 
%title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
title(sprintf('x0 = %.4f,mu1=%.3f,beta1=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,mu1,beta1,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,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.');

NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa))/2;
[YdtDensity,IndexOutYdt,IndexMaxYdt] = MakeDensityFromSimulation_Infiniti_NEW(YYdt,paths,NoOfBins,MaxCutOff );
plot(yydt(wnStart+1:Nn-1),pyydt(wnStart+1:Nn-1),'r',IndexOutYdt(1:IndexMaxYdt),YdtDensity(1:IndexMaxYdt),'g');

%title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
title(sprintf('x0 = %.4f,kappa=%.3f,theta=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,kappa,theta,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));


legend({'Ito-Hermite Path Integral Density','Monte Carlo Path Integral Density'},'Location','northeast')
 
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');



end