function [] = SVMonteCarloDensityWith2DHermitesWilmott()

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

%Volatility and asset equations are
%Vt(0)=V0;
%dVt=(mu1*Vt^beta1 + mu2*Vt^beta2)*dt+sigma0 * Vt^gamma * dz1(t)
%Xt(0)=X0
%dXt=(a*Xt^alpha1+b*Xt^alpha2)*dt+sigma1* Vt^gammaV * Xt^gammaX *dz2(t);
%<dz1(t),dz2(t)>=rho *dt

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

X0=.05;
alpha1=0;
alpha2=1;
kappaX=.850;
thetaX=0.0250;
a=kappaX*thetaX;
b=-kappaX;

rho=-0.0;
sigma1=.8;
gammaV=.5;
gammaX=.7;
%%%%%%%%%%%%%
V0=1.0;
thetaV=.250;
kappaV=1.50;
mu1=kappaV*thetaV;
mu2=-kappaV;
beta1=0;
beta2=1;
gamma=.85;
sigma0=.7;
%%%%%%%%%%%%

dt=.03125;
Tt=32*4;

T=Tt*dt;


%%%%%%%%%%%%%%%%%%%%%%%%5

%rng(29079137, 'twister')
rng(15898837, 'twister')
paths=600000;
V(1:paths)=V0;  %Original process monte carlo.
X(1:paths)=X0;
Random1(1:paths)=0;
Random2(1:paths)=0;

for tt=1:Tt
Random1=randn(size(Random1));
Random2=randn(size(Random2));


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

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


end

X=real(X);
X(X<0)=0;
V(V<0)=0;


SVolMeanAnalytic=thetaV+exp(-kappaV*T)*(V0-thetaV)
SVolMean=sum(V(1:paths))/paths
AssetMeanAnalytic=thetaX+exp(-kappaX*T)*(X0-thetaX)
AssetMean=sum(X(1:paths))/paths

str=input("Look at monte carlo Mean of X and V")

NoOfBinsV=200;
NoOfBinsX=200;

[XYProbDensity2D,XYProbMass2D,XProbMass,VProbMass,IndexMaxX,IndexMaxV,Xn,Vm,Zx,Zv]= Find2DDensityAndCDFFromSimulation_Infiniti(X,V,paths,NoOfBinsX,NoOfBinsV);
%Please look at matlab file of above function for etailed explanation of
%all its input and output parameters.

str=input("Look at sizes of variables----1")

surf(Xn,Vm,XYProbDensity2D)

%hold on

MeanX=sum(X(1:paths))/paths
MeanV=sum(V(1:paths))/paths

str=input("Look at 3D graphs")


%Below calculate Coeffs of 2D Hermite products.
%Oxn is order of Zx hermites. ranges from 1 to 8 in indexing or (0 to 7 in reality)
%Ovm is order of Zv hermites. ranges from 1 to 8 in indexing or (0 to 7 in reality)
%HermiteH(Oxn,Zx(nn,mm)) is a small function that caclulates Oxn order hermite at Zx(nn,mm) value
%HermiteH(Ovm,Zv(mm)) caclulates Ovm order hermite at Zy(mm) value
%Below Eq calculates 2D Coeffs with inner product of coordinates of 2D uniform grid with 
%product of X and V hermites.


C(1:8,1:8)=0;

for Oxn=1:8
    for Ovm=1:8
        for nn=1:(NoOfBinsX+1)
            for mm=1:(NoOfBinsV+1)
                C(Oxn,Ovm)=C(Oxn,Ovm)+1/factorial(Oxn-1)/factorial(Ovm-1).*Xn(nn).*Vm(mm).*HermiteH(Oxn-1,Zx(nn,mm)).*HermiteH(Ovm-1,Zv(mm)).*XYProbMass2D(nn,mm);
                
            end
        end
    end
end
        
C

%If C are 2D hermite Coeffs as calculated in above equation
%Hermite Series representation of 2D random variable would be given by
%double sum equation below as 
%Sum_{n=0}^{n=7} Sum_{m=0}^{m=7} C(n+1,m+1) He(n,Zx) He(m,Zy);

str=input("Look at Coefficients of 2D Hermite products")


%Below we calculate 1D Hermite Series for volatility.
%We only need to calculate the mean (zeroth hermite) aH0 of the volatility
%random variable.


aH0=0;
aH(1:7)=0;
for nn=1:IndexMaxV
aH0  =  aH0+Vm(nn).*VProbMass(nn);    
aH(1)=aH(1)+Vm(nn).*Zv(nn).*VProbMass(nn);
aH(2)=aH(2)+.5*Vm(nn).*(Zv(nn).^2-1).*VProbMass(nn);
aH(3)=aH(3)+1/6*Vm(nn).*(Zv(nn).^3-3*Zv(nn)).*VProbMass(nn);
aH(4)=aH(4)+1/24*Vm(nn).*(Zv(nn).^4-6*Zv(nn).^2+3).*VProbMass(nn);
aH(5)=aH(5)+1/120*Vm(nn).*(Zv(nn).^5-10*Zv(nn).^3+15*Zv(nn)).*VProbMass(nn);
aH(6)=aH(6)+1/720*Vm(nn).*(Zv(nn).^6-15*Zv(nn).^4+45*Zv(nn).^2-15).*VProbMass(nn);
aH(7)=aH(7)+1/720/7*Vm(nn).*(Zv(nn).^7-21*Zv(nn).^5+105*Zv(nn).^3-105*Zv(nn)).*VProbMass(nn);
end


%Below we calculate the Hermite series of asset by integrating over the
%volatility to get marginal series of asset. When we integrate over
%volatility, all volatility hermites(alone or as products with X-hermites)
%in the 2D asset series go to zero(they all have mean zero). 
%So basically all volatility terms go to zero and we are left with 
%one dimensional array containing terms from 2D series where 
%volatility hermites in the product had order zero. This order zero is
%represented by volatility index one in the 2d array of coeffs C(n,m);
%So our 1D asset series is from the line C(1:8,1) where volatilit hermites
%are zero. This is very interesting.

bH0=C(1,1);
bH(1:7)=C(2:8,1);

%Above bH0 is asset Hermite series coeff of zeroth order hermite and
%bH(1:7) represent hermites from 1st to 7th order.

bH0=bH0/aH0;
bH=bH/aH0;

%Above, we need to normalize the asset hermite series with expected value
%of volatility since we integrated over the volatility. Just like in
%Conditioning when we integrate over a variable to get marginal density and
%then we have to divide with the integral of the integrated density as
%well. Similalry, Here we have to divide with the mean/expected value of 
%volatility hermite series(which is integral of the volatility hermite series)
%to normalize our asset equation hermite coefficients.

%Below we convert from hermite series to Z-series of the asset for
%plotting.
[b0,b] = ConvertHCoeffsToZCoeffs(bH0,bH,7);

bH0
bH

b0
b

str=input("Look at 1D Coeffs of X")


dNn=.1;%.1/2;   % Normal density subdivisions width. would change with number of subdivisions
Nn=101;%45*4;  % No of normal density subdivisions
NnMid=((1+Nn)/2)*dNn;
Z(1:Nn)=(((1:Nn)*dNn)-NnMid);



XX(1:Nn)=b0;
for nn=1:7
   XX(1:Nn)=XX(1:Nn)+b(nn)*Z(1:Nn).^nn;
end

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


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



%Below make parzen density of asset in SV equations and also plot the
%analytic density of asset.
 NoOfBins=2400;
 MaxCutOff=100;
  
 [XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti_NEW(X,paths,NoOfBins,MaxCutOff );
 plot(XX(2:Nn-1),pXX(2:Nn-1),'b',IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g');

title('Asset Density(in SV model) From Analytic Hermite Series Calculated from 2D Hermite Series VS Numerical Parzen Density of Asset');%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
legend({'Marginal Asset Density from 2D Hermite Series','Parzen Density From Monte Carlo'},'Location','northeast')

str=input('blue line is density of Asset from Hermite Series, green is numerical density from monte carlo.');




end

function [XYProbDensity2D,XYProbMass2D,XProbMass,YProbMass,IndexMaxX,IndexMaxY,Xn,Ym,Zx,Zy] = Find2DDensityAndCDFFromSimulation_Infiniti(X,Y,Paths,NoOfBinsX,NoOfBinsY)

%This function processes joint observations of two stochastic variables X and Y from 
%common monte carlo paths and makes a two dimensional density of X and Y.

%Explanation of Input Parameters
%X and Y are 1D arrays so that nth element of each makes a joint 
%observation of nth monte carlo path X(n) and Y(n) form joint nth observation of two stochastic variables.
%paths ---- total monte carlo paths.
%NoOfBinsX ---- Uniform Grid size in X for density construction. You can alter
%this. But you will have to make sure to alter grid sizes in main program after
%this function has returned results.
%NoOfBinsY ---- Uniform Grid size in Y for density construction. You can alter
%this.

%Explanation of output Parameters
%XYProbDensity2D --- is 2D numerically constructed Parzen density on input
%grid in X and Y.
%XYProbMass2D, ------ is integrated version of 2D parzen density that gives
%two D probability mass on input grid in X and Y.
%XProbMass  ---- is probability mass of 1D marginal density in X on NoOfBinsX Grid.
%YProbMass  ---- is probability mass of 1D marginal density in Y on NoOfBinsY Grid.
%IndexMaxX ----- This is size of grid in X and corresponds to NoOfBinsX and
%is NoOfBinsX+1. Naturally This is also size of array that returns
%X-values on the X-axis of 2D uniform grid.
%IndexMaxY ----- This is size of grid in X and corresponds to NoOfBinsY and
%is  NoOfBinsY+1. Naturally This is also size of array that returns
%Y-values on Y-axis of the 2D uniform grid.
%Xn--- is an array of size IndexMaxX that contains values of X at the
%center of each grid cell on X-axis on 2D uniform grid.
%Ym--- is an array of size IndexMaxY that contains values of Y at the
%center of each grid cell on Y-axis on 2D uniform grid.
%Zx----- this is a 2D array. This returns the value of standard normal associated with 1D Xn grid for every
%Ym grid cell. Given a certain grid cell Ym, we calculate CDF along all 1D X grid
%cells with the same Ym. This standard normal value is derived from inversion of univariate CDF on X-grid
%or X-axis associated with each Y grid cell. Zx(n,m) is a 2D array.
%Zy----- this is the value of standard normal associated with Ym at center
%of each grid cell on Y-axis. This is derived from inversion of univariate CDF on Y-grid
%or Y-axis. This is 1D array. This is used for vol which is independent of
%X.

Xmin=999999;
Xmin1=999999;
Xmin2=999999;
Xmax=0;
Xmax1=0;
Xmax2=0;

Ymin=999999;
Ymin1=999999;
Ymin2=999999;
Ymax=0;
Ymax1=0;
Ymax2=0;
meanX=0;
meanY=0;
for p=1:Paths

meanX=meanX+X(p)/Paths;
meanY=meanY+Y(p)/Paths;
if(Xmin>real(X(p)))
    Xmin2=Xmin1;
    Xmin1=Xmin;
    Xmin=real(X(p));
end
if(Xmax<real(X(p)))
    Xmax2=Xmax1;
    Xmax1=Xmax;
    Xmax=real(X(p));
end

if(Ymin>real(Y(p)))
    Ymin2=Ymin1;
    Ymin1=Ymin;
    Ymin=real(Y(p));
end
if(Ymax<real(Y(p)))
    Ymax2=Ymax1;
    Ymax1=Ymax;
    Ymax=real(Y(p));
end
end


%Xmin
%Xmin1
%Xmin2

%Xmax
%Xmax1
%Xmax2
%str=input('Look at Xmin and Xmax');
%IndexMax=NoOfBins+1;
BinSizeX=(Xmax2-Xmin2)/NoOfBinsX;
BinSizeY=(Ymax2-Ymin2)/NoOfBinsY;
%IndexMax=floor((Xmax-Xmin)/BinSize+.5)+1
IndexMaxX=floor((Xmax2-Xmin2)/BinSizeX+.5)+1
IndexMaxY=floor((Ymax2-Ymin2)/BinSizeY+.5)+1

Xn(1:IndexMaxX)=Xmin2+(0:(IndexMaxX-1))*BinSizeX;
Ym(1:IndexMaxY)=Ymin2+(0:(IndexMaxY-1))*BinSizeY;

XYProbMass2D(1:IndexMaxX,1:IndexMaxY)=0.0;
XYProbDensity2D(1:IndexMaxX,1:IndexMaxY)=0.0;
XProbMass(1:IndexMaxX)=0;
YProbMass(1:IndexMaxY)=0;

for p=1:Paths
indexX=real(floor(real(X(p)-Xmin2)/BinSizeX+.5)+1);
%indexXp=real(X(p)-Xmin2)/BinSizeX+1;


if(real(indexX)<1)
indexX=1;
end
if(real(indexX)>IndexMaxX)
indexX=IndexMaxX;
end


indexY=real(floor(real(Y(p)-Ymin2)/BinSizeY+.5)+1);
%indexYp=real(Y(p)-Ymin2)/BinSizeY;

if(real(indexY)<1)
indexY=1;
end
if(real(indexY)>IndexMaxY)
indexY=IndexMaxY;
end



%XDensity(index)=XDensity(index)+1.0/Paths/BinSize;
XProbMass(indexX)=XProbMass(indexX)+1.0/Paths;
YProbMass(indexY)=YProbMass(indexY)+1.0/Paths;
XYProbMass2D(indexX,indexY)=XYProbMass2D(indexX,indexY)+1.0/Paths;
XYProbDensity2D(indexX,indexY)=XYProbDensity2D(indexX,indexY)+1.0/Paths/BinSizeX/BinSizeY;


end

%BinSizeX
%BinSizeY

%str=input("Look at BinSizes");

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%Below calculate Zx by inversion of 1D X-CDF associated with each volatility grid cell
%so that resulting array is 2D.
XProbSum(1:IndexMaxX,1:IndexMaxY)=0;
XCDF(1:IndexMaxX,1:IndexMaxY)=0;
Zx(1:IndexMaxX,1:IndexMaxY)=0;

for mm=1:IndexMaxY
    XProbSum(1,mm)=XYProbMass2D(1,mm)*.5;
    if(YProbMass(mm)>0)
        XCDF(1,mm)=XProbSum(1,mm)./YProbMass(mm);
        
    else
        XCDF(1,mm)=0;
    end
    if((XCDF(1,mm)>0) && (XCDF(1,mm)<1))
        Zx(1,mm)=norminv(XCDF(1,mm));
    elseif (XCDF(1,mm)<.0000001)
        Zx(1,mm)=-8;
    elseif (XCDF(1,mm)>.9999999)
        Zx(1,mm)=8;
    end
    for nn=2:IndexMaxX
        XProbSum(nn,mm)=XProbSum(nn-1,mm)+XYProbMass2D(nn-1,mm)*.5+XYProbMass2D(nn,mm)*.5;
        if(YProbMass(mm)>0)
            XCDF(nn,mm)=XProbSum(nn,mm)./YProbMass(mm);
        else
            XCDF(nn,mm)=0;
        end
        if((XCDF(nn,mm)>0) && (XCDF(nn,mm)<1))
            Zx(nn,mm)=norminv(XCDF(nn,mm));
        elseif (XCDF(nn,mm)<.0000001)
            Zx(nn,mm)=-8;
        elseif (XCDF(nn,mm)>.9999999)
            Zx(nn,mm)=8;
        end
    end
end
    
    

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



YCDF(1:IndexMaxY)=0;
Zy(1:IndexMaxY)=0;
YCDF(1)=YProbMass(1)*.5;
Zy(1)=norminv(YCDF(1));
for nn=2:IndexMaxY
    YCDF(nn)=YCDF(nn-1)+YProbMass(nn-1)*.5+YProbMass(nn)*.5;
    Zy(nn)=norminv(YCDF(nn));
end

end

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


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

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==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==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

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


end

function [ValHermite] = HermiteH(n,Z)

if(n==0)
    ValHermite=1;
elseif(n==1)
    ValHermite=Z;
elseif(n==2)
    ValHermite=Z.^2-1;
elseif(n==3)
    ValHermite=Z.^3-3*Z;
elseif(n==4)
    ValHermite=Z.^4-6*Z.^2+3;
elseif(n==5)
    ValHermite=Z.^5-10*Z.^3+15*Z;
elseif(n==6)
    ValHermite=Z.^6-15*Z.^4+45*Z.^2-15; 
elseif(n==7)
    ValHermite=Z.^7-21*Z.^5+105*Z.^3-105*Z;     
end

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=10000;
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