Here is also a quick Case 6 from a few days' earlier program. I have also specified its cumulants under heading case 6 in the program.

.

.

Here are the matlab program files.

.

.

```
function [b0,b] = FindZSeriesGivenCumulantsStandard08A(k1,k2,k3,k4,k5,k6,k7,k8)
mu=0;
sigma=1.0;
%Case1
% k1=0.0;
% k2=(1.0);
% k3=-.45;%1.0;
% k4=1.750;%-.50;
% k5=-3;%-.5;%15;
% k6=8;%10;%100.0;
% k7=-10;
% k8=200;
%Case 2
% k1=0.0;
% k2=(1.0);
% k3=-.40;%1.0;
% k4=1.5;%-.50;
% k5=-3.80;%-.5;%15;
% k6=11.0;%10;%100.0;
% k7=-13.0;
% k8=1500;
%Case 3
% k1=0.0;
% k2=(1.0);
% k3=-.5;%1.0;
% k4=1.150;%-.50;
% k5=-2;%-.5;%15;
% k6=12;%10;%100.0;
% k7=-11;
% k8=1000;
%Case 4
% k1=0.0;
% k2=(1.0);
% k3=-.9;%1.0;
% k4=1.950;%-.50;
% k5=-5;%-.5;%15;
% k6=19.5;%10;%100.0;
% k7=-23;
% k8=3850;
%Case 5
% k1=0.0;
% k2=(1.0);
% k3=-.45;%1.0;
% k4=.850;%-.50;
% k5=-2.15;%-.5;%15;
% k6=5.5;%10;%100.0;
% k7=-10;
% k8=1500;
%Case 6
k1=0.0;
k2=(1.0);
k3=-1.3;%1.0;
k4=2.50;%-.50;
k5=-3.0;%-.5;%15;
k6=8.0;%10;%100.0;
k7=-40;
k8=2500;
[mu,sigma,a0,a1,a2,a3,a4,a5,mu1,rmu1,rmu2,rmu3,rmu4,rmu5,rmu6] = GenerateDensityFromCumulantsX06Simple(mu,sigma,k1,k2,k3,k4,k5,k6)
rmu7=k1*rmu6+6*k2*rmu5+15*k3*rmu4+20*k4*rmu3+15*k5*rmu2+6*k6*rmu1+k7;
rmu8=k1*rmu7+7*k2*rmu6+21*k3*rmu5+35*k4*rmu4+35*k5*rmu3+21*k6*rmu2+7*k7*rmu1+k8;
rmu(1)=rmu1;
rmu(2)=rmu2;
rmu(3)=rmu3;
rmu(4)=rmu4;
rmu(5)=rmu5;
rmu(6)=rmu6;
rmu(7)=rmu7;
rmu(8)=rmu8;
[median] = FindMedianHermiteDensity06(mu1,mu,sigma,a0,a1,a2,a3,a4,a5);
Xm=median;
%str=input('Look at medians');
%Below calculate the density and its derivatives at median.
px=exp(-0.5* ((Xm-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((Xm-mu)/sigma)+ a2* (((Xm-mu)/sigma).^2-1) + ...
a3* ( ((Xm-mu)/sigma).^3 - 3*((Xm-mu)/sigma) )+ ...
a4* (((Xm-mu)/sigma).^4 - 6* ((Xm-mu)/sigma).^2 + 3)+ ...
a5* (((Xm-mu)/sigma).^5 - 10* ((Xm-mu)/sigma).^3 + 15*((Xm-mu)/sigma)));
dpxdX=exp(-0.5* ((Xm-mu)/sigma).^2)/sigma^2/sqrt(2*pi).* -(a0* ((Xm-mu)/sigma)+ a1* (((Xm-mu)/sigma).^2-1) + ...
a2* ( ((Xm-mu)/sigma).^3 - 3*((Xm-mu)/sigma) )+ ...
a3* (((Xm-mu)/sigma).^4 - 6* ((Xm-mu)/sigma).^2 + 3)+ ...
a4* (((Xm-mu)/sigma).^5 - 10* ((Xm-mu)/sigma).^3 + 15*((Xm-mu)/sigma))+ ...
a5* (((Xm-mu)/sigma).^6 - 15* ((Xm-mu)/sigma).^4 + 45*((Xm-mu)/sigma).^2-15));
d2pxdX2=exp(-0.5* ((Xm-mu)/sigma).^2)/sigma^3/sqrt(2*pi).* (a0* (((Xm-mu)/sigma).^2-1) + ...
a1* ( ((Xm-mu)/sigma).^3 - 3*((Xm-mu)/sigma) )+ ...
a2* (((Xm-mu)/sigma).^4 - 6* ((Xm-mu)/sigma).^2 + 3)+ ...
a3* (((Xm-mu)/sigma).^5 - 10* ((Xm-mu)/sigma).^3 + 15*((Xm-mu)/sigma))+ ...
a4* (((Xm-mu)/sigma).^6 - 15* ((Xm-mu)/sigma).^4 + 45*((Xm-mu)/sigma).^2-15) + ...
a5* (((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma)));
d3pxdX3=exp(-0.5* ((Xm-mu)/sigma).^2)/sigma^4/sqrt(2*pi).* -(a0* (((Xm-mu)/sigma).^3 - 3*((Xm-mu)/sigma) )+ ...
a1* (((Xm-mu)/sigma).^4 - 6* ((Xm-mu)/sigma).^2 + 3)+ ...
a2* (((Xm-mu)/sigma).^5 - 10* ((Xm-mu)/sigma).^3 + 15*((Xm-mu)/sigma))+ ...
a3* (((Xm-mu)/sigma).^6 - 15* ((Xm-mu)/sigma).^4 + 45*((Xm-mu)/sigma).^2-15) + ...
a4* (((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma))+ ...
a5* (((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105));
d4pxdX4=exp(-0.5* ((Xm-mu)/sigma).^2)/sigma^5/sqrt(2*pi).* ( a0* (((Xm-mu)/sigma).^4 - 6* ((Xm-mu)/sigma).^2 + 3)+ ...
a1* (((Xm-mu)/sigma).^5 - 10* ((Xm-mu)/sigma).^3 + 15*((Xm-mu)/sigma))+ ...
a2* (((Xm-mu)/sigma).^6 - 15* ((Xm-mu)/sigma).^4 + 45*((Xm-mu)/sigma).^2-15) + ...
a3* (((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma))+ ...
a4* (((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105)+ ...
a5* (((Xm-mu)/sigma).^9 - 36* ((Xm-mu)/sigma).^7 + 378*((Xm-mu)/sigma).^5-1260*((Xm-mu)/sigma).^3+945*((Xm-mu)/sigma)));
d5pxdX5=exp(-0.5* ((Xm-mu)/sigma).^2)/sigma^6/sqrt(2*pi).* -( a0* (((Xm-mu)/sigma).^5 - 10* ((Xm-mu)/sigma).^3 + 15*((Xm-mu)/sigma))+ ...
a1* (((Xm-mu)/sigma).^6 - 15* ((Xm-mu)/sigma).^4 + 45*((Xm-mu)/sigma).^2-15) + ...
a2* (((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma))+ ...
a3* (((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105)+ ...
a4* (((Xm-mu)/sigma).^9 - 36* ((Xm-mu)/sigma).^7 + 378*((Xm-mu)/sigma).^5-1260*((Xm-mu)/sigma).^3+945*((Xm-mu)/sigma))+ ...
a5* (((Xm-mu)/sigma).^10 - 45* ((Xm-mu)/sigma).^8 + 630*((Xm-mu)/sigma).^6-3150*((Xm-mu)/sigma).^4+4725*((Xm-mu)/sigma).^2- 945 ));
d6pxdX6=exp(-0.5* ((Xm-mu)/sigma).^2)/sigma^7/sqrt(2*pi).* (a0* (((Xm-mu)/sigma).^6 - 15* ((Xm-mu)/sigma).^4 + 45*((Xm-mu)/sigma).^2-15) + ...
a1* (((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma))+ ...
a2* (((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105)+ ...
a3* (((Xm-mu)/sigma).^9 - 36* ((Xm-mu)/sigma).^7 + 378*((Xm-mu)/sigma).^5-1260*((Xm-mu)/sigma).^3+945*((Xm-mu)/sigma))+ ...
a4* (((Xm-mu)/sigma).^10 - 45* ((Xm-mu)/sigma).^8 + 630*((Xm-mu)/sigma).^6-3150*((Xm-mu)/sigma).^4+4725*((Xm-mu)/sigma).^2- 945 ) + ...
a5* (((Xm-mu)/sigma).^11 - 55* ((Xm-mu)/sigma).^9 + 990*((Xm-mu)/sigma).^7-6930*((Xm-mu)/sigma).^5+17325*((Xm-mu)/sigma).^3-10395*((Xm-mu)/sigma) ));
%Below calculate standard normal density and its derivatives at median.
Zm=0;
pz=exp(-0.5* (Zm.^2))/sqrt(2*pi);
dpzdZ=-Zm*exp(-0.5* (Zm.^2))/sqrt(2*pi);
d2pzdZ2=(Zm^2-1)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
d3pzdZ3=-(Zm^3-3*Zm)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
d4pzdZ4=(Zm^4-6*Zm^2+3)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
d5pzdZ5=-(Zm^5-10*Zm^3+15*Zm)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
d6pzdZ6=(Zm^6-15*Zm^4+45*Zm.^2-15)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
dXdZ=pz/px;
d2XdZ2=(dpzdZ-(dpxdX *dXdZ^2))/px
d3XdZ3=(d2pzdZ2- (2* dpxdX* dXdZ *d2XdZ2+dXdZ *(dXdZ^2* d2pxdX2+dpxdX *d2XdZ2)))/px
d4XdZ4=(d3pzdZ3-(3* d2XdZ2* (dXdZ^2* d2pxdX2+dpxdX *d2XdZ2)+3* dpxdX *dXdZ *d3XdZ3+dXdZ* (3 *dXdZ* d2pxdX2* d2XdZ2+dXdZ^3* d3pxdX3+dpxdX* d3XdZ3)))/px
d5XdZ5=(d4pzdZ4- (6 *(dXdZ^2 *d2pxdX2+dpxdX *d2XdZ2) *d3XdZ3+4 *d2XdZ2* (3* dXdZ *d2pxdX2 *d2XdZ2+dXdZ^3* d3pxdX3+dpxdX* d3XdZ3)+4* dpxdX* dXdZ* d4XdZ4+dXdZ *(3* d2pxdX2* d2XdZ2^2+6* dXdZ^2* d2XdZ2* d3pxdX3+4 *dXdZ* d2pxdX2* d3XdZ3+dXdZ^4* d4pxdX4+dpxdX* d4XdZ4)))/px
d6XdZ6=(d5pzdZ5- (10* d3XdZ3 *(3 *dXdZ* d2pxdX2* d2XdZ2+dXdZ^3* d3pxdX3+dpxdX *d3XdZ3)+10* (dXdZ^2* d2pxdX2+dpxdX *d2XdZ2)* d4XdZ4+5* d2XdZ2* (3* d2pxdX2* d2XdZ2^2+6* dXdZ^2* d2XdZ2* d3pxdX3+4 *dXdZ* d2pxdX2* d3XdZ3+dXdZ^4* d4pxdX4+dpxdX* d4XdZ4)+5 *dpxdX* dXdZ* d5XdZ5+dXdZ *(15* dXdZ* d2XdZ2^2* d3pxdX3+10* d2pxdX2* d2XdZ2* d3XdZ3+10* dXdZ^2 *d3pxdX3* d3XdZ3+10* dXdZ^3* d2XdZ2* d4pxdX4+5* dXdZ* d2pxdX2* d4XdZ4+dXdZ^5* d5pxdX5+dpxdX* d5XdZ5)))/px
d7XdZ7=(d6pzdZ6- (20 *(3* dXdZ* d2pxdX2* d2XdZ2+dXdZ^3 *d3pxdX3+dpxdX *d3XdZ3)* d4XdZ4+15* d3XdZ3* (3 *d2pxdX2 *d2XdZ2^2+6* dXdZ^2* d2XdZ2* d3pxdX3+4 *dXdZ* d2pxdX2* d3XdZ3+dXdZ^4* d4pxdX4+dpxdX* d4XdZ4)+15* (dXdZ^2* d2pxdX2+dpxdX* d2XdZ2)* d5XdZ5+6* d2XdZ2* (15* dXdZ* d2XdZ2^2* d3pxdX3+10* d2pxdX2* d2XdZ2 *d3XdZ3+10 *dXdZ^2 *d3pxdX3 *d3XdZ3+10 *dXdZ^3 *d2XdZ2* d4pxdX4+5 *dXdZ *d2pxdX2 *d4XdZ4+dXdZ^5 *d5pxdX5+dpxdX* d5XdZ5)+6* dpxdX* dXdZ *d6XdZ6+dXdZ* (15 *d2XdZ2^3 *d3pxdX3+60 *dXdZ *d2XdZ2 *d3pxdX3* d3XdZ3+10 *d2pxdX2 *d3XdZ3^2+45 *dXdZ^2* d2XdZ2^2 *d4pxdX4+20* dXdZ^3 *d3XdZ3* d4pxdX4+15 *d2pxdX2* d2XdZ2 *d4XdZ4+15 *dXdZ^2* d3pxdX3 *d4XdZ4+15 *dXdZ^4* d2XdZ2* d5pxdX5+6 *dXdZ *d2pxdX2* d5XdZ5+dXdZ^6* d6pxdX6+dpxdX *d6XdZ6)))/px
%Below calculate iniotial guess coefficients of Z-series to be input to
%Newton method.
c0=median;
c1=dXdZ;
c2=1/2*1*d2XdZ2;
c3=1/6*1*d3XdZ3;
c4=1/24*1*d4XdZ4;
c5=1/120*1*d5XdZ5;
c6=1/720*1*d6XdZ6;
c7=1/720/7*d7XdZ7;
%earlier program worked only because odd coefficients that are related to
%even moments were reasonablyclose to true values. Even coefficients
%responsible for odd moments were usually grossly off. Below I try to make
%sure that sign of even coefficients matches the sign of odd moments. This
%is not completely definitive but usually right.
cc0=c0;
cc(1)=c1;
if(sign(k3)~=sign(c2))
cc(2)=-1*c2;
else
cc(2)=c2;
end
cc(3)=c3;
if(sign(k5)~=sign(c4))
cc(4)=-1*c4;
else
cc(4)=c4;
end
cc(5)=c5;
if(sign(k7)~=sign(c6))
cc(6)=-1*c6;
else
cc(6)=c6;
end
cc(7)=c7;
cc0=-(cc(2)+cc(4)*3+cc(6)*15);
SeriesOrder=7;
NMoments=8;
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
rmu
str=input('Look at moments--before')
%Below is the function that finds a guess for Newton method using heuristic
%bisection.
[cc0,cc] = ImproveNewtonGuess02(cc0,cc,rmu,25);
%str=input('Look at initial match of moments for the Newton guess');
%Below specify target cumulants.
C(1:NMoments+1)=0;
C(1)=k1;
C(2)=k2;
C(3)=k3;
C(4)=k4;
C(5)=k5;
C(6)=k6;
C(7)=k7;
C(8)=k8;
%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)
da(1,1)=cc0;
da(2:SeriesOrder+1,1)=cc(1:SeriesOrder);
NoOfCumulants=NMoments
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments);
ObjBest=0;
%ObjBest=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
% abs((rmu(5)-Moments(5))^(1/2.5))+abs((rmu(6)-Moments(6))^(1/3.0))+abs((rmu(7)-Moments(7))^(1/3.5))+ ...
% abs((rmu(8)-Moments(8))^(1/4.0));
ObjBest=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
abs((rmu(5)-Moments(5))^(1/2.0))+abs((rmu(6)-Moments(6))^(1/2.0))+abs((rmu(7)-Moments(7))^(1/2.0))+ ...
abs((rmu(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) ))
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);
[IsValidFlag] = CheckIsValidDensity(b0,b)
[Moments] = CalculateMomentsOfZSeries(b0,b,SeriesOrder,NMoments);
ObjNew=0;
%ObjNew=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
% abs((rmu(5)-Moments(5))^(1/2.5))+abs((rmu(6)-Moments(6))^(1/3.0))+abs((rmu(7)-Moments(7))^(1/3.5))+ ...
% abs((rmu(8)-Moments(8))^(1/4.0));
ObjNew=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
abs((rmu(5)-Moments(5))^(1/2.0))+abs((rmu(6)-Moments(6))^(1/2.0))+abs((rmu(7)-Moments(7))^(1/2.0))+ ...
abs((rmu(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.
rmu
str=input('Look at comparison of moments calibration');
%Below we do calculations to construct the density from our newly derived
%Z-series.
dNn=.2/2; % Normal density subdivisions width. would change with number of subdivisions
Nn=90; % No of normal density subdivisions
NnMid=4.0;
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);
w(1:Nn)=b0;
for nn=1:SeriesOrder
w(1:Nn)=w(1:Nn)+b(nn)*Z(1:Nn).^nn;
end
wnStart=1;
Dfw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
Dfw(nn) = (w(nn + 1) - w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
end
pw(1:Nn)=0;
for nn = wnStart+1:Nn-1
pw(nn) = (normpdf(Z(nn),0, 1))/abs(Dfw(nn));
end
%
% %Below is the construction of density of our originial analytic density for
% %plotting purposes
%
dMm=.02/10;
Mm=1250*5;
X(1:Mm)=-5+dMm+dMm*(1:Mm);
%
% %Below is the newly constructed density so that its output central moments
% %match the central moments input by the user.
% pdf=exp(-0.5* ((X-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X-mu)/sigma)+ a2* (((X-mu)/sigma).^2-1) + ...
% a3* ( ((X-mu)/sigma).^3 - 3*((X-mu)/sigma) )+ ...
% a4* (((X-mu)/sigma).^4 - 6* ((X-mu)/sigma).^2 + 3)+ ...
% a5* (((X-mu)/sigma).^5 - 10* ((X-mu)/sigma).^3 + 15*((X-mu)/sigma))+ ...
% a6* (((X-mu)/sigma).^6 - 15* ((X-mu)/sigma).^4 + 45*((X-mu)/sigma).^2-15))+ ...
% a7*(((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma));
%
% %Below is standard Gaussian PDF.
pdfG=exp(-0.5* ((X-mu)/sigma).^2)/sigma/sqrt(2*pi);
%
plot(w(wnStart+1:Nn-1),pw(wnStart+1:Nn-1),'r',X(1:Mm),pdfG(1:Mm),'b');
title(sprintf('Fm1 = %.4f,m1=%.4f,Fm2=%.2f,m2=%.2f,Fm3=%.2f,m3 =%.2f,Fm4=%.2f,m4=%.2f,Fm5=%.2f,m5=%.2f,Fm6=%.1f,m6=%.1f,Fm7=%.1f,m7=%.1f,Fm8=%.1f,m8=%.1f',Moments(1),rmu1,Moments(2),rmu2,Moments(3),rmu3,Moments(4),rmu4,Moments(5),rmu5,Moments(6),rmu6,Moments(7),rmu7,Moments(8),rmu8));
%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
legend({'Newton Density','Analytical Density','Standard Gaussian Density'},'Location','northeast')
cc0
b0
cc
b
end
```

.

.

.

.

```
function [cc0Best,ccBest] = ImproveNewtonGuess02(cc0,cc,rmu,Iterations)
%This function finds a good guess for Newton method by heuristic bisection.
SeriesOrder=7;
NMoments=8;
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments);
%ObjectiveBest=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
% abs((rmu(5)-Moments(5))^(1/2.5))+abs((rmu(6)-Moments(6))^(1/3.0))+abs((rmu(7)-Moments(7))^(1/3.5))+ ...
% abs((rmu(8)-Moments(8))^(1/4.0));
ObjectiveBest=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
abs((rmu(5)-Moments(5))^(1/2.0))+abs((rmu(6)-Moments(6))^(1/2.0))+abs((rmu(7)-Moments(7))^(1/2.0))+ ...
abs((rmu(8)-Moments(8))^(1/2.0));
ccBest=cc;
cc0Best=cc0;
r1=.5;%bisection weight one
r2=.5;%bisection weight two
for nn=1:Iterations
if(nn>1)
%r1=r1+.025;
%r2=r2-.025;
end
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
Ratio2=(r1*rmu(2)+r2*Moments(2))/Moments(2);
cc(1)=cc(1)*sqrt(Ratio2);
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
Ratio3=(r1*rmu(3)+r2*Moments(3))/Moments(3);
cc(2)=cc(2)*(Ratio3);
cc0=-(cc(2)+3*cc(4)+15*cc(6));
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
Ratio4=(r1*rmu(4)+r2*Moments(4))/Moments(4);
cc(3)=cc(3)*(Ratio4);
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
Ratio5=(r1*rmu(5)+r2*Moments(5))/Moments(5);
cc(4)=cc(4)*(Ratio5);
cc0=-(cc(2)+3*cc(4)+15*cc(6));
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
Ratio6=(r1*rmu(6)+r2*Moments(6))/Moments(6);
cc(5)=cc(5)*(Ratio6);
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
Ratio7=(r1*rmu(7)+r2*Moments(7))/Moments(7);
cc(6)=cc(6)*(Ratio7);
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
Ratio8=(r1*rmu(8)+r2*Moments(8))/Moments(8);
cc(7)=cc(7)*(Ratio8);
cc0=-(cc(2)+3*cc(4)+15*cc(6));
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
[IsValidFlag] = CheckIsValidDensity(cc0,cc)
%ObjectiveNew=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
% abs((rmu(5)-Moments(5))^(1/2.5))+abs((rmu(6)-Moments(6))^(1/3.0))+abs((rmu(7)-Moments(7))^(1/3.5))+ ...
% abs((rmu(8)-Moments(8))^(1/4.0));
ObjectiveNew=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
abs((rmu(5)-Moments(5))^(1/2.0))+abs((rmu(6)-Moments(6))^(1/2.0))+abs((rmu(7)-Moments(7))^(1/2.0))+ ...
abs((rmu(8)-Moments(8))^(1/2.0));
if((ObjectiveNew<ObjectiveBest) &&( IsValidFlag))
ObjectiveBest=ObjectiveNew;
ccBest=cc;
cc0Best=cc0;
end
end
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
rmu
str=input('Look at moments--In function---After');
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 [mu,sigma,a0,a1,a2,a3,a4,a5,mu1,rmu1,rmu2,rmu3,rmu4,rmu5,rmu6] = GenerateDensityFromCumulantsX06Simple(mu,sigma,k1,k2,k3,k4,k5,k6)
%In this program, you specify the first six cumulants through input parametersand the program creates a density based on these
%cumulantsents.
%The algorithm takes a base gaussian density that has it mean mu and
%sq-root of variance given by sigma.
%This mean and sigma can be different from first and sq-root of second
%central moment(variance) of the target density that we specified as mu1 and mu2 earlier.
%Below is the construction of grid to calculate pdf, numerically integrate the pdf
%to find cdf. And to numerically calculate the six moments on the constructed density
%so they can be matched with input moments.
dNn=.02/10;
Nn=1250*10;
X(1:Nn)=-10+dNn+dNn*(1:Nn);
%Below, we calculate first six raw-moments from the input cumulants.
%These raw moments would be used to calculate the coefficients of hermite
%polynomials so that output central moments calculated on new density match
%the input central moments.
mu1=k1;
rmu1=k1;
rmu2=k1*rmu1+k2;
rmu3=k1*rmu2+2*k2*rmu1+k3;
rmu4=k1*rmu3+3*k2*rmu2+3*k3*rmu1+k4;
rmu5=k1*rmu4+4*k2*rmu3+6*k3*rmu2+4*k4*rmu1+k5;
rmu6=k1*rmu5+5*k2*rmu4+10*k3*rmu3+10*k4*rmu2+5*k5*rmu1+k6;
%rmu7=k1*rmu6+6*k2*rmu5+15*k3*rmu4+20*k4*rmu3+15*k5*rmu2+6*k6*rmu1+k7;
%rmu8=k1*rmu7+7*k2*rmu6+21*k3*rmu5+35*k4*rmu4+35*k5*rmu3+21*k6*rmu2+7*k7*rmu1+k8;
%Below is base gaussian pdf.
pdfG=exp(-0.5* ((X-mu)/sigma).^2)/sigma/sqrt(2*pi);
%Below I calculate the coefficients of polynomials making the density so
%that the resulting density has the same as specified raw moments.
%This follows my description of method 4 in the post 1491 given at
%https://forum.wilmott.com/viewtopic.php?f=4&t=99702&p=871093#p871047
MOrder=6+6
SeriesOrder=6;
SeriesOrder=6;
[EZ] = CalculateGaussianMoments(MOrder)
%str=input('Look at Z Expectations')
CMatrix(1,1:SeriesOrder)=EZ(1:SeriesOrder);
CMatrix(2,1:SeriesOrder)=EZ(2:SeriesOrder+1);
CMatrix(3,1:SeriesOrder)=EZ(3:SeriesOrder+2);
CMatrix(4,1:SeriesOrder)=EZ(4:SeriesOrder+3);
CMatrix(5,1:SeriesOrder)=EZ(5:SeriesOrder+4);
CMatrix(6,1:SeriesOrder)=EZ(6:SeriesOrder+5);
%CMatrix(7,1:SeriesOrder)=EZ(7:SeriesOrder+6);
%CMatrix(8,1:SeriesOrder)=EZ(8:SeriesOrder+7);
Mu(1,1)=1.0;
Mu(2,1)=rmu1;
Mu(3,1)=rmu2;
Mu(4,1)=rmu3;
Mu(5,1)=rmu4;
Mu(6,1)=rmu5;
%Mu(7,1)=rmu6;
%Mu(8,1)=rmu7;
Coeff=CMatrix\Mu;
Coeff
%str=input('Look at Mu results');
%Below I convert from simple polynomials to hermite polynomials so that
%resulting density remains identical with the initial density.
[CoeffaH] = ConvertZCoeffsToHCoeffsAB(Coeff,SeriesOrder)
a0=CoeffaH(1);
a1=CoeffaH(2);
a2=CoeffaH(3);
a3=CoeffaH(4);
a4=CoeffaH(5);
a5=CoeffaH(6);
a6=0;%CoeffaH(7);
a7=0;%CoeffaH(8);
a8=0;
CoeffaH(7)=0;
CoeffaH(8)=0;
%str=input('Look at Mu results---2');
%Below is the newly constructed density so that its output raw moments
%match the raw moments input by the user.
%Sometimes this density is negative at places but we still use it to derive
%an initial guess for construction of Z-series Newton density.
pdf=exp(-0.5* ((X-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X-mu)/sigma)+ a2* (((X-mu)/sigma).^2-1) + ...
a3* ( ((X-mu)/sigma).^3 - 3*((X-mu)/sigma) )+ ...
a4* (((X-mu)/sigma).^4 - 6* ((X-mu)/sigma).^2 + 3)+ ...
a5* (((X-mu)/sigma).^5 - 10* ((X-mu)/sigma).^3 + 15*((X-mu)/sigma))+ ...
a6* (((X-mu)/sigma).^6 - 15* ((X-mu)/sigma).^4 + 45*((X-mu)/sigma).^2-15)+ ...
a7* (((X-mu)/sigma).^7 - 21* ((X-mu)/sigma).^5 + 105*((X-mu)/sigma).^3-105*((X-mu)/sigma))+ ...
a8* (((X-mu)/sigma).^8 - 28* ((X-mu)/sigma).^6 + 210*((X-mu)/sigma).^4-420*((X-mu)/sigma).^2+105));
%Below cdf is calculated by doing numerical integration of the density. I
%have replaced integration with summations.
cdf0(1:Nn)=0;
for nn=2:Nn
cdf0(nn)=cdf0(nn-1)+exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2/2* (((X(nn)-mu)/sigma).^2-1) + ...
a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15) + ...
a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
end
%Below is numerically calculated value of first raw moment. N in mu1N
%stands for numerical.
rmu1N=0;
for nn=2:Nn
rmu1N=rmu1N+X(nn) .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2/2* (((X(nn)-mu)/sigma).^2-1) + ...
a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15)+ ...
a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
end
%Below is numerically calculated value of second raw moment.
rmu2N=0;
for nn=2:Nn
rmu2N=rmu2N+(X(nn)-0*rmu1N).^2 .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2* (((X(nn)-mu)/sigma).^2-1) + ...
a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15)+ ...
a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
end
%Below is numerically calculated value of third raw moment.
rmu3N=0;
for nn=2:Nn
rmu3N=rmu3N+(X(nn)-0*rmu1N).^3 .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2* (((X(nn)-mu)/sigma).^2-1) + ...
a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15)+ ...
a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
end
%Below is numerically calculated value of fourth raw moment.
rmu4N=0;
for nn=2:Nn
rmu4N=rmu4N+(X(nn)-0*rmu1N).^4 .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2* (((X(nn)-mu)/sigma).^2-1) + ...
a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15)+ ...
a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
end
%Below is numerically calculated value of fifth raw moment.
rmu5N=0;
for nn=2:Nn
rmu5N=rmu5N+(X(nn)-0*rmu1N).^5 .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2* (((X(nn)-mu)/sigma).^2-1) + ...
a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15) + ...
a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
end
% %Below is numerically calculated value of sixth raw moment.
% rmu6N=0;
% for nn=2:Nn
% rmu6N=rmu6N+(X(nn)-0*rmu1N).^6 .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2* (((X(nn)-mu)/sigma).^2-1) + ...
% a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
% a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
% a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
% a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15) + ...
% a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
% a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
% end
%
% rmu7N=0;
% for nn=2:Nn
% rmu7N=rmu7N+(X(nn)-0*rmu1N).^7 .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2* (((X(nn)-mu)/sigma).^2-1) + ...
% a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
% a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
% a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
% a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15) + ...
% a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
% a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
% end
%
% rmu8N=0;
% for nn=2:Nn
% rmu8N=rmu8N+(X(nn)-0*rmu1N).^8 .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2* (((X(nn)-mu)/sigma).^2-1) + ...
% a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
% a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
% a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
% a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15) + ...
% a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
% a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
% end
%Below is analytically calculated cdf.
cdf(1:Nn)=normcdf(((X(1:Nn)-mu)/sigma))*a0+exp(-0.5* ((X(1:Nn)-mu)/sigma).^2)/sqrt(2*pi).* (- a1 - a2* ((X(1:Nn)-mu)/sigma) - ...
a3* ( ((X(1:Nn)-mu)/sigma).^2 - 1 )- ...
a4* (((X(1:Nn)-mu)/sigma).^3 - 3* ((X(1:Nn)-mu)/sigma)) - ...
a5* (((X(1:Nn)-mu)/sigma).^4 - 6* ((X(1:Nn)-mu)/sigma).^2 + 3) - ...
a6* (((X(1:Nn)-mu)/sigma).^5 - 10* ((X(1:Nn)-mu)/sigma).^3 + 15*((X(1:Nn)-mu)/sigma)));
pdf02=exp(-0.5* ((X-mu)/sigma).^2)/sigma/sqrt(2*pi).* (CoeffaH(1)+ CoeffaH(2)* ((X-mu)/sigma)+ CoeffaH(3)* (((X-mu)/sigma).^2-1) + ...
CoeffaH(4)* ( ((X-mu)/sigma).^3 - 3*((X-mu)/sigma) )+ ...
CoeffaH(5)* (((X-mu)/sigma).^4 - 6* ((X-mu)/sigma).^2 + 3)+ ...
CoeffaH(6)* (((X-mu)/sigma).^5 - 10* ((X-mu)/sigma).^3 + 15*((X-mu)/sigma))+ ...
CoeffaH(7)* (((X-mu)/sigma).^6 - 15* ((X-mu)/sigma).^4 + 45*((X-mu)/sigma).^2-15)+ ...
CoeffaH(8)* (((X-mu)/sigma).^7 - 21* ((X-mu)/sigma).^5 + 105*((X-mu)/sigma).^3-105*((X-mu)/sigma)));%+ ...
%a8/720/56* (((X-mu)/sigma).^8 - 28* ((X-mu)/sigma).^6 + 210*((X-mu)/sigma).^4-420*((X-mu)/sigma).^2+105));
plot(X(1:Nn),cdf0(1:Nn),'g',X(1:Nn),cdf(1:Nn),'r',X(1:Nn),pdf(1:Nn),'b',X(1:Nn),pdfG(1:Nn),'m',X(1:Nn),pdf02(1:Nn),'k');
legend({'Numerically Calculated CDF','Analytically Calculated CDF','Newly Constructed Moment Matched PDF','Base Gaussian PDF'},'Location','northeast')
%Below program checks if density becomes slightly negative.
IsNegative=0;
for nn=2:Nn
if(pdf(nn)<0)
IsNegative=1;
end
end
%cdf
%Below compare the input raw moments rmu's with moments calculated on
%newly constructed density denoted by rmu#N
rmu1N
rmu1
rmu2N
rmu2
rmu3N
rmu3
rmu4N
rmu4
rmu5N
rmu5
%rmu6N
%rmu6
%rmu7N
%rmu7
%rmu8N
%rmu8
%str=input('Look at match of moments');
%Below are coefficients on hermite polynomials found out so that output
%moments in newly constructed density match input central moment values.
a1
a2
a3
a4
a5
%a6
%a7
%a8
IsNegative
%str=input('Look if the density is negative');
end
```

.

.

.

```
function [median] = FindMedianHermiteDensity06(mu1,mu,sigma,a0,a1,a2,a3,a4,a5)
%Uses Newton-Raphson to calculate the median.
X=mu1;
cdf=normcdf(((X-mu)/sigma))*a0+exp(-0.5* ((X-mu)/sigma).^2)/sqrt(2*pi).* (- a1 - a2* ((X-mu)/sigma) - ...
a3* ( ((X-mu)/sigma).^2 - 1 )- ...
a4* (((X-mu)/sigma).^3 - 3* ((X-mu)/sigma)) - ...
a5* (((X-mu)/sigma).^4 - 6* ((X-mu)/sigma).^2 + 3));% - ...
%a6/720* (((X-mu)/sigma).^5 - 10* ((X-mu)/sigma).^3 + 15*((X-mu)/sigma)));
nn=0;
while (abs(cdf-.5)> .00000000001)
nn=nn+1;
pdf=exp(-0.5* ((X-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X-mu)/sigma)+ a2* (((X-mu)/sigma).^2-1) + ...
a3* ( ((X-mu)/sigma).^3 - 3*((X-mu)/sigma) )+ ...
a4* (((X-mu)/sigma).^4 - 6* ((X-mu)/sigma).^2 + 3)+ ...
a5* (((X-mu)/sigma).^5 - 10* ((X-mu)/sigma).^3 + 15*((X-mu)/sigma)));%+ ...
%a6/720* (((X-mu)/sigma).^6 - 15* ((X-mu)/sigma).^4 + 45*((X-mu)/sigma).^2-15));
X=X-(cdf-.5)/pdf;
cdf=normcdf(((X-mu)/sigma))*a0+exp(-0.5* ((X-mu)/sigma).^2)/sqrt(2*pi).* (- a1 - a2* ((X-mu)/sigma) - ...
a3* ( ((X-mu)/sigma).^2 - 1 )- ...
a4* (((X-mu)/sigma).^3 - 3* ((X-mu)/sigma)) - ...
a5* (((X-mu)/sigma).^4 - 6* ((X-mu)/sigma).^2 + 3));% - ...
%a6/720* (((X-mu)/sigma).^5 - 10* ((X-mu)/sigma).^3 + 15*((X-mu)/sigma)));
cdf
end
nn
cdf
mu1
median=X
% dNn=.002/100;
% Nn=1250*1000;
% X(1:Nn)=-10+dNn+dNn*(1:Nn);
%
% Mm=round((median+10-dNn)/dNn);
%
%
%
% cdf0(1:Nn)=0;
% for nn=2:Mm
% cdf0(nn)=cdf0(nn-1)+exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (1+ a1* ((X(nn)-mu)/sigma)+ a2/2* (((X(nn)-mu)/sigma).^2-1) + ...
% a3/6* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
% a4/24* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
% a5/120* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
% a6/720* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15))*dNn;
% end
%
%
% cdf0(Mm)
%
% str=input('Look at median cdf');
end
```

.

.

.

```
function [F,dF] = CalculateCumulantsAndDerivativesFromMoments_0(C,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;
a0=0;% ---1
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);
k1=aa0+a(2);
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;
% du1(1)=1;%----2
% du2(1)=2*EXZ(2,1);
% du3(1)=3*EXZ(3,1);
% du4(1)=4*EXZ(4,1);
du1(1)=0;%----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);
k5=u5-5*u4*u1-10*u3*u2+20*u3*u1^2+30*u2^2*u1-60*u2*u1^3+24*u1^5;
du5(1)=0;
for mm=2:SeriesOrder+1
du5(mm)=5*EXZ(5,mm);
end
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;
du6(1)=0;
for mm=2:SeriesOrder+1
du6(mm)=6*EXZ(6,mm);
end
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
k7= u7- 7 *u1* u6 + 42* u1^2* u5 - 21* u2 *u5- 210* u1^3 *u4+ ...
210* u1*u2* u4 - 35* u3 *u4 + 840 *u1^4* u3 + 210* u2^2* u3 + ...
210*12* u1^3* u2^2 - 1260* u1^2 * u2* u3 - 630 *u1* u2^3 + 140*u1* u3^2 - 2520* u1^5* u2 + 720* u1^7 ;
%str=input('Look at k7 and k71')
du7(1)=0;
for mm=2:SeriesOrder+1
du7(mm)=7*EXZ(7,mm);
end
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;
du8(1)=0;
for mm=2:SeriesOrder+1
du8(mm)=8*EXZ(8,mm);
end
end
% 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);
% du5(mm)=5*EXZ(5,mm);
% du6(mm)=6*EXZ(6,mm);
% du7(mm)=7*EXZ(7,mm);
% du8(mm)=8*EXZ(8,mm);
% end
dk1(1)=1.0;
dk1(2)=0.0;
dk1(3)=1.0;
dk1(4)=0.0;
if(SeriesOrder>=4)
dk1(5)=3.0;
end
if(SeriesOrder>=5)
dk1(6)=0.0;
end
if(SeriesOrder>=6)
dk1(7)=15.0;
end
if(SeriesOrder>=7)
dk1(8)=0.0;
end
if(SeriesOrder>=8)
dk1(9)=105.0;
end
% dk2(1)=0;
% dk3(1)=0;
% dk4(1)=0;
% dk5(1)=0;
% dk6(1)=0;
% dk7(1)=0;
% dk8(1)=0;
for mm=2:SeriesOrder+1
%dk1(mm)=du1(mm);
dk2(mm)=-2.* u1.*du1(mm)+du2(mm);
dk3(mm)=du3(mm)-3*du1(mm)*u2-3*u1*du2(mm)+2*3*u1^2*du1(mm);
dk4(mm)=du4(mm)-4*du3(mm)*u1-4*u3*du1(mm)-3*2*u2*du2(mm)+12*du2(mm)*u1^2+12*2*u2*u1*du1(mm)-6*4*u1^3*du1(mm);
if(NMoments>=5)
dk5(mm)=du5(mm)-5*du4(mm)*u1-5*u4*du1(mm)-10*du3(mm)*u2-10*u3*du2(mm)+20*du3(mm)*u1^2+20*2*u3*u1*du1(mm)+ ...
30*2*u2*du2(mm)*u1+30*u2^2*du1(mm)-60*du2(mm)*u1^3-60*3*u2*u1^2*du1(mm)+24*5*u1^4*du1(mm);
end
if(NMoments>=6)
dk6(mm)=du6(mm)-6*du5(mm)*u1-6*u5*du1(mm)-15*du4(mm)*u2-15*u4*du2(mm)+30*du4(mm)*u1^2+30*2*u4*u1*du1(mm)-10*2*u3*du3(mm)+ ...
120*du3(mm)*u2*u1-120*u3*du2(mm)*u1-120*u3*u2*du1(mm)-120*du3(mm)*u1^3-120*3*u3*u1^2*du1(mm)+30*3*u2^2*du2(mm) - ...
270*2*u2*du2(mm)*u1^2+270*2*u2^2*u1*du1(mm)+360*du2(mm)*u1^4+360*4*u2*u1^3*du1(mm)-120*6*u1^5*du1(mm);
end
if(NMoments>=7)
dk7(mm)= du7(mm)- 7 *du1(mm)* u6- 7 *u1* du6(mm) + 42*2* du1(mm)* u5+ 42* u1^2* du5(mm) ...
- 21* du2(mm) *u5- 21* u2 *du5(mm)- 210* 3*u1^2*du1(mm) *u4- 210* u1^3 *du4(mm) + ...
210* du1(mm)*u2* u4+210* u1*du2(mm)* u4+210* u1*u2* du4(mm) - 35* du3(mm) *u4- 35* u3 *du4(mm) ...
+ 840 *4*u1^3*du1(mm)* u3 + 840 *u1^4* du3(mm)+ 210* 2*u2*du2(mm)* u3 + 210* u2^2* du3(mm)+ ...
210*12* 3*u1^2*du1(mm)* u2^2 +210*12* u1^3* 2*u2*du2(mm) - 1260* 2*u1*du1(mm) * u2* u3- 1260* u1^2 * du2(mm)* u3 ...
- 1260* u1^2 * u2* du3(mm) - 630 *du1(mm)* u2^3- 630 *u1* 3*u2^2*du2(mm) ...
+140*du1(mm)* u3^2+140*u1* 2*u3*du3(mm) - 2520* 5*u1^4*du1(mm)* u2- 2520* u1^5* du2(mm) + 720* 7*u1^6*du1(mm);
%dk7(mm)= du7(mm)- 7 *du1(mm)* u6 + 42* u1^2* u5 - 21* u2 *u5- 210* u1^3 *u4+ ...
% 210* u1*u2* u4 - 35* u3 *u4 + 840 *u1^4* u3 + 210* u2^2* u3 + ...
% 210*12* u1^3* u2^2 - 1260* u1^2 * u2* u3 - 630 *u1* u2^3 + 140*u1* u3^2 - 2520* u1^5* u2 + 720* u1^7
%dk7(mm)= du7(mm) - du1(mm)* u6- u1* du6(mm) - 6* dk2(mm)* u5 - 6 *k2 *du5(mm) - 15* dk3(mm)* u4- 15* k3* du4(mm) ...
% - 20* dk4(mm)* u3- 20* k4* du3(mm) -15* dk5(mm)* u2-15* k5* du2(mm) - 6* dk6(mm)* u1- 6 *k6* du1(mm);
%str=input('Look at comparison of derivatives of 7th cumulant');
end
if(NMoments>=8)
dk8(mm)= du8(mm) - du1(mm)* u7- u1* du7(mm) - 7* dk2(mm)* u6- 7* k2* du6(mm) - 21* dk3(mm)* u5- 21* k3* du5(mm) ...
-35* dk4(mm)* u4-35* k4* du4(mm) - 35* dk5(mm)* u3- 35* k5* du3(mm) - 21* dk6(mm)* u2- 21* k6 *du2(mm) ...
- 7* dk7(mm)* u1- 7 *k7* du1(mm);
end
end
F(1,1)=k1-C(1);
F(2,1)=k2-C(2);
F(3,1)=k3-C(3);
F(4,1)=k4-C(4);
if(NMoments>=5)
F(5,1)=k5-C(5);
end
if(NMoments>=6)
F(6,1)=k6-C(6);
end
if(NMoments>=7)
F(7,1)=k7-C(7);
end
if(NMoments>=8)
F(8,1)=k8-C(8);
end
for mm=1:SeriesOrder+1
dF(1,mm)=dk1(mm);
dF(2,mm)=dk2(mm);
dF(3,mm)=dk3(mm);
dF(4,mm)=dk4(mm);
if(NMoments>=5)
dF(5,mm)=dk5(mm);
end
if(NMoments>=6)
dF(6,mm)=dk6(mm);
end
if(NMoments>=7)
dF(7,mm)=dk7(mm);
end
if(NMoments>=8)
dF(8,mm)=dk8(mm);
end
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
```

.

.

.

You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal