Clear[t,y,Zi,Wi,Yi,Y,y0,i,i1,p,n,ZAns];
n=1;
Array[y,n+1,0];
y[1]:= y[0]^2;
Array[y0,n,0];
y0[0]=1;
p=8;
Y:=y0[0];
For[k=1,k<n,k++,Y=Y+y0[k] *t^k/k!];
Array[Zi,n+p-1,0];
Array[Wi,n+p-1,0];
Array[Yi,n+p-1,0];
Array[ti,n+p-1,0];
Zi[n]=(y[n]/.t-> ti[n]);
For [i=n,i<p+n-1,i++,(Zi[i+1]=0;For [k=0,k<n,k++,(Zi[i+1]=Zi[i+1]+D[Zi[i],y[k]]*(y[k+1]/.t-> ti[i+1]));]);];
For[i=n,i<= p+n-1,i++,(Wi[i]=Zi[i];For[k=0,k<n,k++,(Wi[i]=(Wi[i]/.y[k]-> y0[k]))];)];
For[i=n,i<= p+n-1,i++,(Yi[i]=Wi[i];For[i1=i,i1>=1,i1--,(Yi[i1-1]=Integrate[Yi[i1],{ti[i1],0,ti[i1-1]}])];Y=Y+(Yi[0]/.ti[0]-> t);)]
ZAns=Collect[Y,t,Simplify]//PolynomialForm[#,TraditionalOrder -> False]&

Clear[t,y,Zi,Wi,Yi,Y,y0,i,i1,p,n,ZAns];   This line is handy to clear variables especially when you want to solve for more than one ODE.
n=1;                                      n is order of the ODE. For 1st order ODE n=1; 
Array[y,n+1,0];                           Define a workhorse variable.
y[1]:= y[0]^2;                            Specify the ODE. y[1]= dy/dt. y[0]=y[t] and similarly for higher order ODE's y[2]=d2y/dt2
Array[y0,n,0];
y0[0]=1;                                  Specify initial conditions. y0[0]=y(0) similalry for 2nd order ODE, you will also have to specify y0[1]=dy(0)/dt and so on.
p=8;                                       Number of powers in expansion of series.
Y:=y0[0];                                  Allocate first initial condition.
For[k=1,k<n,k++,Y=Y+y0[k] *t^k/k!];       Allocation of higher order initial conditions.
Array[Zi,n+p-1,0];
Array[Wi,n+p-1,0];
Array[Yi,n+p-1,0];
Array[ti,n+p-1,0];
Zi[n]=(y[n]/.t-> ti[n]);
For [i=n,i<p+n-1,i++,(Zi[i+1]=0;For [k=0,k<n,k++,(Zi[i+1]=Zi[i+1]+D[Zi[i],y[k]]*(y[k+1]/.t-> ti[i+1]));]);];
For[i=n,i<= p+n-1,i++,(Wi[i]=Zi[i];For[k=0,k<n,k++,(Wi[i]=(Wi[i]/.y[k]-> y0[k]))];)];
For[i=n,i<= p+n-1,i++,(Yi[i]=Wi[i];For[i1=i,i1>=1,i1--,(Yi[i1-1]=Integrate[Yi[i1],{ti[i1],0,ti[i1-1]}])];Y=Y+(Yi[0]/.ti[0]-> t);)];     This is integration loop which iteratively integrates starting from zero. replace zero with initial time when it is different from zero.
ZAns=Collect[Y,t,Simplify]//PolynomialForm[#,TraditionalOrder -> False]&

function [] = FPESeriesCoeffsFromMomentsAndVarianceAdditionAdaptStep()

%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 simulated the transformed 
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coordinates.


dt=.125/2/2*1;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.
Tt=128/2/2*2;%*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:3)=dt/4;

%Below order for monte carlo
OrderA=4;  %
OrderM=4;  %

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



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

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


%StepSize Management. decrease the step size if MomentsIncreaseRatio is larger
    %than  MaxMomentRatio and increase the step size if 
    %MomentsIncreaseRatio is smaller than  MinMomentRatio
    %Look at the code in the body and feel free to play with stepsize
    %variables
MaxMomentRatio=1.01;
MinMomentRatio=1.002;
MaxStepSize=1/128;
MinStepSize=1/512/16;



%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);


NnMid=((1+Nn)/2)*dNn;
Z(1:Nn)=(((1:Nn)*dNn)-NnMid);
%Z(1:Nn)=(((1:Nn)-20.5)*dNn-NnMid);
Z
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);
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);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This block is analytics for higher order monte carlo. This is explained in
%my thread on wilmott.com

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Below analytics to calculate the Z-series of SDE at different time steps. 
%We use the principle of addition of cumulants to find the resulting
%density of SDE at each step.
ExpnOrder=7;
SeriesOrder=7;
NMoments=8;
wnStart=1;
tic
%Tx is running time.
%T is terminal time.
%tt is time index.
Tx=0;
tt=0;
while(Tx<T)
    tt=tt+1;
   % t1=(tt-1)*dt;
   % t2=(tt)*dt;
    if(tt==1)
 
        [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),8);
 
        c0=wMu0dt+w0;
        c(1)=sigma0*sqrt(ds(tt))
        c(2:10)=0.0;
        w0=c0;
        
        Tx=Tx+ds(tt);
    
    elseif(tt==2) 

    %dt=ds(tt);    
        
    %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),8);
    [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),8);
     [b0,b] = CalculateDriftbCoeffs08(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);
     c(6)=b(6);
     c(7)=b(7);
     
     Tx=Tx+ds(tt);       
            
    else
        
    
    %below wMu0dta is drift and dwMu0dtdwa are its derivatives.
    %wMu1dt is its first derivative and dwMu1dtdw are derivatives of first
    %derivative.


     %[wMu0dta,dwMu0dtdwa,wMu1dt,dwMu1dtdw] = BesselDriftAndDerivatives08(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,10);
     
     %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),8);
     [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),8);
     [b0,b] = CalculateDriftbCoeffs08(wMu0dt,dwMu0dtdw,c,SeriesOrder);

%Below Calculate original terms(associated with first hermite and second
%hermite polynomial in volatility)and its eight derivatives    
    [wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),8);
    [wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),8);

%Below Calculate Z-Series expansions of volatility terms associated with 
%first hermite polynomial and 2nd hermite polynomial.
    [g0,g] = CalculateDriftbCoeffs08(wVol0dt,dwVol0dtdw,c,SeriesOrder);
    [h0,h] = CalculateDriftbCoeffs08(wVol2dt,dwVol2dtdw,c,SeriesOrder);
%Principal term in first hermite was not included in expansion. It is 
%added here
    g0=g0+sigma0*sqrt(ds(tt));
    
    %Below function calculates Moments and Cumulants of diffusion term of
    %the SDE
    [Moments,kk] = CalculateVolTermMoments(g0,g,h0,h,SeriesOrder,NMoments);
    
    
 %Below drift series is linearly added to SDE series to get new intermediate 
 %SDE series.     
    c0=c0+b0;
    c(1:7)=c(1:7)+b(1:7);
    
 %Below calculate cumulants and moments of intermediate SDE series   
    [CC] = CalculateCumulants(c0,c,SeriesOrder,NMoments);
    [MM1] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
    
    
    %Below Calculate Cumulants of target density by arithmetic addition of
    %first eight cumulants of diffusion term and intermediate SDE term.
    k1=CC+kk;
    
    %Below Calculate Standardized Cumulants of the target density. 
    
    k(1)=0;
    k(2)=1;
    k(3)=k1(3)/k1(2).^1.5;
    k(4)=k1(4)/k1(2).^2.0;
    k(5)=k1(5)/k1(2).^2.5;
    k(6)=k1(6)/k1(2).^3.0;
    k(7)=k1(7)/k1(2).^3.5;
    k(8)=k1(8)/k1(2).^4.0;
    
    
    %Below Calculate Standardized Moments of target density.
    rmu(1)=k(1);
    rmu(2)=k(1)*rmu(1)+k(2);
    rmu(3)=k(1)*rmu(2)+2*k(2)*rmu(1)+k(3);
    rmu(4)=k(1)*rmu(3)+3*k(2)*rmu(2)+3*k(3)*rmu(1)+k(4);    
    rmu(5)=k(1)*rmu(4)+4*k(2)*rmu(3)+6*k(3)*rmu(2)+4*k(4)*rmu(1)+k(5);
    rmu(6)=k(1)*rmu(5)+5*k(2)*rmu(4)+10*k(3)*rmu(3)+10*k(4)*rmu(2)+5*k(5)*rmu(1)+k(6);
    rmu(7)=k(1)*rmu(6)+6*k(2)*rmu(5)+15*k(3)*rmu(4)+20*k(4)*rmu(3)+15*k(5)*rmu(2)+6*k(6)*rmu(1)+k(7);
    rmu(8)=k(1)*rmu(7)+7*k(2)*rmu(6)+21*k(3)*rmu(5)+35*k(4)*rmu(4)+35*k(5)*rmu(3)+21*k(6)*rmu(2)+7*k(7)*rmu(1)+k(8);

    %We provide initial guess for calibration procedure
    
    if(tt==3)
        a0Guess=(c0-k1(1))/sqrt(CC(2));
        aGuess(1)=c(1)/sqrt(CC(2));
        aGuess(2:7)=c(2:7)/sqrt(CC(2));
    else
        a0Guess=(c0-k1(1))/sqrt(CC(2));
        aGuess(1:7)=(c(1:7))/sqrt(CC(2));
    end
        
    %We calculate coefficients of Standardized target density
    [c0,c] = CalculateDensityFromStandardizedMoments(rmu,a0Guess,aGuess);
    
    %From standardized Z-seriesCoefficients of target density, we change to
    %actual non-standardized Z-series coefficients.
    c0=c0*sqrt(k1(2))+k1(1);
    c=c*sqrt(k1(2));
 
    %Actual Raw Moments after the Calculation of target density.
    [MM2] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
    %Calculate the Ratio of Moments between intermediate density and the
    %new target density
    MomentsIncreaseRatio=max(MM2(2:8)./MM1(2:8));
    
    %Below decrease the step size if MomentsIncreaseRatio is larger
    %than  MaxMomentRatio and increase the step size if 
    %MomentsIncreaseRatio is smaller than  MinMomentRatio

    IncreaseRatioFlag=0;
    DecreaseRatioFlag=0;
    if(MomentsIncreaseRatio>MaxMomentRatio)
        if (ds(tt)>MinStepSize)
            DecreaseRatioFlag=1;
        end
    elseif (MomentsIncreaseRatio<MinMomentRatio)   
        if (ds(tt)<MaxStepSize)
            IncreaseRatioFlag=1;
        end
    end
    
    if(DecreaseRatioFlag==1)
        ds(tt+1)=ds(tt)/2;
    
    elseif(IncreaseRatioFlag==1)
        ds(tt+1)=ds(tt)*2;
    else
        ds(tt+1)=ds(tt);
    end    
    
    %Tx is running time.
    Tx=Tx+ds(tt);
    
    %At last step make sure that step size is appropriate.
    if(Tx+ds(tt+1)>T)
        ds(tt+1)=T-Tx;
    end
    

    %Assign the median of density to w0
    w0=c0;
    
    end

end

ttFirst=tt;
wnStart=1;


%Below w which is the SDE in bessel coordinates and it is calulated from
%Z-series in Bessel coordinates.
w(1:Nn)=c0;
for nn=1:ExpnOrder
    w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^nn;
end
Flag=0;
for nn=Nn-1:-1:1
    if((w(nn)<0)&&(Flag==0))
        wnStart=nn+1;
        Flag=1;
    end
end


%Below we directly transform the density of SDE in bessel coordinates into
%density of SDE in original coordinates.

yy(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
%yy0=((1-gamma).*w0).^(1/(1-gamma));

%Below we make calculations of density in original coordinates by doing 
%change of probability derivative with respect to gaussian density.
Dfyy(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    
    Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    
    %Change of variable derivative for densities
end
Dfyy(Nn)=Dfyy(Nn-1);
Dfyy(1)=Dfyy(2);


pyy(1:Nn)=0;
for nn = wnStart:Nn
    
    pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
end


toc

%Below we do calculations for change of Z-series coefficients from
%coordinates in Bessel to Z-series coefficients in original coordinates of
%the SDE.
%The new series is expanded around c0 which is median of the density in
%Bessel coordinates.
   Bmu=c0;%+c(2)+3*c(4)+15*c(6)
   Bmu
   
%Calculate the first seven derivatives of transformation from Bessel
%coordinates to original coordinates of SDE at the median value Bmu=c0.
   
  yB0=((1-gamma).*Bmu).^(1/(1-gamma));
  dyB0(1)=((1-gamma).*Bmu).^(1/(1-gamma)-1);
  dyB0(2)=(1/(1-gamma)-1).*((1-gamma).*Bmu).^(1/(1-gamma)-2)*(1-gamma);
  dyB0(3)=(1/(1-gamma)-1).*(1/(1-gamma)-2).*((1-gamma).*Bmu).^(1/(1-gamma)-3)*(1-gamma)^2;
  dyB0(4)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*((1-gamma).*Bmu).^(1/(1-gamma)-4)*(1-gamma)^3;
  dyB0(5)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*((1-gamma).*Bmu).^(1/(1-gamma)-5)*(1-gamma)^4;
  dyB0(6)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*((1-gamma).*Bmu).^(1/(1-gamma)-6)*(1-gamma)^5;
  dyB0(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).*Bmu).^(1/(1-gamma)-7)*(1-gamma)^6;
%  dyB0(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).*Bmu).^(1/(1-gamma)-8)*(1-gamma)^7;
%  dyB0(9)=(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/(1-gamma)-8)*((1-gamma).*Bmu).^(1/(1-gamma)-9)*(1-gamma)^8;
%  dyB0(10)=(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/(1-gamma)-8)*(1/(1-gamma)-9)*((1-gamma).*Bmu).^(1/(1-gamma)-10)*(1-gamma)^9;

c(8:10)=0;

%Below finally calculate the Z-series coefficients from Bessel coordinates 
%to original coordinates using derivatives calculated at previous step.
   [d0,d] = CalculateDriftbCoeffs08(yB0,dyB0,c,7);
 
%Below calculate SDE variable in original coordinates after expanding the 
%Z-series coefficients in original coordinates.   
Yw(1:Nn)=d0;
for nn=1:7
    Yw(1:Nn)=Yw(1:Nn)+d(nn)*Z(1:Nn).^nn;
end

%Below calculate the density from Z-series Coefficients in original
%coordinates. We want to see how exact they are since we would need to use
%them in analytics later when calculating integrals of variance for
%stochastic volatility models.

DfYw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    
    DfYw(nn) = (Yw(nn + 1) - Yw(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    
    %Change of variable derivative for densities
end
DfYw(Nn)=DfYw(Nn-1);
DfYw(1)=DfYw(2);


pYw(1:Nn)=0;
for nn = wnStart:Nn
    
    pYw(nn) = (normpdf(Z(nn),0, 1))/abs(DfYw(nn));
end

%AnalyticMean1 is calculated from direct transformation of density in Bessel coordinates
AnalyticMean1=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from transformed density
%AnalyticMean2 is calculated from Z-series coordinates of density in Original coordinates
AnalyticMean2=d0+d(2)+3*d(4)+15*d(6) %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
%yB0

%theta1=1;

%Monte Carlo starts here. For details see my thread on wilmott.com 
%where I gave details of higher order monte carlo.

rng(29079137, 'twister')
paths=50000;
YY(1:paths)=x0;  %Original process monte carlo.
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;

 
    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);
    
    
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 
 
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(Yw(wnStart+1:Nn-1),pYw(wnStart+1:Nn-1),'b',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');

 
dtAvg=T/ttFirst; 
 
title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dtAvg=%.5f,M=%.4f,TM=%.4f', x0,theta,kappa,gamma,sigma0,T,dtAvg,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
legend({'SDE density From Z-series in Original Coordinates','SDE density Transformed Directly From Bessel Coordinates','Monte Carlo Density'},'Location','northeast')
 




str=input('red line is density of SDE from Addition of Cumulants, green is monte carlo.');


plot((wnStart+1:Nn-1)*dNn-NnMid,Yw(wnStart+1:Nn-1),'b',(wnStart+1:Nn-1)*dNn-NnMid,yy(wnStart+1:Nn-1),'r');

title('SDE Stochastic Variable as a Function of Standard Normal');%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
legend({'SDE Variable From Z-series in Original Coordinates','SDE Variable Transformed Directly From Bessel Coordinates'},'Location','northwest')

str=input('red line is density of SDE from addition of Cumulants method, green is monte carlo.');



end

function [Y] = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma)


%In the coefficient calculation program which calculates Y(l1,l2,l3,l4), 
%I have used four levels of looping each for relevant expansion order. 
%The first loop takes four values and second loop takes 16 values and 
%third loop takes 64 values and so on. And then each coefficient 
%term can be individually calculated while carefully accounting 
%for path dependence. 
%So for example in a nested loop structure 
%m1= 1:mDim
% m2=1:mDim
% m3=1:mDim
%    l(m1)=l(m1)+1;
%    l(m2)=l(m2)+1;
%    l(m3)=l(m3)+1;

%in the above looping loop takes values from one to four with one 
%indicating the first drift term, two indicating the second drift term 
%and three indicating quadratic variation term and 
%four indicating the volatility term. And with this looping structure 
%we can So in the above looping m1=1 would mean that all terms are 
%descendent of first drift term and m2=4 would mean that all terms are 
%descendent of first drift term on first expansion order and descendent 
%of volatility term on the second order and so we can keep track of path 
%dependence perfectly. 
%And on each level, we individually calculate the descendent terms. While 
%keeping track of path dependence and calculating the coefficients with 
%careful path dependence consideration, we update the appropriate element 
%in our polynomial like expansion coefficient array 

%explaining the part of code
%m1= 1:mDim
% m2=1:mDim
% m3=1:mDim
%    l(m1)=l(m1)+1;
%    l(m2)=l(m2)+1;
%    l(m3)=l(m3)+1;
%Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2);

%Here l(1) denotes l1 but written as l(1) so it can be conveniently 
%updated with the loop variable when the loop variable takes value one 
%indicating first drift term . And l(2) could be conveniently updated when 
%the loop variable takes value equal to two indicating second 
%drift term and so on.
%Here is the part of code snippet for that

%for m1=1:mDim
%    l(1)=1;
%    l(2)=1;
%    l(3)=1;
%    l(4)=1;
%    l(m1)=l(m1)+1;
%CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
%    CoeffDX2 = CoeffDX1 - 1;
%    ArrIndex0=m1;
%    ArrIndex=(m1-1)*mDim;
%    Coeff1st=Y1(ArrIndex0)*CoeffDX1;
%    Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
%    Y2(1+ArrIndex)=Coeff1st;
%    Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
%    Y2(2+ArrIndex)=Coeff1st;
%    Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
%    Y2(3+ArrIndex)=Coeff2nd;
%    Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,2,n1(m1));
%    Y2(4+ArrIndex)=Coeff1st;
%    Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,3,n1(m1));

 
%The first four lines update the array indices according to the parent term. 
%And then CoeffDX1 and CoeffDX2 are calculated according to algebraic exponents on parent terms.    
%ArrIndex0=m1; calculates the array index of the parent term 
%And ArrIndex=(m1-1)*mDim; calculates the array index of the descendent terms
%And coefficient of the drift and volatility descendent terms is 
%calculated by multiplying the coefficient of the parent term by 
%Coeff1st=Y1(ArrIndex0)*CoeffDX1;
%And coefficient of the quadratic variation descendent terms is 
%calculated by multiplying the coefficient of the parent term by 
%Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
%And then each of the four descendent terms are updated with Coeff1st 
%if they are drift or volatility descendent terms or Coeff2nd if 
%they are quadratic variation descendent terms.
%Here Y1 indicates the temporary coefficient array with parent terms on 
%first level. And Y2 denotes temporary coefficient array with parent terms 
%on second level and Y3 indicates temporary coefficient array with parent terms on third level.

[IntegralCoeff,IntegralCoeffdt,IntegralCoeffdz] = ComputeIntegralCoeffs();

n1(1)=2;
n1(2)=2;
n1(3)=2;
n1(4)=3;
%n1(1), n1(2), n1(3) are drift and quadratic variation term variables 
%and take a value equal to 2 indicating a dt integral.
%n1(4) is volatility term variable and indicates a dz-integral by taking a
%value of 3. 

mDim=4; % four descendent terms in each expansion
Y(1:5,1:5,1:5,1:5)=0;
Y1(1:mDim)=0;
Y2(1:mDim*mDim)=0;
Y3(1:mDim*mDim*mDim)=0;


%First Ito-hermite expansion level starts here. No loops but four
%descendent terms.
l(1)=1;
l(2)=1;
l(3)=1;
l(4)=1;
CoeffDX1 = alpha;
CoeffDX2 = CoeffDX1 - 1;
Coeff1st=CoeffDX1;
Coeff2nd=.5*CoeffDX1*CoeffDX2;

Y1(1)=Coeff1st;
Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2);
Y1(2)=Coeff1st;
Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2);
Y1(3)=Coeff2nd;
Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,2);
Y1(4)=Coeff1st;
Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,3);

%Second Ito-hermite expansion level starts. It has a loop over four parent
%terms and there are four descendent terms for each parent term.
%The coefficient terms are then lumped in a polynomial-like expansion
%array of coefficients.
for m1=1:mDim
    l(1)=1;
    l(2)=1;
    l(3)=1;
    l(4)=1;
    l(m1)=l(m1)+1;
    CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
    CoeffDX2 = CoeffDX1 - 1;
    ArrIndex0=m1;
    ArrIndex=(m1-1)*mDim;
    Coeff1st=Y1(ArrIndex0)*CoeffDX1;
    Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
    Y2(1+ArrIndex)=Coeff1st;
    Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
    Y2(2+ArrIndex)=Coeff1st;
    Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
    Y2(3+ArrIndex)=Coeff2nd;
    Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,2,n1(m1));
    Y2(4+ArrIndex)=Coeff1st;
    Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,3,n1(m1));
    %Third Ito-hermite expansion level starts and it is a nested loop with
    %a total of sixteen parents and each parent takes four descendent
    %terms.
    for m2=1:mDim
        l(1)=1;
        l(2)=1;
        l(3)=1;
        l(4)=1;
        l(m1)=l(m1)+1;
        l(m2)=l(m2)+1;
        CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
        CoeffDX2=CoeffDX1-1;            
        ArrIndex0=(m1-1)*mDim+m2;
        ArrIndex=((m1-1)*mDim+(m2-1))*mDim;
        Coeff1st=Y2(ArrIndex0)*CoeffDX1;
        Coeff2nd=Y2(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
        Y3(1+ArrIndex)=Coeff1st;
        Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,2,n1(m2),n1(m1));
        Y3(2+ArrIndex)=Coeff1st;
        Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,2,n1(m2),n1(m1));
        Y3(3+ArrIndex)=Coeff2nd;
        Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,2,n1(m2),n1(m1));
        Y3(4+ArrIndex)=Coeff1st;
        Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,3,n1(m2),n1(m1));
        %fourht Ito-hermite expansion level starts and it is a triply-nested loop with
        %a total of sixteen parents and each parent takes four descendent
        %terms. We then lump the terms in a relatively sparse polynomial 
        %like expansion coefficient array that has smaller number of
        %non-zero terms.

        for m3=1:mDim
            l(1)=1;
            l(2)=1;
            l(3)=1;
            l(4)=1;
            l(m1)=l(m1)+1;
            l(m2)=l(m2)+1;
            l(m3)=l(m3)+1;
            CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
            CoeffDX2=CoeffDX1-1;
            ArrIndex0=((m1-1)*mDim+(m2-1))*mDim+m3;
            %ArrIndex=(((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim;
            Coeff1st=Y3(ArrIndex0)*CoeffDX1;
            Coeff2nd=Y3(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
            %Y4(1+ArrIndex)=Coeff1st;
            Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(2,n1(m3),n1(m2),n1(m1));
            %Y4(2+ArrIndex)=Coeff1st;
            Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(2,n1(m3),n1(m2),n1(m1));
            %Y4(3+ArrIndex)=Coeff2nd;
            Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(2,n1(m3),n1(m2),n1(m1));
            %Y4(4+ArrIndex)=Coeff1st;
            Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(3,n1(m3),n1(m2),n1(m1));
         end
    end
end
        
        
        

end

function [IntegralCoeff0,IntegralCoeffdt,IntegralCoeffdz] = ComputeIntegralCoeffs()


IntegralCoeff(1:3,1:3,1:3,1:3)=0;
IntegralCoeffdt(1:3,1:3,1:3,1:3)=0;
IntegralCoeffdz(1:3,1:3,1:3,1:3)=0;

IntegralCoeff0(1:3,1:3,1:3,1:3)=0;

IntegralCoeff0(1,1,1,2)=1;
IntegralCoeff0(1,1,1,3)=1;

IntegralCoeff0(1,1,2,2)=1/2;
IntegralCoeff0(1,1,3,2)=1-1/sqrt(3);
IntegralCoeff0(1,1,2,3)=1/sqrt(3);
IntegralCoeff0(1,1,3,3)=1/2;

IntegralCoeff0(1,2,2,2)=1/6;
IntegralCoeff0(1,3,2,2)=(1-1/sqrt(3))*1/2*(1-1/sqrt(5));
IntegralCoeff0(1,2,3,2)=1/sqrt(3)/2*(1-1/sqrt(5));
IntegralCoeff0(1,3,3,2)=1/2*(1-sqrt(2)/2);
IntegralCoeff0(1,2,2,3)=1/2*1/sqrt(5);
IntegralCoeff0(1,3,2,3)=(1-1/sqrt(3))*1/sqrt(2)*1/2;
IntegralCoeff0(1,2,3,3)=1/sqrt(3)/sqrt(2)*1/(2);
IntegralCoeff0(1,3,3,3)=1/6;

IntegralCoeff0(2,2,2,2)=1/24;
IntegralCoeff0(2,2,3,2)=1/2*1/sqrt(5)*1/3*(1-1/sqrt(7));
IntegralCoeff0(2,3,2,2)=1/sqrt(3)*1/2*(1-1/sqrt(5))*1/3*(1-1/sqrt(7));
IntegralCoeff0(3,2,2,2)=(1-1/sqrt(3))*1/2*(1-1/sqrt(5))*1/3*(1-1/sqrt(7));
IntegralCoeff0(3,3,2,2)=1/2*(1-sqrt(2)/2)*1/2*(1-sqrt(2)/sqrt(6));
IntegralCoeff0(2,3,3,2)=1/sqrt(3)*1/sqrt(2)*1/2*1/2*(1-sqrt(2)/sqrt(6));
IntegralCoeff0(3,2,3,2)=(1-1/sqrt(3))*1/sqrt(2)*1/2*1/2*(1-sqrt(2)/sqrt(6));
IntegralCoeff0(3,3,3,2)=1/6*(1-sqrt(3)/sqrt(5));

IntegralCoeff0(2,2,2,3)=1/6*1/sqrt(7);
IntegralCoeff0(2,2,3,3)=1/2*1/sqrt(5)*1/sqrt(2)*1/sqrt(6);
IntegralCoeff0(2,3,2,3)=1/sqrt(3)*1/2*(1-1/sqrt(5))*1/sqrt(2)*1/sqrt(6);
IntegralCoeff0(3,2,2,3)=(1-1/sqrt(3))*1/2*(1-1/sqrt(5))*1/sqrt(2)*1/sqrt(6);
IntegralCoeff0(3,3,2,3)=1/2*(1-sqrt(2)/2)*1/sqrt(3)*1/sqrt(5);
IntegralCoeff0(2,3,3,3)=1/sqrt(3)*1/sqrt(2)*1/sqrt(4)*1/sqrt(3)*(1/sqrt(5));
IntegralCoeff0(3,2,3,3)=(1-1/sqrt(3))*1/sqrt(2)*1/2*1/sqrt(3)*1/sqrt(5);
IntegralCoeff0(3,3,3,3)=1/24;


%Can also be calculated with algorithm below.
%here IntegralCoeffdt indicates the coefficients of a dt-integral.
%This dt-integral means that you are calculating the Ito-hermite expansion
%of  Integral X(t)^alpha dt with X(t) dynamics given by the SDE 
%here IntegralCoeffdz indicates the coefficients of a dz-integral.
%This dz-integral means that you are calculating the Ito-hermite expansion
%of  Integral X(t)^alpha dz(t) with X(t) dynamics given by the SDE 

%IntegralCoeff is is associated with expansion of X(t)^alpha.
%IntegralCoeff below is not returned and IntegralCoeff0 manually calculated
%above is returned but both are the same.

l0(1:2)=1;

for m4=1:2
    l0(1)=1;
    l0(2)=1;
    
    %IntegralCoeff4(m4,1,1,1)=1;
    %IntegralCoeff4(m4,1,1,1)=1;
    %1 is added to m4 since index 1 stands for zero, 2 for one and three
    %for two.
    IntegralCoeff(1,1,1,m4+1)=1;
    l0(m4)=l0(m4)+1;
    IntegralCoeffdt(1,1,1,m4+1)=IntegralCoeff(1,1,1,m4+1)* ...
        1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
    IntegralCoeffdz(1,1,1,m4+1)= IntegralCoeff(1,1,1,m4+1)* ...
        1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
    for m3=1:2
        l0(1)=1;
        l0(2)=1;
        l0(m4)=l0(m4)+1;
        
        if(m3==1)
            IntegralCoeff(1,1,m4+1,m3+1)=1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
        end
        if(m3==2)
            IntegralCoeff(1,1,m4+1,m3+1)= 1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
        end
        l0(m3)=l0(m3)+1;
        %IntegralCoeff(1,1,m4+1,m3+1)=IntegralCoeff4(m4,m3,1,1);
        IntegralCoeffdt(1,1,m4+1,m3+1)=IntegralCoeff(1,1,m4+1,m3+1)* ...
            1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
        IntegralCoeffdz(1,1,m4+1,m3+1)= IntegralCoeff(1,1,m4+1,m3+1)* ...
            1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
        for m2=1:2
            l0(1)=1;
            l0(2)=1;
            l0(m4)=l0(m4)+1;
            l0(m3)=l0(m3)+1;
            
            if(m2==1)
                IntegralCoeff(1,m4+1,m3+1,m2+1)=IntegralCoeff(1,1,m4+1,m3+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
            end
            if(m2==2)
                IntegralCoeff(1,m4+1,m3+1,m2+1)= IntegralCoeff(1,1,m4+1,m3+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
            end
            l0(m2)=l0(m2)+1;
            %IntegralCoeff(1,m4+1,m3+1,m2+1)=IntegralCoeff4(m4,m3,m2,1);
            IntegralCoeffdt(1,m4+1,m3+1,m2+1)=IntegralCoeff(1,m4+1,m3+1,m2+1)* ...
                1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
            IntegralCoeffdz(1,m4+1,m3+1,m2+1)= IntegralCoeff(1,m4+1,m3+1,m2+1)* ...
                1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);

            for m1=1:2
                l0(1)=1;
                l0(2)=1;
                l0(m4)=l0(m4)+1;
                l0(m3)=l0(m3)+1;
                l0(m2)=l0(m2)+1;
                if(m1==1)
                    IntegralCoeff(m4+1,m3+1,m2+1,m1+1)=IntegralCoeff(1,m4+1,m3+1,m2+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                end
                if(m1==2)
                    IntegralCoeff(m4+1,m3+1,m2+1,m1+1)= IntegralCoeff(1,m4+1,m3+1,m2+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                end
                l0(m1)=l0(m1)+1;
                %IntegralCoeff(m4+1,m3+1,m2+1,m1+1)=IntegralCoeff4(m4,m3,m2,m1);
                IntegralCoeffdt(m4+1,m3+1,m2+1,m1+1)=IntegralCoeff(m4+1,m3+1,m2+1,m1+1)* ...
                    1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                IntegralCoeffdz(m4+1,m3+1,m2+1,m1+1)= IntegralCoeff(m4+1,m3+1,m2+1,m1+1)* ...
                    1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
            end
        end
    end
end

end

function [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder)


NoOfTerms=19;
%NoOfTerms=9;

YqCoeffa(1:NoOfTerms)=0.0;
%Fp1
Fp1=Fp1/(1-gamma);
%Fp1
%gamma
YqCoeffa(1)=YqCoeff0(1,1,2,1)*dt;%*(1-gamma)^Fp1(1,1,2,1)*dt;
YqCoeffa(2)=YqCoeff0(1,2,1,1)*dt;%*(1-gamma)^Fp1(1,2,1,1)*dt;
YqCoeffa(3)=YqCoeff0(2,1,1,1)*dt;%*(1-gamma)^Fp1(2,1,1,1)*dt;
YqCoeffa(4)=YqCoeff0(1,1,3,1)*dt^2;%*(1-gamma)^Fp1(1,1,3,1)*dt^2;
YqCoeffa(5)=YqCoeff0(1,2,2,1)*dt^2;%*(1-gamma)^Fp1(1,2,2,1)*dt^2;
YqCoeffa(6)=YqCoeff0(2,1,2,1)*dt^2;%*(1-gamma)^Fp1(2,1,2,1)*dt^2;
YqCoeffa(7)=YqCoeff0(1,3,1,1)*dt^2;%*(1-gamma)^Fp1(1,3,1,1)*dt^2;
YqCoeffa(8)=YqCoeff0(2,2,1,1)*dt^2;%*(1-gamma)^Fp1(2,2,1,1)*dt^2;
YqCoeffa(9)=YqCoeff0(3,1,1,1)*dt^2;%*(1-gamma)^Fp1(3,1,1,1)*dt^2;
YqCoeffa(10)=YqCoeff0(1,1,4,1)*dt^3;%;*(1-gamma)^Fp1(1,1,4,1)*dt^3;
YqCoeffa(11)=YqCoeff0(1,2,3,1)*dt^3;%;*(1-gamma)^Fp1(1,2,3,1)*dt^3;
YqCoeffa(12)=YqCoeff0(2,1,3,1)*dt^3;%*(1-gamma)^Fp1(2,1,3,1)*dt^3;
YqCoeffa(13)=YqCoeff0(1,3,2,1)*dt^3;%*(1-gamma)^Fp1(1,3,2,1)*dt^3;
YqCoeffa(14)=YqCoeff0(2,2,2,1)*dt^3;%*(1-gamma)^Fp1(2,2,2,1)*dt^3;
YqCoeffa(15)=YqCoeff0(3,1,2,1)*dt^3;%*(1-gamma)^Fp1(3,1,2,1)*dt^3;
YqCoeffa(16)=YqCoeff0(1,4,1,1)*dt^3;%*(1-gamma)^Fp1(1,4,1,1)*dt^3;
YqCoeffa(17)=YqCoeff0(2,3,1,1)*dt^3;%*(1-gamma)^Fp1(2,3,1,1)*dt^3;
YqCoeffa(18)=YqCoeff0(3,2,1,1)*dt^3;%*(1-gamma)^Fp1(3,2,1,1)*dt^3;
YqCoeffa(19)=YqCoeff0(4,1,1,1)*dt^3;%*(1-gamma)^Fp1(4,1,1,1)*dt^3;

Fp2(1)=Fp1(1,1,2,1);
Fp2(2)=Fp1(1,2,1,1);
Fp2(3)=Fp1(2,1,1,1);
Fp2(4)=Fp1(1,1,3,1);
Fp2(5)=Fp1(1,2,2,1);
Fp2(6)=Fp1(2,1,2,1);
Fp2(7)=Fp1(1,3,1,1);
Fp2(8)=Fp1(2,2,1,1);
Fp2(9)=Fp1(3,1,1,1);
Fp2(10)=Fp1(1,1,4,1);
Fp2(11)=Fp1(1,2,3,1);
Fp2(12)=Fp1(2,1,3,1);
Fp2(13)=Fp1(1,3,2,1);
Fp2(14)=Fp1(2,2,2,1);
Fp2(15)=Fp1(3,1,2,1);
Fp2(16)=Fp1(1,4,1,1);
Fp2(17)=Fp1(2,3,1,1);
Fp2(18)=Fp1(3,2,1,1);
Fp2(19)=Fp1(4,1,1,1);

%YqCoeffa
%Fp2
%str=input('Look at numbers');
wMu0dt0=0;
dwMu0dt(1:ExpnOrder)=0.0;

wMu0dt=0;
dwMu0dtdw(1:ExpnOrder)=0.0;

for mm=1:NoOfTerms

    wMu0dt0=YqCoeffa(mm).*((1-gamma)*w0).^Fp2(mm);
    for nn=1:ExpnOrder
        if(nn==1)
            dwMu0dt(nn)=wMu0dt0*Fp2(mm)/w0;
        else
        dwMu0dt(nn)=dwMu0dt(nn-1)*(Fp2(mm)-(nn-1))/w0;
        end
    end
    wMu0dt=wMu0dt+wMu0dt0;
    for nn=1:ExpnOrder
        dwMu0dtdw(nn)=dwMu0dtdw(nn)+dwMu0dt(nn);
    end
        
    
end



end

function [b0,b] = CalculateDriftbCoeffs08(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

if(SeriesOrder>=9)
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));
end


if(SeriesOrder>=10)
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

end

function [wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder)




% c1(wnStart:Nn)=YqCoeff0(1,1,1,2).*yy(wnStart:Nn).^Fp1(1,1,1,2).*sqrt(dt)+ ...
%     (YqCoeff0(1,1,2,2).*yy(wnStart:Nn).^Fp1(1,1,2,2)+YqCoeff0(1,2,1,2).*yy(wnStart:Nn).^Fp1(1,2,1,2)+ ...
%     YqCoeff0(2,1,1,2).*yy(wnStart:Nn).^Fp1(2,1,1,2)).*dt^1.5+ ...
%     (YqCoeff0(1,1,3,2).*yy(wnStart:Nn).^Fp1(1,1,3,2)+YqCoeff0(1,2,2,2).*yy(wnStart:Nn).^Fp1(1,2,2,2)+ ...
%     YqCoeff0(2,1,2,2).*yy(wnStart:Nn).^Fp1(2,1,2,2)+YqCoeff0(1,3,1,2).*yy(wnStart:Nn).^Fp1(1,3,1,2)+ ...
%     YqCoeff0(2,2,1,2).*yy(wnStart:Nn).^Fp1(2,2,1,2)+YqCoeff0(3,1,1,2).*yy(wnStart:Nn).^Fp1(3,1,1,2)).*dt^2.5;

NoOfTerms=9;%excluding dt^3 terms
YqCoeffa(1:NoOfTerms)=0.0;

Fp1=Fp1/(1-gamma);
YqCoeffa(1)=YqCoeff0(1,1,2,2)*(1-gamma)^Fp1(1,1,2,2)*dt^1.5;
YqCoeffa(2)=YqCoeff0(1,2,1,2)*(1-gamma)^Fp1(1,2,1,2)*dt^1.5;
YqCoeffa(3)=YqCoeff0(2,1,1,2)*(1-gamma)^Fp1(2,1,1,2)*dt^1.5;
YqCoeffa(4)=YqCoeff0(1,1,3,2)*(1-gamma)^Fp1(1,1,3,2)*dt^2.5;
YqCoeffa(5)=YqCoeff0(1,2,2,2)*(1-gamma)^Fp1(1,2,2,2)*dt^2.5;
YqCoeffa(6)=YqCoeff0(2,1,2,2)*(1-gamma)^Fp1(2,1,2,2)*dt^2.5;
YqCoeffa(7)=YqCoeff0(1,3,1,2)*(1-gamma)^Fp1(1,3,1,2)*dt^2.5;
YqCoeffa(8)=YqCoeff0(2,2,1,2)*(1-gamma)^Fp1(2,2,1,2)*dt^2.5;
YqCoeffa(9)=YqCoeff0(3,1,1,2)*(1-gamma)^Fp1(3,1,1,2)*dt^2.5;

Fp2(1)=Fp1(1,1,2,2);
Fp2(2)=Fp1(1,2,1,2);
Fp2(3)=Fp1(2,1,1,2);
Fp2(4)=Fp1(1,1,3,2);
Fp2(5)=Fp1(1,2,2,2);
Fp2(6)=Fp1(2,1,2,2);
Fp2(7)=Fp1(1,3,1,2);
Fp2(8)=Fp1(2,2,1,2);
Fp2(9)=Fp1(3,1,1,2);


wVol0dt0=0;
dwVol0dt(1:ExpnOrder)=0.0;

wVol0dt=0;
dwVol0dtdw(1:ExpnOrder)=0.0;

for mm=1:NoOfTerms

    wVol0dt0=YqCoeffa(mm).*w0.^Fp2(mm);
    for nn=1:ExpnOrder
        if(nn==1)
            dwVol0dt(nn)=wVol0dt0*Fp2(mm)*1/w0;
        else
        dwVol0dt(nn)=dwVol0dt(nn-1)*(Fp2(mm)-(nn-1))/w0;
        end
    end
    wVol0dt=wVol0dt+wVol0dt0;
    for nn=1:ExpnOrder
        dwVol0dtdw(nn)=dwVol0dtdw(nn)+dwVol0dt(nn);
    end
        
    
end



end

function [wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder)


% c2(wnStart:Nn)=YqCoeff0(1,1,1,3).*yy(wnStart:Nn).^Fp1(1,1,1,3) *dt + ...
%     (YqCoeff0(1,1,2,3).*yy(wnStart:Nn).^Fp1(1,1,2,3)+YqCoeff0(1,2,1,3).*yy(wnStart:Nn).^Fp1(1,2,1,3)+ ...
%     YqCoeff0(2,1,1,3).*yy(wnStart:Nn).^Fp1(2,1,1,3)).*dt^2+ ...
%     (YqCoeff0(1,1,3,3).*yy(wnStart:Nn).^Fp1(1,1,3,3)+YqCoeff0(1,2,2,3).*yy(wnStart:Nn).^Fp1(1,2,2,3)+ ...
%     YqCoeff0(2,1,2,3).*yy(wnStart:Nn).^Fp1(2,1,2,3) + YqCoeff0(1,3,1,3).*yy(wnStart:Nn).^Fp1(1,3,1,3)+ ...
%     YqCoeff0(2,2,1,3).*yy(wnStart:Nn).^Fp1(2,2,1,3)+YqCoeff0(3,1,1,3).*yy(wnStart:Nn).^Fp1(3,1,1,3)).*dt^3;

Fp1=Fp1/(1-gamma);
 
NoOfTerms=9;%excluding dt^3 terms
YqCoeffa(1:NoOfTerms)=0.0;
YqCoeffa(1)=YqCoeff0(1,1,2,3)*(1-gamma)^Fp1(1,1,2,3)*dt^2;
YqCoeffa(2)=YqCoeff0(1,2,1,3)*(1-gamma)^Fp1(1,2,1,3)*dt^2;
YqCoeffa(3)=YqCoeff0(2,1,1,3)*(1-gamma)^Fp1(2,1,1,3)*dt^2;
YqCoeffa(4)=YqCoeff0(1,1,3,3)*(1-gamma)^Fp1(1,1,3,3)*dt^3;
YqCoeffa(5)=YqCoeff0(1,2,2,3)*(1-gamma)^Fp1(1,2,2,3)*dt^3;
YqCoeffa(6)=YqCoeff0(2,1,2,3)*(1-gamma)^Fp1(2,1,2,3)*dt^3;
YqCoeffa(7)=YqCoeff0(1,3,1,3)*(1-gamma)^Fp1(1,3,1,3)*dt^3;
YqCoeffa(8)=YqCoeff0(2,2,1,3)*(1-gamma)^Fp1(2,2,1,3)*dt^3;
YqCoeffa(9)=YqCoeff0(3,1,1,3)*(1-gamma)^Fp1(3,1,1,3)*dt^3;


Fp2(1)=Fp1(1,1,2,3);
Fp2(2)=Fp1(1,2,1,3);
Fp2(3)=Fp1(2,1,1,3);
Fp2(4)=Fp1(1,1,3,3);
Fp2(5)=Fp1(1,2,2,3);
Fp2(6)=Fp1(2,1,2,3);
Fp2(7)=Fp1(1,3,1,3);
Fp2(8)=Fp1(2,2,1,3);
Fp2(9)=Fp1(3,1,1,3);


wVol2dt0=0;
dwVol2dt(1:ExpnOrder)=0.0;

wVol2dt=0;
dwVol2dtdw(1:ExpnOrder)=0.0;

for mm=1:NoOfTerms

    wVol2dt0=YqCoeffa(mm).*w0.^Fp2(mm);
    for nn=1:ExpnOrder
        if(nn==1)
            dwVol2dt(nn)=wVol2dt0*Fp2(mm)*1/w0;
        else
        dwVol2dt(nn)=dwVol2dt(nn-1)*(Fp2(mm)-(nn-1))/w0;
        end
    end
    wVol2dt=wVol2dt+wVol2dt0;
    for nn=1:ExpnOrder
        dwVol2dtdw(nn)=dwVol2dtdw(nn)+dwVol2dt(nn);
    end
        
    
end



end

function [Moments,k] = CalculateVolTermMoments(a0,a,b0,b,SeriesOrder,NMoments)


[EA,EB,EAB] = CalculateSeriesProducts(a0,a,b0,b,SeriesOrder,NMoments);


c0=0;
c(1)=1;
c(2)=0;
d0=-1;
d(1)=0;
d(2)=1;

[EZ1,EZ2,EZ1Z2] = CalculateSeriesProducts(c0,c,d0,d,2,8);


Moments(1)=EA(1)*EZ1(1)+EB(1)*EZ2(1);
Moments(2)=EA(2)*EZ1(2)+2*EAB(1,1)*EZ1Z2(1,1)+EB(2)*EZ2(2);
Moments(3)=EA(3)*EZ1(3)+3*EAB(2,1)*EZ1Z2(2,1)+3*EAB(1,2)*EZ1Z2(1,2)+EB(3)*EZ2(3);
Moments(4)=EA(4)*EZ1(4)+4*EAB(3,1)*EZ1Z2(3,1)+6*EAB(2,2)*EZ1Z2(2,2)+4*EAB(1,3)*EZ1Z2(1,3)+EB(4)*EZ2(4);
Moments(5)=EA(5)*EZ1(5)+5*EAB(4,1)*EZ1Z2(4,1)+10*EAB(3,2)*EZ1Z2(3,2)+10*EAB(2,3)*EZ1Z2(2,3)+5*EAB(1,4)*EZ1Z2(1,4)+EB(5)*EZ2(5);
Moments(6)=EA(6)*EZ1(6)+6*EAB(5,1)*EZ1Z2(5,1)+15*EAB(4,2)*EZ1Z2(4,2)+20*EAB(3,3)*EZ1Z2(3,3)+15*EAB(2,4)*EZ1Z2(2,4)+6*EAB(1,5)*EZ1Z2(1,5)+EB(6)*EZ2(6);
Moments(7)=EA(7)*EZ1(7)+7*EAB(6,1)*EZ1Z2(6,1)+21*EAB(5,2)*EZ1Z2(5,2)+35*EAB(4,3)*EZ1Z2(4,3)+35*EAB(3,4)*EZ1Z2(3,4)+21*EAB(2,5)*EZ1Z2(2,5)+7*EAB(1,6)*EZ1Z2(1,6)+EB(7)*EZ2(7);
Moments(8)=EA(8)*EZ1(8)+8*EAB(7,1)*EZ1Z2(7,1)+28*EAB(6,2)*EZ1Z2(6,2)+56*EAB(5,3)*EZ1Z2(5,3)+70*EAB(4,4)*EZ1Z2(4,4)+56*EAB(3,5)*EZ1Z2(3,5)+28*EAB(2,6)*EZ1Z2(2,6)+8*EAB(1,7)*EZ1Z2(1,7)+EB(8)*EZ2(8);



k(1)=0;
k(2)=Moments(2);
k(3)=Moments(3);
k(4)=Moments(4)-3*Moments(2).^2;
k(5)=Moments(5)-10*Moments(3)*Moments(2);
k(6)=Moments(6)-15*Moments(4)*Moments(2)-10*Moments(3).^2+30*Moments(2).^3;
k(7)=Moments(7)- 21* Moments(2) *Moments(5)- 35* Moments(3) *Moments(4) + 210* Moments(2).^2* Moments(3);
k(8)=Moments(8)- 7 *k(2)* Moments(6) - 21* k(3)* Moments(5) - 35* k(4)* Moments(4) - ... 
  35* k(5)* Moments(3) - 21* k(6)* Moments(2);


end

%function [u1,u2,u3,u4,u5,u6] = CalculateCumulants(a0,a,SeriesOrder,NMoments)
function [CC] = CalculateCumulants(a0,a,SeriesOrder,NMoments)
%[CC] = CalculateCumulants(a0,a,SeriesOrder,NoOfCumulants);

%a(1)=a1;
%a(2)=a2;
%a(3)=a3;
%a(4)=a4;
%a(5)=a5;
%a(6)=a6;
%a(7)=a7;



aa0=a0;
a0=0;

EZ(1)=0;
EZ(2)=1;
for nn=3:NMoments*SeriesOrder+2
    if rem(nn,2)==1
        EZ(nn)=0;
    else
        EZ(nn)=EZ(nn-2)*(nn-1);
        EZ(nn);
    end
end
EZ;

EXZ(1,1)=1;
%for pp1=1:NZterms
%        EXZ(1,pp1+1)=EZ(pp1);
%end


a(SeriesOrder+1:60)=0;
b0=a0;
b=a;

for mm=1:NMoments
    if(mm>1)
        [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*mm);
        b(SeriesOrder*mm+1:60)=0;
    end
   % b0
   % b
%str=input('Look at numbers')    
    EXZ(mm+1,1)=b0;
    for pp2=1:SeriesOrder*mm
        
            EXZ(mm+1,1)=EXZ(mm+1,1)+b(pp2).*EZ(pp2);
    end
end

%u1=EXZ(2,1);

if(SeriesOrder>=4)
u1=a(2)+3*a(4);
end
if(SeriesOrder>=6)
u1=u1+15*a(6);
end


u2=EXZ(3,1);
u3=EXZ(4,1);
u4=EXZ(5,1);

if(SeriesOrder>=4)
k1=aa0+a(2)+3*a(4);
end
if(SeriesOrder>=6)
k1=k1+15*a(6);
end
k2=u2-u1^2;
k3=u3-3*u2*u1+2*u1^3;
k4=u4-4*u3*u1-3*u2^2+12*u2*u1^2-6*u1^4;

CC(1)=k1;
CC(2)=k2;
CC(3)=k3;
CC(4)=k4;

if(NMoments>=5)
u5=EXZ(6,1);
k5=u5-5*u4*u1-10*u3*u2+20*u3*u1^2+30*u2^2*u1-60*u2*u1^3+24*u1^5;
CC(5)=k5;
end
if(NMoments>=6)
u6=EXZ(7,1);
k6=u6-6*u5*u1-15*u4*u2+30*u4*u1^2-10*u3^2+120*u3*u2*u1-120*u3*u1^3+30*u2^3 - ...
    270*u2^2*u1^2+360*u2*u1^4-120*u1^6;
CC(6)=k6;
end
if(NMoments>=7)
u7=EXZ(8,1);
k7=u7-u1*u6-6*k2*u5-15*k3*u4-20*k4*u3-15*k5*u2-6*k6*u1;
CC(7)=k7;
end
if(NMoments>=8)
u8=EXZ(9,1);
k8 = u8 - u1* u7 - 7 *k2* u6 - 21* k3* u5 - 35* k4* u4 - ... 
  35* k5* u3 - 21* k6* u2 - 7 *k7* u1;
CC(8)=k8;
end


end

function [Moments] = CalculateMomentsOfZSeries(a0,a,SeriesOrder,NMoments)


%aa0=a0;
%a0=0;% ---1

EZ(1)=0;
EZ(2)=1;
%LogEZ(2)=0;
for nn=3:NMoments*SeriesOrder
    if rem(nn,2)==1
        EZ(nn)=0;
    else
        EZ(nn)=EZ(nn-2)*(nn-1);
        %LogEZ(nn)=log(EZ(nn-2))+log(nn-1);
        %EZ(nn);
    end
end
%EZ
%EZ(1:30)
%str=input('Look at numbers')
a(SeriesOrder+1:SeriesOrder*NMoments+1)=0;
b0=a0;
b=a;


for mm=1:NMoments
    if(mm>1)
        %[b0,b] =SeriesProductLogarithmic(a0,a,b0,b,SeriesOrder*mm);
        
        [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*mm);
        
        %b0-bb0
        
        %b-bb
        
        %str=input('Look at two products')
        
        
        b(SeriesOrder*mm+1:SeriesOrder*NMoments+1)=0;
    end
   % Logb=log(abs(b));
   % Signb=sign(b);
    EXZ(mm,1)=b0;
    for pp2=1:SeriesOrder*mm
        if rem(pp2,2)==0
            %EXZ(mm,1)=EXZ(mm,1)+exp(Logb(pp2)+LogEZ(pp2)).*Signb(pp2);
            EXZ(mm,1)=EXZ(mm,1)+b(pp2).*EZ(pp2);
 %           b(pp2).*EZ(pp2)
 %           exp(Logb(pp2)+LogEZ(pp2)).*Signb(pp2)
 %           mm
 %           pp2
            %str=input('Look at moment values')
        end
    end
end

    
for nn=1:NMoments
    Moments(nn)=EXZ(nn,1);
end



end

function [b0,b] = CalculateDensityFromStandardizedMoments(Moments0,b0Guess,bGuess)

% cmu(1)=0;
% cmu(2)=1;
% cmu(3)=Moments0(3)/Moments0(2).^1.5;
% cmu(4)=Moments0(4)/Moments0(2).^2.0;
% cmu(5)=Moments0(5)/Moments0(2).^2.50;
% cmu(6)=Moments0(6)/Moments0(2).^3.0;
% cmu(7)=Moments0(7)/Moments0(2).^3.50;
% cmu(8)=Moments0(8)/Moments0(2).^4.0;

cmu=Moments0;

%[c0,c] =PreSmoothingGuess02(cmu);
%cc=c;
%cc0=-(cc(2)+cc(4)*3+cc(6)*15);

%cc=bGuess;
%cc0=b0Guess;

%[cc0,cc] =PreSmoothingGuess02(cmu);
[cc0,cc] =PreSmoothingGuess02NewFromGuess(cmu,bGuess);

SeriesOrder=7;
NMoments=8;
NZterms=8;




[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
cmu

%str=input('Look at moments--before')



%Below calculate F, the defect vector from target cumulants. And calculate
%dF which is the matrix consisting of derivatives of defect vector from
%target cumulants.

%[F,dF] = CalculateCumulantsAndDerivativesFromMoments_0(C,cc0,cc,SeriesOrder,NMoments,8)
[F,dF] = CalculateMomentsAndDerivatives_0(cmu,cc0,cc,SeriesOrder,NZterms,NMoments);
da(1,1)=cc0;
da(2:SeriesOrder+1,1)=cc(1:SeriesOrder);



[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments);
ObjBest=0;


    
   ObjBest=100000*(abs(cmu(1)-Moments(1)))+abs(cmu(2)-Moments(2))+abs((cmu(3)-Moments(3))^(1/1.5))+abs((cmu(4)-Moments(4))^(1/2.0)) + ...
    abs((cmu(5)-Moments(5))^(1/2.0))+abs((cmu(6)-Moments(6))^(1/2.0))+abs((cmu(7)-Moments(7))^(1/2.0))+ ...
    abs((cmu(8)-Moments(8))^(1/2.0));


b0Best=cc0;
bBest(1:SeriesOrder)=cc(1:SeriesOrder);

nn=0;
while((nn<100)&&((abs(F(1,1))>.000000000001) || (abs(F(2,1))>.000000000001) || (abs(F(3,1))>.000000000001) || (abs(F(4,1))>.00000000001)|| (abs(F(5,1))>.00000000001)|| (abs(F(6,1))>.00000000001)|| (abs(F(7,1))>.00000000001)|| (abs(F(8,1))>.00000000001)  ))

nn=nn+1;
%Below Newton matrix equation to improve the Z-series coefficients guess at previous step.



da=da-dF\F;

%b0=median
b0=da(1,1);
b(1:SeriesOrder)=da(2:SeriesOrder+1,1);

%[F,dF] = CalculateCumulantsAndDerivativesFromMoments_0(C,b0,b,SeriesOrder,SeriesOrder,NoOfCumulants);
[F,dF] = CalculateMomentsAndDerivatives_0(cmu,b0,b,SeriesOrder,NZterms,NMoments);
[IsValidFlag] = CheckIsValidDensity(b0,b);

[Moments] = CalculateMomentsOfZSeries(b0,b,SeriesOrder,NMoments);
ObjNew=0;

    ObjNew=100000*(abs(cmu(1)-Moments(1)))+abs(cmu(2)-Moments(2))+abs((cmu(3)-Moments(3))^(1/1.5))+abs((cmu(4)-Moments(4))^(1/2.0)) + ...
    abs((cmu(5)-Moments(5))^(1/2.0))+abs((cmu(6)-Moments(6))^(1/2.0))+abs((cmu(7)-Moments(7))^(1/2.0))+ ...
    abs((cmu(8)-Moments(8))^(1/2.0));

  
  if((ObjBest>ObjNew) &&( IsValidFlag))
      
       ObjBest=ObjNew;
       b0Best=b0;
       bBest(1:SeriesOrder)=b(1:SeriesOrder);
   end

da(1,1)=b0;
da(2:SeriesOrder+1,1)=b(1:SeriesOrder);


end
  b0=b0Best;%Best;
  b(1:SeriesOrder)=bBest;%Best(1:SeriesOrder);
  
  %Above are the best chosen Z-series coefficients based on the objective function.
  %But my objective function can be improved.
  
[Moments] = CalculateMomentsOfZSeries(b0,b,SeriesOrder,NMoments)
%Above are raw moments calculated from the resulting Z-series from Newton
%Method

%Below are raw moments that were input to calibration.
cmu
%str=input('Look at comparison of moments calibration--00');
%b0
%b
nn
%str=input('Look at comparison of moments calibration--0');


b0=b0*sqrt(Moments0(2));
b=b*sqrt(Moments0(2));
b0;
b;
%str=input('Look at comparison of moments calibration--1');
end

function [XDensity,IndexOut,IndexMax] = MakeDensityFromSimulation_Infiniti_NEW(X,Paths,NoOfBins,MaxCutOff )
%Processes monte carlo paths to return a series Xdensity as a function of IndexOut. IndexMax is the maximum value of index.
% 


Xmin=0;
Xmax=0;
for p=1:Paths
if(X(p)>MaxCutOff)
X(p)=MaxCutOff;
end
%if(X(p)<0)
%X(p)=0;
%end
if(Xmin>real(X(p)))
Xmin=real(X(p));
end
if(Xmax<real(X(p)))
Xmax=real(X(p));
end
end

%IndexMax=NoOfBins+1;
BinSize=(Xmax-Xmin)/NoOfBins;
%IndexMax=floor((Xmax-Xmin)/BinSize+.5)+1
IndexMax=floor((Xmax-Xmin)/BinSize+.5)+1

XDensity(1:IndexMax)=0.0;

for p=1:Paths
index=real(floor(real(X(p)-Xmin)/BinSize+.5)+1);
if(real(index)<1)
index=1;
end
if(real(index)>IndexMax)
index=IndexMax;
end

XDensity(index)=XDensity(index)+1.0/Paths/BinSize;

end

IndexOut(1:IndexMax)=Xmin+(0:(IndexMax-1))*BinSize;

end
 

function [c0,c] = PreSmoothingGuess02NewFromGuess(cmu,cin)


Mul=1.0;
SeriesOrder=7;
NMoments=8;

c0=0;
c(1)=1.0;

%At this point, we have set first two coefficients of Z-series. Rest of the coefficients are zero.    
    
if(cmu(3)==0)
    c(2)=cin(2);
else
    if(cmu(3)>0)
        c(2)=.001;
    elseif(cmu(3)<0)
        c(2)=-.001;
    end
end


for nn=1:5*Mul
    [Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
    
    c(2)=c(2)*(cmu(3))/Moments(3);
    c0=-c(2);
end
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
c=c/sqrt(Moments(2));
c0=c0/sqrt(Moments(2));
%Above we have standardized the coefficients so that 2nd central moment is
%one.
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
Moments2=Moments;
%At this point first three coefficients have been set and the rest are
%zero.
%Moments2 above has been calculated with only c0, c1, c2 non-zero. All
%higher coefficients upto c7 are zero. Particularly in its calculation c3
%is zero.


%if(cmu(4)-Moments2(4)>0)
%    c(3)=.0001;
%elseif(cmu(4)-Moments2(4)<0)
%    c(3)=-.0001;
%end

c(3)=cin(3);

for nn=1:5*Mul
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%please recall that in calculation of Moments2, c3 is zero. We are taking
%slope for fourth moment from where the value of c3 is zero.
Ratio4=(cmu(4)-Moments2(4))/(Moments(4)-Moments2(4));
c(3)=c(3)*Ratio4;


if(cmu(4)<Moments2(4))
    c(3)=abs(c(3))*-1;
end
 %cmu(4)
 %Moments2(4)
 %Moments(4)
 %Ratio4
 %c(3)
% str=input('Look at fourth moment');
% 


end
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%Moments


%str=input('Look at moments--inside my function---M4');
c=c/sqrt(Moments(2));
c0=c0/sqrt(Moments(2));
%Above we have standardized the coefficients so that 2nd central moment is
%one.
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
Moments2=Moments;
%Moments2 above has been calculated with only c0, c1, c2,c3 non-zero. All
%higher coefficients upto c7 are zero. Particularly in its calculation c4
%is zero from where we will calculate slope in next loop..
%if(cmu(5)-Moments2(5)>0)
%    c(4)=.0001;
%elseif(cmu(5)-Moments2(5)<0)
%    c(4)=-.0001;
%end
    
c(4)=cin(4);

for nn=1:5*Mul
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%please recall that in calculation of Moments2, c4 is zero. We are taking
%slope for fifth moment from where the value of c4 is zero.
Ratio5=(cmu(5)-Moments2(5))/(Moments(5)-Moments2(5));
c(4)=c(4)*Ratio5;
c0=-c(2)-3*c(4);





end
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
c=c/sqrt(Moments(2));
c0=c0/sqrt(Moments(2));
%Above we have standardized the coefficients so that 2nd central moment is
%one.    
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);    
Moments2=Moments;    
%Moments2 above has been calculated with only c0, c1, c2,c3,c4 non-zero. All
%higher coefficients upto c7 are zero. Particularly in its calculation c5
%is zero from where we will calculate slope in next loop..    

%if(cmu(6)-Moments2(6)>0)
%    c(5)=.0001;
%elseif(cmu(6)-Moments2(6)<0)
%    c(5)=-.0001;
%end


c(5)=cin(5);

for nn=1:5*Mul
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%please recall that in calculation of Moments2, c5 is zero. We are taking
%slope for sixth moment from where the value of c5 is zero.
Ratio6=(cmu(6)-Moments2(6))/(Moments(6)-Moments2(6));
c(5)=c(5)*Ratio6;

end
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
c=c/sqrt(Moments(2));
c0=c0/sqrt(Moments(2));
%Above we have standardized the coefficients so that 2nd central moment is
%one.
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
Moments2=Moments;
%Moments2 above has been calculated with only c0, c1, c2,c3,c4,c5 non-zero. All
%higher coefficients upto c7 are zero. Particularly in its calculation c6
%is zero from where we will calculate slope in next loop..   

% c(6)=0;
% if(cmu(7)-Moments2(7)>0)
%     c(6)=.0001;
% 
% elseif(cmu(7)-Moments2(7)<0)
%     c(6)=-.0001;
% end


c(6)=cin(6);
    
for nn=1:5*Mul
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%please recall that in calculation of Moments2, c6 is zero. We are taking
%slope for seventh moment from where the value of c6 is zero.
Ratio7=(cmu(7)-Moments2(7))/(Moments(7)-Moments2(7));
c(6)=c(6)*Ratio7;
c0=-c(2)-3*c(4)-15*c(6);

if(cmu(6)<Moments2(6))
    c(5)=abs(c(5))*-1;
end
% cmu(6)
% Moments2(6)
% Moments(6)
% Ratio6
% c(5)
% str=input('Look at sixth moment');


end
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
c=c/sqrt(Moments(2));
c0=c0/sqrt(Moments(2));
%Above we have standardized the coefficients so that 2nd central moment is
%one.
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
Moments2=Moments;
%Moments2 above has been calculated with only c0, c1, c2,c3,c4,c5,c6 non-zero. All
%higher coefficients upto c7 are zero. Particularly in its calculation c7
%is zero from where we will calculate slope in next loop..
% c(7)=0;
% if(cmu(8)-Moments2(8)>0)
%     c(7)=.00005;
% 
% elseif(cmu(8)-Moments2(8)<0)
%     c(7)=-.00005;
% end


c(7)=cin(7);



for nn=1:5*Mul
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%please recall that in calculation of Moments2, c7 is zero. We are taking
%slope for eighth moment from where the value of c7 is zero.

Ratio8=(cmu(8)-Moments2(8))/(Moments(8)-Moments2(8));
c(7)=c(7)*(abs(Ratio8))*sign(Ratio8);
if(cmu(8)<Moments2(8))
    c(7)=abs(c(7))*-1;

% cmu(8)
% Moments2(8)
% Moments(8)
% Ratio8
% c(7)
% str=input('Look at eighth moment');
end
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
c=c/sqrt(Moments(2));
c0=c0/sqrt(Moments(2));
%Above we have standardized the coefficients so that 2nd central moment is
%one.
end

function [F,dF] = CalculateMomentsAndDerivatives_0(Ms,a0,a,SeriesOrder,NZterms,NMoments)
%[F,dF] = CalculateCumulantsAndDerivativesFromMoments(C,a0,a,SeriesOrder,SeriesOrder,NoOfCumulants);

%a(1)=a1;
%a(2)=a2;
%a(3)=a3;
%a(4)=a4;
%a(5)=a5;
%a(6)=a6;
%a(7)=a7;



%aa0=a0;
if(NMoments>8)
a0=-(a(2)+3*a(4)+15*a(6)+105*a(8));% ---1
end
if(NMoments<=8)
a0=-(a(2)+3*a(4)+15*a(6));% ---1
end

EZ(1)=0;
EZ(2)=1;
for nn=3:NMoments*SeriesOrder+NZterms+2
    if rem(nn,2)==1
        EZ(nn)=0;
    else
        EZ(nn)=EZ(nn-2)*(nn-1);
        EZ(nn);
    end
end
EZ;

EXZ(1,1)=1;
for pp1=1:NZterms
        EXZ(1,pp1+1)=EZ(pp1);
end


a(SeriesOrder+1:NMoments*SeriesOrder+1)=0;
b0=a0;
b=a;

for mm=1:NMoments
    if(mm>1)
        [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*mm);
        b(SeriesOrder*mm+1:NMoments*SeriesOrder+1)=0;
    end
   % b0
   % b
%str=input('Look at numbers')    
    EXZ(mm+1,1)=b0;
    for pp2=1:SeriesOrder*mm
        
            EXZ(mm+1,1)=EXZ(mm+1,1)+b(pp2).*EZ(pp2);
    end
    for pp1=1:NZterms
        EXZ(mm+1,pp1+1)=b0.*EZ(pp1);
        for pp2=1:SeriesOrder*mm
        
            EXZ(mm+1,pp1+1)=EXZ(mm+1,pp1+1)+b(pp2).*EZ(pp2+pp1);
        end
    end
end

%u1=EXZ(2,1);
u2=EXZ(3,1);
u3=EXZ(4,1);
u4=EXZ(5,1);

u1=a0+a(2);
if(SeriesOrder>=4)
u1=a0+a(2)+3*a(4);
end
if(SeriesOrder>=6)
u1=u1+15*a(6);
end
if(SeriesOrder>=8)
u1=u1+105*a(8);
end


%k2=u2-u1^2;
%k3=u3-3*u2*u1+2*u1^3;
%k4=u4-4*u3*u1-3*u2^2+12*u2*u1^2-6*u1^4;


du1(1)=1;%----2
% du1(2)=0;%----2
% du1(3)=1;%----2
% du1(4)=0;%----2
% du1(5)=1;%----2
% du1(6)=0;%----2
% du1(7)=15;%----2
% du1(8)=0;%----2
% du1(9)=105;%----2
% du1(10)=0;%----2
 du2(1)=2*EXZ(2,1);
 du3(1)=3*EXZ(3,1);
 du4(1)=4*EXZ(4,1);


%du1(1)=1;%----2
%du2(1)=0;
%du3(1)=0;
%du4(1)=0;

for mm=2:SeriesOrder+1
du1(mm)=EXZ(1,mm);
du2(mm)=2*EXZ(2,mm);
du3(mm)=3*EXZ(3,mm);
du4(mm)=4*EXZ(4,mm);

end


if(NMoments>=5)
u5=EXZ(6,1);


%du5(1)=5*EXZ(5,1);

for mm=1:SeriesOrder+1
du5(mm)=5*EXZ(5,mm);
end

end

if(NMoments>=6)
u6=EXZ(7,1);
%du6(1)=0;

for mm=1:SeriesOrder+1
du6(mm)=6*EXZ(6,mm);
end

end


if(NMoments>=7)
u7=EXZ(8,1);
%str=input('Look at k7 and k71')

%du7(1)=0;

for mm=1:SeriesOrder+1
du7(mm)=7*EXZ(7,mm);
end

end
if(NMoments>=8)

u8=EXZ(9,1);
for mm=1:SeriesOrder+1
du8(mm)=8*EXZ(8,mm);
end

end

if(NMoments>=9)

u9=EXZ(10,1);
for mm=1:SeriesOrder+1
du9(mm)=9*EXZ(9,mm);
end

end
if(NMoments>=10)

u10=EXZ(11,1);
for mm=1:SeriesOrder+1
du10(mm)=10*EXZ(10,mm);
end

end




 F(1,1)=u1-Ms(1);
 F(2,1)=u2-Ms(2);
 F(3,1)=u3-Ms(3);
 F(4,1)=u4-Ms(4);
 if(NMoments>=5)
 F(5,1)=u5-Ms(5);
 end
 if(NMoments>=6)
 F(6,1)=u6-Ms(6);
 end
 if(NMoments>=7)
 F(7,1)=u7-Ms(7);
 end
 if(NMoments>=8)
 F(8,1)=u8-Ms(8);
 end
 if(NMoments>=9)
 F(9,1)=u9-Ms(9);
 end
 if(NMoments>=10)
 F(10,1)=u10-Ms(10);
 end
 
     for mm=1:SeriesOrder+1
        dF(1,mm)=du1(mm);
        dF(2,mm)=du2(mm);
        dF(3,mm)=du3(mm);
        dF(4,mm)=du4(mm);
        if(NMoments>=5)
        dF(5,mm)=du5(mm);
        end
        if(NMoments>=6)
        dF(6,mm)=du6(mm);
        end
        if(NMoments>=7)
        dF(7,mm)=du7(mm);
        end
        if(NMoments>=8)
        dF(8,mm)=du8(mm);
        end
        if(NMoments>=9)
        dF(9,mm)=du9(mm);
        end
        if(NMoments>=10)
        dF(10,mm)=du10(mm);
        end
     end

end

function [EA,EB,EAB] = CalculateSeriesProducts(a0,a,b0,b,SeriesOrder,NMoments)


EA(1:NMoments)=0;
EB(1:NMoments)=0;
EAB(1:NMoments,1:NMoments)=0;
EA(1)=a0+a(2);
EB(1)=b0+b(2);
if(SeriesOrder>=4)
EA(1)=EA(1)+3*a(4);
EB(1)=EB(1)+3*b(4);

end
if(SeriesOrder>=6)
EA(1)=EA(1)+15*a(6);
EB(1)=EB(1)+15*b(6);
end


EZ(1)=0;
EZ(2)=1;
for nn=3:NMoments*SeriesOrder
    if rem(nn,2)==1
        EZ(nn)=0;
    else
        EZ(nn)=EZ(nn-2)*(nn-1);
        EZ(nn);
    end
end
EZ;



aa(1:NMoments,1:SeriesOrder*NMoments)=0;

aa0(1)=a0;
aa(1,1:SeriesOrder)=a(:);

bb(1:NMoments,1:SeriesOrder*NMoments)=0;
bb0(1)=b0;
bb(1,1:SeriesOrder)=b(:);




for mm=1:NMoments-1
        [ab0,ab] =SeriesProductHigher(a0,a,aa0(mm),aa(mm,:),SeriesOrder*(mm+1));
        aa0(mm+1)=ab0;
        aa(mm+1,1:(mm+1)*SeriesOrder)=ab(1:(mm+1)*SeriesOrder);
        
        [ac0,ac] =SeriesProductHigher(b0,b,bb0(mm),bb(mm,:),SeriesOrder*(mm+1));
        bb0(mm+1)=ac0;
        bb(mm+1,1:(mm+1)*SeriesOrder)=ac(1:(mm+1)*SeriesOrder);

        EA(mm+1)=aa0(mm+1);
        EB(mm+1)=bb0(mm+1);
        for pp=1:SeriesOrder*(mm+1)
        
            EA(mm+1)=EA(mm+1)+aa(mm+1,pp).*EZ(pp);
            EB(mm+1)=EB(mm+1)+bb(mm+1,pp).*EZ(pp);
            
        end
end
        %b(SeriesOrder*mm+1:NMoments*SeriesOrder+1)=0;
aa0;
aa;

for mm=1:NMoments
    for nn=1:NMoments-mm
        dd0=aa0(mm);
        dd=aa(mm,1:SeriesOrder*NMoments);
        ee0=bb0(nn);
        ee=bb(nn,1:SeriesOrder*NMoments);
        
        [ab0,ab] =SeriesProductHigher(dd0,dd,ee0,ee,SeriesOrder*(mm+nn));
        AB0(mm,nn)=ab0;
        AB(mm,nn,1:SeriesOrder*(mm+nn))=ab(1:SeriesOrder*(mm+nn));
        
        EAB(mm,nn)=AB0(mm,nn);
        for pp=1:SeriesOrder*(mm+nn)
        
            EAB(mm,nn)=EAB(mm,nn)+AB(mm,nn,pp).*EZ(pp);
        end
        
        
    end
        
end

end

function [c0,c] =SeriesProductHigher(a0,a,b0,b,SeriesOrder)

a(size(a,2)+1:SeriesOrder)=0;
b(size(b,2)+1:SeriesOrder)=0;
c0=a0*b0;
c(1:SeriesOrder)=0.0;
%a
%b
%str=input('Look at numbers');
for nn=1:SeriesOrder
    c(nn)=c(nn)+a0*b(nn);
    c(nn)=c(nn)+b0*a(nn);
    for kk=1:nn-1
        c(nn)=c(nn)+a(kk)*b(nn-kk);
    end
end

end

function [IsValidFlag] = CheckIsValidDensity(cc0,cc)
%Roughly checks if the density is valid. It is not a definitive check but
%just a rough one but mostly works.

Z=2.0;
Xp1=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;

Z=3.0;
Xp2=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;

Z=4.65;
Xp3=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;

Z=-2.0;
Xn1=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;

Z=-3.0;
Xn2=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;

Z=-4.65;
Xn3=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;

IsValidFlag=0;
if( ((Xp2>Xp1) && (Xp3>Xp2)) &&(Xp1>Xn1) && ((Xn2<Xn1) && (Xn3<Xn2)))
    IsValidFlag=1;
end
    

end

function [c0,c] =SeriesProduct(a0,a,b0,b,SeriesOrder)


c0=a0*b0;
c(1:SeriesOrder)=0.0;
for nn=1:SeriesOrder
    c(nn)=c(nn)+a0*b(nn);
    c(nn)=c(nn)+b0*a(nn);
    for kk=1:nn-1
        c(nn)=c(nn)+a(kk)*b(nn-kk);
    end
end

end