Friends here are the matlab functions required for Z-series density calibration to first eight raw moments.

.

.

```
function [c0,c,FittedMoments] = CalculateZSeriesDensityFromRawMoments02(rMu)
%The above function takes an array of first eight raw moments as input.
%And returns Z-series Coeffs upto 7th power of Z as in
%Y=c0+c(1)*Z+c(2)*Z^2+c(3)*Z^3+c(4)*Z^4+c(5)*Z^5+c(6)*Z^6+c(7)*Z^7
%Fitted Moments are raw moments of the Z-series density and usually are exactly
%the same as input raw moments rMu(1:8).i.e FittedMoments(1:8)==rMu(1:8)
%with high precision.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%If you do not have valid moments from a density available, please
%uncomment the following input moments to run the program.
rMu(1)=0.000012499863483*10^5;
rMu(2)=0.000021197093700*10^5;
rMu(3)=0.000050274347857*10^5;
rMu(4)=0.000171368317209*10^5;
rMu(5)=0.000856094970826*10^5;
rMu(6)=0.006320436615909*10^5;
rMu(7)=0.068698344472325*10^5;
rMu(8)=1.084082326087001*10^5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
SeriesOrder=7; %Z-series till 7th power. hardcoded here cannot chagne.
mOrder=8; %Number of raw moments in calibration. hardcoded here cannot chagne.
NMoments=mOrder;
%Find central Moments from raw moments
[Mu1,cMu] = ConvertRawMomentsToCentralMoments(rMu,mOrder);
%Find standardized Moments from central moments
sMu(1)=0;
sMu(2)=1;
sMu(3)=cMu(3)/cMu(2).^1.5;
sMu(4)=cMu(4)/cMu(2).^2.0;
sMu(5)=cMu(5)/cMu(2).^2.5;
sMu(6)=cMu(6)/cMu(2).^3.0;
sMu(7)=cMu(7)/cMu(2).^3.5;
sMu(8)=cMu(8)/cMu(2).^4.0;
iter=5; %Number of iterations in the calibration.
%Please experiment and play around.
bGuess(2:7)=0; %We have assigned initial values of coefficients
bGuess(1)=1; %other than c(1) to zero. We are starting with no prior information.
[c0,c] = PreSmoothingGuessAdvancedFromGuessBestNew(sMu,bGuess,iter);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Comment this block if you do not want multidimensional Newton pass
%Even though preSmoothing guess is very close to optimal, slight
%inaccuracies can change the shape of probabioity density so this
%Multidimensional Newton is usually required.
[F,dF] = CalculateMomentsAndDerivatives_0(sMu,c0,c,SeriesOrder,SeriesOrder,NMoments);
da(1,1)=c0;
da(2:SeriesOrder+1,1)=c(1:SeriesOrder);
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%Replace with your own more intelligent objective function if you like.
ObjBest=100000*(abs(sMu(1)-Moments(1)))+abs(sMu(2)-Moments(2))+abs((sMu(3)-Moments(3))^(1/1.5))+abs((sMu(4)-Moments(4))^(1/2.0)) + ...
abs((sMu(5)-Moments(5))^(1/2.0))+abs((sMu(6)-Moments(6))^(1/2.0))+abs((sMu(7)-Moments(7))^(1/2.0))+ ...
abs((sMu(8)-Moments(8))^(1/2.0));
b0Best=c0;
bBest(1:SeriesOrder)=c(1:SeriesOrder);
nn=0;
while((nn<20)&&((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(sMu,b0,b,SeriesOrder,SeriesOrder,NMoments);
[IsValidFlag] = CheckIsValidDensity(b0,b);
[Moments] = CalculateMomentsOfZSeries(b0,b,SeriesOrder,NMoments);
ObjNew=100000*(abs(sMu(1)-Moments(1)))+abs(sMu(2)-Moments(2))+abs((sMu(3)-Moments(3))^(1/1.5))+abs((sMu(4)-Moments(4))^(1/2.0)) + ...
abs((sMu(5)-Moments(5))^(1/2.0))+abs((sMu(6)-Moments(6))^(1/2.0))+abs((sMu(7)-Moments(7))^(1/2.0))+ ...
abs((sMu(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
c0=b0Best;%Best;
c(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.
%Multidimensional Newton root search ends. Comment till here if you like to
%exclude this
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Make sure you do not comment the lines below
c0=c0*sqrt(cMu(2)); %convert from standardized to central moments
c=c*sqrt(cMu(2)); %convert from standardized to central moments
c0=c0+Mu1; %convert from central to raw moments by adding mean.
[FittedMoments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Please Uncomment this if you want to see a comparison of input moments and
%fitted moments
rMu
FittedMoments
str=input('Look at comparison of fitted moments with input moments');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%COMMENT lines below if you do not want to see the density of new Z-series random variable with
%fitted moments.
%First calculate normal random variable on a grid below
dNn=.1/2; % Normal density subdivisions width. would change with number of subdivisions
Nn=45*4; % No of normal density subdivisions
NnMid=((1+Nn)/2)*dNn;
Z(1:Nn)=(((1:Nn)*dNn)-NnMid);
%Z
%str=input('Look at Z');
%Now calculate Y as a function of normal random variable.
Y(1:Nn)=c0;
for nn=1:7
Y(1:Nn)=Y(1:Nn)+c(nn)*Z(1:Nn).^nn;
end
%Now take change of densities derivative of Y with respect to normal
DfY(1:Nn)=0;
for nn=2:Nn-1
DfY(nn) = (Y(nn + 1) - Y(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
DfY(Nn)=DfY(Nn-1);
DfY(1)=DfY(2);
%Now calculate the density of Y from density of normal random variable
%using change of probability derivative.
pY(1:Nn)=0;
for nn = 1:Nn
pY(nn) = (normpdf(Z(nn),0, 1))/abs(DfY(nn));
end
plot(Y(1:Nn),pY(1:Nn),'b');
title('Density of the Z-Series Variable with Fitted Moments');
%str=input('Look at the density of Z-series random variable')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
```

.

.

.

```
function [c0,c] = PreSmoothingGuessAdvancedFromGuessBestNew(cmu,cin,iter)
%Mul=1.0;
SeriesOrder=7;
NMoments=8;
EZ(1)=0;
EZ(2)=1;
for nn=3:NMoments*SeriesOrder+SeriesOrder+2
if rem(nn,2)==1
EZ(nn)=0;
else
EZ(nn)=EZ(nn-2)*(nn-1);
EZ(nn);
end
end
c0=-(cin(2)+3*cin(4)+15*cin(6)); %This is for zero mean condition.
c=cin;
MaxIter=10;
for nn=1:iter %increase the iterations over coefficients if neeeded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
tolerance=.0000001*cmu(3);
mm=0;
while ((abs(M3-cmu(3)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(2)=c(2)-(M3-cmu(3))/dc2;
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,3,2,EZ);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
tolerance=.0000001*cmu(4);
mm=0;
while ((abs(M4-cmu(4)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(3)=c(3)-(M4-cmu(4))/dc3;
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,4,2,EZ);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
tolerance=.0000001*cmu(5);
mm=0;
while ((abs(M5-cmu(5)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(4)=c(4)-(M5-cmu(5))/dc4;
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,5,2,EZ);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
tolerance=.0000001*cmu(6);
mm=0;
while ((abs(M6-cmu(6)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(5)=c(5)-(M6-cmu(6))/dc5;
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,6,2,EZ);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M7,dc6] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,7,EZ);
tolerance=.0000001*cmu(7);
mm=0;
while ((abs(M7-cmu(7)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(6)=c(6)-(M7-cmu(7))/dc6;
[M7,dc6] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,7,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,7,2,EZ);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M8,dc7] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,8,EZ);
tolerance=.0000001*cmu(8);
mm=0;
while ((abs(M8-cmu(8)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(7)=c(7)-(M8-cmu(8))/dc7;
[M8,dc7] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,8,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,8,2,EZ);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
```

.

.

.

```
function [EnMoment,dEnMomentdCoeff] = CalculateParticularMomentAndDerivativeOfItsCoeff(a0,a,SeriesOrder,nMoment,EZ)
% if(SeriesOrder>=8)
% a0=-(a(2)+3*a(4)+15*a(6)+105*a(8));% ---1
% end
% if(SeriesOrder<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:nMoment*SeriesOrder+1+SeriesOrder)=0;
b0=a0;
b=a;
for mm=2:nMoment
[b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*mm+SeriesOrder+1);
b(SeriesOrder*mm+1:nMoment*SeriesOrder+SeriesOrder+1)=0;
if(mm==nMoment-1)
dEnMomentdCoeff=b0.*EZ(nMoment-1);
for pp2=1:SeriesOrder*nMoment-1
dEnMomentdCoeff=dEnMomentdCoeff+b(pp2).*EZ(pp2+nMoment-1);
end
dEnMomentdCoeff=dEnMomentdCoeff*nMoment;
end
if(mm==nMoment)
EnMoment=b0;
for pp2=1:SeriesOrder*nMoment
EnMoment=EnMoment+b(pp2).*EZ(pp2);
end
end
end
end
```

.

.

.

```
function [c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,nMoment,iter,EZ)
MaxIter=7;
SeriesOrder=7;
c0=-(c(2)+3*c(4)+15*c(6));
if(nMoment==3)
for nn=1:iter %increase the iterations over coefficients if neeeded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
tolerance=.0000001*cmu(3);
mm=0;
while ((abs(M3-cmu(3)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(2)=c(2)-(M3-cmu(3))/dc2;
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
if(nMoment==4)
for nn=1:iter %increase the iterations over coefficients if neeeded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
tolerance=.0000001*cmu(3);
mm=0;
while ((abs(M3-cmu(3)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(2)=c(2)-(M3-cmu(3))/dc2;
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
tolerance=.0000001*cmu(4);
mm=0;
while ((abs(M4-cmu(4)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(3)=c(3)-(M4-cmu(4))/dc3;
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(nMoment==5)
for nn=1:iter %increase the iterations over coefficients if neeeded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
tolerance=.0000001*cmu(3);
mm=0;
while ((abs(M3-cmu(3)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(2)=c(2)-(M3-cmu(3))/dc2;
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
tolerance=.0000001*cmu(4);
mm=0;
while ((abs(M4-cmu(4)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(3)=c(3)-(M4-cmu(4))/dc3;
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
tolerance=.0000001*cmu(5);
mm=0;
while ((abs(M5-cmu(5)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(4)=c(4)-(M5-cmu(5))/dc4;
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(nMoment==6)
for nn=1:iter %increase the iterations over coefficients if neeeded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
tolerance=.0000001*cmu(5);
mm=0;
while ((abs(M3-cmu(3)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(2)=c(2)-(M3-cmu(3))/dc2;
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
tolerance=.0000001*cmu(5);
mm=0;
while ((abs(M4-cmu(4)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(3)=c(3)-(M4-cmu(4))/dc3;
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
tolerance=.0000001*cmu(5);
mm=0;
while ((abs(M5-cmu(5)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(4)=c(4)-(M5-cmu(5))/dc4;
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
tolerance=.0000001*cmu(6);
mm=0;
while ((abs(M6-cmu(6)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(5)=c(5)-(M6-cmu(6))/dc5;
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(nMoment==7)
for nn=1:iter %increase the iterations over coefficients if neeeded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
tolerance=.0000001*cmu(3);
mm=0;
while ((abs(M3-cmu(3)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(2)=c(2)-(M3-cmu(3))/dc2;
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
tolerance=.0000001*cmu(4);
mm=0;
while ((abs(M4-cmu(4)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(3)=c(3)-(M4-cmu(4))/dc3;
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
tolerance=.0000001*cmu(5);
mm=0;
while ((abs(M5-cmu(5)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(4)=c(4)-(M5-cmu(5))/dc4;
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
tolerance=.0000001*cmu(6);
mm=0;
while ((abs(M6-cmu(6)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(5)=c(5)-(M6-cmu(6))/dc5;
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M7,dc6] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,7,EZ);
tolerance=.0000001*cmu(7);
mm=0;
while ((abs(M7-cmu(7)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(6)=c(6)-(M7-cmu(7))/dc6;
[M7,dc6] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,7,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(nMoment==8)
for nn=1:iter %increase the iterations over coefficients if neeeded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
tolerance=.0000001*cmu(3);
mm=0;
while ((abs(M3-cmu(3)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(2)=c(2)-(M3-cmu(3))/dc2;
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
tolerance=.0000001*cmu(4);
mm=0;
while ((abs(M4-cmu(4)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(3)=c(3)-(M4-cmu(4))/dc3;
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
tolerance=.0000001*cmu(5);
mm=0;
while ((abs(M5-cmu(5)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(4)=c(4)-(M5-cmu(5))/dc4;
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
tolerance=.0000001*cmu(6);
mm=0;
while ((abs(M6-cmu(6)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(5)=c(5)-(M6-cmu(6))/dc5;
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M7,dc6] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,7,EZ);
tolerance=.0000001*cmu(7);
mm=0;
while ((abs(M7-cmu(7)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(6)=c(6)-(M7-cmu(7))/dc6;
[M7,dc6] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,7,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M8,dc7] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,8,EZ);
tolerance=.0000001*cmu(8);
mm=0;
while ((abs(M8-cmu(8)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(7)=c(7)-(M8-cmu(8))/dc7;
[M8,dc7] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,8,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
end
```

.

.

.

```
function [SecondMoment] = CalculateSecondMomentC7(a0,a)
SecondMoment=a0^2+a(1)^2 +2* a0* a(2) +3*(a(2).^2+2* a(1).* a(3) +2* a0* a(4))+ ...
15*(a(3).^2 +2 *a(2).* a(4) +2 *a(1).* a(5)+2* a0* a(6))+105*(a(4).^2 +2* a(3).* a(5)+2* a(2).* a(6) +2* a(1).* a(7) )+ ...
945*(a(5).^2 +2* a(4).* a(6) +2 *a(3).* a(7) )+(10395)*(a(6).^2+2* a(5).* a(7) )+(135135)*a(7).^2;
end
```

.

.

.

```
function [mu1,cMu] = ConvertRawMomentsToCentralMoments(Mu,mOrder)
%cu
%mu
% %mOrder
%str=input('Look at numbers');
mu1=Mu(1);
for mm=2:mOrder
cMu(mm)=0;
for jj=0:mm
if(jj==0)
cMu(mm)=cMu(mm)+ (-1)^(mm-jj).*factorial(mm)/factorial(jj)/factorial(mm-jj)*1*mu1.^(mm-jj);
else
cMu(mm)=cMu(mm)+ (-1)^(mm-jj).*factorial(mm)/factorial(jj)/factorial(mm-jj)*Mu(jj)*mu1.^(mm-jj);
end
end
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 [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
```

.

.

.

```
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 [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
```

.

.

.

When you run the program, you will see a comparison of input moments (rMu) and fitted moments.

rMu =

1.0e+05 *

Columns 1 through 7

0.000012499863483 0.000021197093700 0.000050274347857 0.000171368317209 0.000856094970826 0.006320436615909 0.068698344472325

Column 8

1.084082326087001

FittedMoments =

1.0e+05 *

Columns 1 through 7

0.000012499863483 0.000021197093700 0.000050274347857 0.000171368317209 0.000856094970826 0.006320436615909 0.068698344472325

Column 8

1.084082326087000

Look at comparison of fitted moments with input moments.

And you will see the following graph of the density.

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