Breakthrough in the theory of stochastic differential equations and their simulation

Amin · July 25th, 2022, 7:36 am

With standardized moments, i.e k1=0, and k2=1, you can only specify upto k3=+/-.4, if k4=0, k5=0, k6=0. Similarly you can only specify k4=+1.5 when other cumulants are zero. You can increase k4 to 4 or 5 but for that you would have to judiciously increase k6. Similarly you can increase k3 to +/-2 but for that you need to have enough of k4, k5, and k6.You can slightly decrease K4 to -.5 while keeping k3,k5 equal to zero but you have to specify a small positive value of k6=2.5(for k4=-.5).

For significant skew without significant kurtosis, either you have to transform this density just like we did from Bessel coordinates to original CEV coordinates. Its Z-series can be transformed to Z-series of the target transformed density. Or you would have to take a gamma density(that inhenrently takes large skews) as base density and then try to compare its first few derivatives to gaussian density to find the derivatves dnXdZn that would be used as a guess for Newton method in order to directly find the Z-series of gamma type densities.

Here are the matlab program files used to create above densities.

function [] = FindZSeriesGivenCumulantsPolynomialSimple06()


mu=0;
sigma=1.0;
k1=0.0;
k2=(1.0);
k3=-0.0;%1.0;
k4=-1;%-.50;
k5=0;%-.5;%15;
k6=10;%10;%100.0;



[mu,sigma,a0,a1,a2,a3,a4,a5,mu1,rmu1,rmu2,rmu3,rmu4,rmu5,rmu6] = GenerateDensityFromCumulantsX06Simple(mu,sigma,k1,k2,k3,k4,k5,k6)

[median] = FindMedianHermiteDensity06(mu1,mu,sigma,a0,a1,a2,a3,a4,a5);
a6=0;
a7=0;
a8=0;


%[median] = FindMedianHermiteDensity08Simple(mu1,mu,sigma,a0,a1,a2,a3,a4,a5,a6,a7,a8)

Xm=median

%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))+ ...
    a6* (((Xm-mu)/sigma).^6 - 15* ((Xm-mu)/sigma).^4 + 45*((Xm-mu)/sigma).^2-15)+ ...
    a7*(((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma))+ ...
    a8*(((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105));

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) + ...
    a6* (((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma))+...
    a7*(((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105)+ ...
    a8*(((Xm-mu)/sigma).^9 - 36* ((Xm-mu)/sigma).^7 + 378*((Xm-mu)/sigma).^5-1260*((Xm-mu)/sigma).^3+945*((Xm-mu)/sigma)));

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))+ ...
    a6* (((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105) + ...
    a7*(((Xm-mu)/sigma).^9 - 36* ((Xm-mu)/sigma).^7 + 378*((Xm-mu)/sigma).^5-1260*((Xm-mu)/sigma).^3+945*((Xm-mu)/sigma))+ ...
    a8/720/56*(((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));


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)+ ...
    a6* (((Xm-mu)/sigma).^9 - 36* ((Xm-mu)/sigma).^7 + 378*((Xm-mu)/sigma).^5-1260*((Xm-mu)/sigma).^3+945*((Xm-mu)/sigma)) + ...
    a7*(((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)+ ...
    a8*(((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)));


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))+ ...
    a6* (((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 ) + ...
    a7*(((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))+ ...
    a8*(((Xm-mu)/sigma).^12 - 66* ((Xm-mu)/sigma).^10 + 1485*((Xm-mu)/sigma).^8-13860*((Xm-mu)/sigma).^6+51975*((Xm-mu)/sigma).^4-62370.*((Xm-mu)/sigma).^2+10395));




dpxdX
d2pxdX2
d3pxdX3
d4pxdX4

%str=input('Look at derivatives results')

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


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 

%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

cc0=c0;
cc(1)=c1;
cc(2)=c2;
cc(3)=c3;
cc(4)=c4;
cc(5)=c5;

SeriesOrder=5;
NMoments=6;


[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
rmu1
rmu2
rmu3
rmu4
rmu5
rmu6


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



%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,6)

da(1,1)=cc0;
da(2:SeriesOrder+1,1)=cc(1:SeriesOrder);


NoOfCumulants=NMoments

ObjBest=0;
 for mm=1:NoOfCumulants
     if(C(mm)~=0)
     ObjBest=ObjBest+F(mm,1).^2/(abs(C(mm)));
     end
 end

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

  ObjNew=0;
  for mm=1:NoOfCumulants
      if(C(mm)~=0)
      ObjBest=ObjBest+F(mm,1).^2/(abs(C(mm)));
      end
  end
  
  if(ObjBest<ObjNew)
       
   else
       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.
  
  
  
  [F,dF] = CalculateCumulantsAndDerivativesFromMoments_0(C,b0,b,SeriesOrder,SeriesOrder,NoOfCumulants);
  %[F,dF,uu] = CalculateCumulantsAndDerivativesFromMoments_0A(C,b0,b,SeriesOrder,SeriesOrder,NoOfCumulants);
 ObjLast=0;
 for mm=1:NoOfCumulants
     if(C(mm)~=0)
     ObjLast=ObjLast+F(mm,1).^2/(abs(C(mm)));
     end
 end
ObjLast
  
nn

F


[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.
rmu1
rmu2
rmu3
rmu4
rmu5
rmu6
%str=input('Look at new moments');


%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),pdf(1:Mm),'g',X(1:Mm),pdfG(1:Mm),'b');

title(sprintf('Fm1 = %.2f,m1=%.2f,Fm2=%.2f,m2=%.2f,Fm3=%.2f,m3 =%.2f,Fm4=%.2f,m4=%.2f,Fm5=%.2f,m5=%.2f,Fm6=%.2f,m6=%.3f,b0=%.4f,b1=%.4f,b2=%.4f,b3=%.4f,b4=%.4f,b5=%.4f',Moments(1),rmu1,Moments(2),rmu2,Moments(3),rmu3,Moments(4),rmu4,Moments(5),rmu5,Moments(6),rmu6,b0,b(1),b(2),b(3),b(4),b(5)));
%,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')


b0
b




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;
[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 [EZ] = CalculateGaussianMoments(MOrder)
%Because matlab does not start indexing from zero, The indexing has been
%moved ahead by one number.

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

.
.
The function below finds an equivalent hermite representation of an algebraic expression in powers of X. Indexing is moved ahead by one point because matlab does not start indexing at zero.
.

function [aH] = ConvertZCoeffsToHCoeffsAB(a,SeriesOrder)

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


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


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



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


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



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

.
.
.

Amin · July 25th, 2022, 7:48 am

You can also find Z-series for densities with fourth cumulant smaller than zero. But as I mentioned in previous post, you have to specify some positive value for sixth cumulant.
You can specify some small skew with negative fourth cumulant.
Here are graphs of a few such densities
.
.

Amin · July 25th, 2022, 8:42 am

Later tomorrow, I will post a simple program that will show how to take arbitrary moments (that are meaningful and coherent) , standardize them, find Z-series of standardized density and then convert Z-series of Standardized density to Z-series of original density with initially specified arbitrary moments. It is very simple but I think a small program that does this will be very helpful.

Amin · July 25th, 2022, 12:37 pm

Friends, I will post my general program tomorrow. But briefly I want to mention that standardization used for central moments is exactly the same as standardization used for Z-series coefficients (we have to do some very simple accounting for mean). So once we specify a general density and calculate its central moments, we simply need to standardize the central moments and turn standardized moments into cumulants to fit the density. Then we can retrieve the original density's raw moments perfectly after scaling our Z-series with square-root of variance i.e. the same ratio that was used to standardize the density. So Z-series of original density is simply obtained by multiplying the Z-series of standardized density with square-root of variance (we have to make minor accounting for mean). This can become obvious when we see that nth raw moment is simply expectation of nth power of Z-series.

Amin · July 26th, 2022, 8:50 am

Friends, here is the general program which takes input raw moments and calculates its Z-series and also draws its graph. I have provided a small routine that takes a general mean and a general variance and standardized cumulants to create valid raw moments. If your density is surely valid, you can skip this routine and directly specify the raw moments. For reasonable values of raw moments that result in good standardized cumulants, the program should work reasonably well.
I have added comments to the program.
Here are the matlab functions.
.
.

function [] = CalculateZSeriesGeneralDensity06A()

%In this program, I first generate a valid density using arbitrary mean and variance and 
%standardized cumulants. This is done in order to make sure that input
%density is a vlaid density and standardized cumulants are not too off the
%possible values.
%This general density is generated with following command
%[Mu] = GenerateValidDensityRawMoments(mean,var,kk3,kk4,kk5,kk6)
%If you have a valid density, you can directly substitute Mu array with
%first six raw moments and skip this command and earlier part of the
%program.

%mean and variance of the general density.
mean=10;
var=3.0;

%Below are standardized cumulants that would be used to generate the general density
%Cumulants take the same standardization ratio as for respective central
%moments.
%Valid Density generation program takes these standardized cumulants and
%converts them to absolute cumulants using powers of var and then
%calculates raw moments from them.
kk1=0.0;
kk2=(1.0);
kk3=-.4;%1.0;
kk4=2;%-.50;
kk5=-3;%-.5;%15;
kk6=20;%10;%100.0;

%Mu below is the array carrying first six raw moments. you can directly
%replace it witn input raw moments if you are sure that your density is
%valid.
[Mu] = GenerateValidDensityRawMoments(mean,var,kk3,kk4,kk5,kk6);

%Below function converts raw moments to central moments and mean.
[mu1,cMu] = ConvertRawMomentsToCentralMoments(Mu,6);

%Below function converts central moments to absolute cumulants. First
%cumulant is kept at zero.
 K1=0;
 K2=cMu(2);
 K3=cMu(3);
 K4=cMu(4)-3*cMu(2).^2;
 K5=cMu(5)-10*cMu(3)*cMu(2);
 K6=cMu(6)-15*cMu(4)*cMu(2)-10*cMu(3).^2+30*cMu(2).^3;

 %Below we convert absolute cumulants to standardized cumulants.
 k1=0;
 k2=1;
 k3=K3/K2^1.5;
 k4=K4/K2^2.0;
 k5=K5/K2^2.5;
 k6=K6/K2^3.0;
 
%  k1
%  k2
%  
%  k3
%  kk3
%  
%  k4
%  kk4
%  
%  k5
%  kk5
%  
%  k6
%  kk6
%  
% str=input('Look at comparison of cumulants'); 

%Below is the program to find standardized Z-series coefficients from
%standardized cumulants.
[b0,b] = FindZSeriesGivenCumulantsStandard06(k1,k2,k3,k4,k5,k6);

%Below we convert standardized Z-series coefficients to Z-series
%coefficients that correspond to input raw moments.

a0=b0*sqrt(var);
a(1)=b(1)*sqrt(var);
a(2)=b(2)*sqrt(var);
a(3)=b(3)*sqrt(var);
a(4)=b(4)*sqrt(var);
a(5)=b(5)*sqrt(var);

%We add the mean to first coefficient. We had done earlier calculations
%with zero mean.
a0=a0+mean;

%Below we calculate the raw moments from Z-series coefficients.
SeriesOrder=5;
NMoments=6;
[Moments] = CalculateMomentsOfZSeries(a0,a,SeriesOrder,NMoments);

%Below we visually compare the Moments array that has raw momwnts calculated from Z-series
%coefficients with Mu array that contains raw moments input to the program.
Moments
Mu

str=input('Look at comparison of raw moments of Z-series density and input raw moments');

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

%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)=a0;
  for nn=1:SeriesOrder
      w(1:Nn)=w(1:Nn)+a(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


plot(w(wnStart+1:Nn-1),pw(wnStart+1:Nn-1),'r');

title(sprintf('Fm1 = %.2f,m1=%.2f,Fm2=%.2f,m2=%.2f,Fm3=%.2f,m3 =%.2f,Fm4=%.2f,m4=%.2f,Fm5=%.2f,m5=%.2f,Fm6=%.2f,m6=%.3f',Moments(1),Mu(1),Moments(2),Mu(2),Moments(3),Mu(3),Moments(4),Mu(4),Moments(5),Mu(5),Moments(6),Mu(6)));
%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
legend({'Newton Z-Series Density'},'Location','northeast')

a0
a
str=input('Look at Z-series Coefficients of General density with input raw Moments');


end

.
.
.

function [Mu] = GenerateValidDensityRawMoments(mean,var,k3,k4,k5,k6)
%This function takes mean,var and higher standardized cumulants. I am using
%standardized cumulants as input since it would be easier to realize when
%we are stretching the input density too much when we deal with increasing
%standardized cumulants.
%mean in parameters is true mean.
%var in parameters is true variance. not standardized.
%k3 is third standardized cumulant. 
%I checked first six Cumulants take the same standardization ratio
%as standardized central moments.
%k4 is fourth standardized cumulant.
%k5 is fifth standardized cumulant.
%k6 is sixth standardized cumulant.

mu1=mean;


k2=1;
cmu1=0;
cmu2=1;
cmu3=k3;
cmu4=3*k2^2+k4;
cmu5=10*k3*k2+k5;
cmu6=15*k2*k4+10*k3^2+15*k2^3+k6;

 cMu(1)=0;
 cMu(2)=cmu2*var;
 cMu(3)=cmu3*var^1.5;
 cMu(4)=cmu4*var^2;
 cMu(5)=cmu5*var^2.5;
 cMu(6)=cmu6*var^3;

 [Mu] = ConvertCentralMomentsToRawMoments(cMu,mu1,6);

%[mu1,cMu1] = ConvertRawMomentsToCentralMoments(Mu,6)
 
% cMu
% cMu1
 
 
end

,
,
,

function [Mu] = ConvertCentralMomentsToRawMoments(cu,mu,mOrder)

%cu
%mu
%mOrder
%str=input('Look at numbers');
Mu(1)=mu;
for mm=2:mOrder
    Mu(mm)=0;
    for jj=0:mm
        if(jj==0)
            Mu(mm)=Mu(mm)+factorial(mm)/factorial(jj)/factorial(mm-jj)*1*mu.^(mm-jj);
        else
            Mu(mm)=Mu(mm)+factorial(mm)/factorial(jj)/factorial(mm-jj)*cu(jj)*mu.^(mm-jj);
        end
    end
end



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 [b0,b] = FindZSeriesGivenCumulantsStandard06(k1,k2,k3,k4,k5,k6)


mu=0;
sigma=1.0;
% k1=0.0;
% k2=(1.0);
% k3=0;%1.0;
% k4=-1;%-.50;
% k5=0;%-.5;%15;
% k6=10;%10;%100.0;



[mu,sigma,a0,a1,a2,a3,a4,a5,mu1,rmu1,rmu2,rmu3,rmu4,rmu5,rmu6] = GenerateDensityFromCumulantsX06Simple(mu,sigma,k1,k2,k3,k4,k5,k6)

[median] = FindMedianHermiteDensity06(mu1,mu,sigma,a0,a1,a2,a3,a4,a5);
a6=0;
a7=0;
a8=0;


%[median] = FindMedianHermiteDensity08Simple(mu1,mu,sigma,a0,a1,a2,a3,a4,a5,a6,a7,a8)

Xm=median

%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))+ ...
    a6* (((Xm-mu)/sigma).^6 - 15* ((Xm-mu)/sigma).^4 + 45*((Xm-mu)/sigma).^2-15)+ ...
    a7*(((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma))+ ...
    a8*(((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105));

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) + ...
    a6* (((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma))+...
    a7*(((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105)+ ...
    a8*(((Xm-mu)/sigma).^9 - 36* ((Xm-mu)/sigma).^7 + 378*((Xm-mu)/sigma).^5-1260*((Xm-mu)/sigma).^3+945*((Xm-mu)/sigma)));

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))+ ...
    a6* (((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105) + ...
    a7*(((Xm-mu)/sigma).^9 - 36* ((Xm-mu)/sigma).^7 + 378*((Xm-mu)/sigma).^5-1260*((Xm-mu)/sigma).^3+945*((Xm-mu)/sigma))+ ...
    a8/720/56*(((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));


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)+ ...
    a6* (((Xm-mu)/sigma).^9 - 36* ((Xm-mu)/sigma).^7 + 378*((Xm-mu)/sigma).^5-1260*((Xm-mu)/sigma).^3+945*((Xm-mu)/sigma)) + ...
    a7*(((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)+ ...
    a8*(((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)));


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))+ ...
    a6* (((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 ) + ...
    a7*(((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))+ ...
    a8*(((Xm-mu)/sigma).^12 - 66* ((Xm-mu)/sigma).^10 + 1485*((Xm-mu)/sigma).^8-13860*((Xm-mu)/sigma).^6+51975*((Xm-mu)/sigma).^4-62370.*((Xm-mu)/sigma).^2+10395));




dpxdX
d2pxdX2
d3pxdX3
d4pxdX4

%str=input('Look at derivatives results')

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


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 

%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

cc0=c0;
cc(1)=c1;
cc(2)=c2;
cc(3)=c3;
cc(4)=c4;
cc(5)=c5;

SeriesOrder=5;
NMoments=6;


[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
rmu1
rmu2
rmu3
rmu4
rmu5
rmu6


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



%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,6)

da(1,1)=cc0;
da(2:SeriesOrder+1,1)=cc(1:SeriesOrder);


NoOfCumulants=NMoments

ObjBest=0;
 for mm=1:NoOfCumulants
     if(C(mm)~=0)
     ObjBest=ObjBest+F(mm,1).^2/(abs(C(mm)));
     end
 end

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

  ObjNew=0;
  for mm=1:NoOfCumulants
      if(C(mm)~=0)
      ObjBest=ObjBest+F(mm,1).^2/(abs(C(mm)));
      end
  end
  
  if(ObjBest<ObjNew)
       
   else
       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.
  
  
  
  [F,dF] = CalculateCumulantsAndDerivativesFromMoments_0(C,b0,b,SeriesOrder,SeriesOrder,NoOfCumulants);
  %[F,dF,uu] = CalculateCumulantsAndDerivativesFromMoments_0A(C,b0,b,SeriesOrder,SeriesOrder,NoOfCumulants);
 ObjLast=0;
 for mm=1:NoOfCumulants
     if(C(mm)~=0)
     ObjLast=ObjLast+F(mm,1).^2/(abs(C(mm)));
     end
 end
%ObjLast
  
%nn

%F


[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.
rmu1
rmu2
rmu3
rmu4
rmu5
rmu6
%str=input('Look at new moments');


%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),pdf(1:Mm),'g',X(1:Mm),pdfG(1:Mm),'b');
% 
% title(sprintf('Fm1 = %.2f,m1=%.2f,Fm2=%.2f,m2=%.2f,Fm3=%.2f,m3 =%.2f,Fm4=%.2f,m4=%.2f,Fm5=%.2f,m5=%.2f,Fm6=%.2f,m6=%.3f,b0=%.4f,b1=%.4f,b2=%.4f,b3=%.4f,b4=%.4f,b5=%.4f',Moments(1),rmu1,Moments(2),rmu2,Moments(3),rmu3,Moments(4),rmu4,Moments(5),rmu5,Moments(6),rmu6,b0,b(1),b(2),b(3),b(4),b(5)));
% %,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')
% 
% 
% b0
% b




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

.
.
.
There are other functions but I have posted them yesterday and you can copy them from my yesterday's post.

Here is the graph of the density that is generated by the above program as output.
.
.

Amin · July 26th, 2022, 10:29 am

I want to write this post especially to thank friends who made it possible for me to work on my research by protesting to mind control agencies about my persecution. Several years ago when I was largely on my own, it was extremely difficult for me to consistently work on my research. At that time, there would be drugging of water and food on a large scale every second or third week. And I would lose all momentum and would have to spend all my time trying to get good food or water or trying to recover from the effect of drugs that I had taken in with drugged food. It would be impossible to totally avoid the mind control drugs in food and water and I would usually be hit again and again with bad food and sometimes it would take quite a bit to recover from the effect of those drugs.
However the things did improve largely on average since past few years when friends noticed my research and protested to mind control agencies about my inhuman treatment. There were some periods of time when I would be able to get good food or water continuously for several months and I was able to utilize those periods to concentrate on my research. Though my persecution continues and still intensifies once in a while even now, overall conditions are much better than they used to be several years ago.
I am really grateful to friends for all the efforts they made to end my persecution. If it were not for protests and efforts of good people, I would not have been able to do any good research and I really owe it to good people because of whose efforts I was able to work reasonably successfully on interesting problems. I will take this opportunity to thank good Americans again and want to tell them that their good wishes, sympathies and efforts to end my human misery really means a lot to me and I am thankful to them for helping me at such a crucial point in my life when things had become extremely tough.

Amin · July 30th, 2022, 8:05 am

Friends, sorry that I have not posted for the past whole week. I have an extreme urgency to be able to make some money. Some money I receive from my parents has been the same for past ten years but prices of everything have increased by 5-6 times and it is almost impossible to survive. So I started working again on my market trading research programs.
I have previously mentioned that I had some ideas about how to initially calibrate the distribution to six moments and then add further moments one by one through optimization and Newton method. Though I am sure friends could have better ideas on that, I would still post my thoughts with equations in a day or two just to stimulate thought. I think the whole thing eventually has to go to Artificial Intelligence where standardized cumulants are input and Z-series coefficients are output and then output of AI is refined by a few quick iterations of Newton method.
I also want to tell friends that my persecution still continues. Mind control agents do not want to make a big deal about it but still continue to drug water at a lot of places where it would be likely for me to get water. They do not want to let me have a full mind like I have had in past few weeks and they are concerned about somehow controlling me. Water in various neighborhoods had improved but recently it has become very difficult for me to buy good bottled water in neighborhoods around. Once I take water from a public water filtration plant or some other random water-cooler outside of some random house, I know that I would not be able to get good water from the same place again. So I have to keep changing places where to get water. One good thing is that still they are not drugging the ground water otherwise it would be completely impossible for me to get good water easily.
Also when they know there is some food at home that I can possibly take, they usually drug it so I have to be very careful with taking food from home.
I have been mostly managing to avoid getting hit by being very careful but there is some urgency among mind control services that they do not want me to continue with my full brain for longer time and they are trying their best to drug my food somewhere and take some neurotransmitters from my brain that could take another week or more for me to recover and by that time, they would be prepared to do something else.
Another thing that I really want to mention is that they also continue to use gases especially when I am working on my research and this is very very annoying but mind control agents still use gases sometimes though I do not think they used gas during past week when I was sleeping.
I really want to request friends to please try to stop mind control agencies from persecuting me further and ask them to let me live with my brain intact.

Amin · July 30th, 2022, 12:01 pm

Friends, In the previous post I said that they were not using any gas on me while I would sleep. I was just trying to sleep to have a nap for an hour and a half but they started to use gas in my room. Probably they were emboldened because I had just mentioned that they were not doing anything like that and therefore they thought that I would simply stand it and not write anything about it. I really want to tell friends that they are doing almost everything they can to take good neurotransmitters out of my brain again. Please force mind control agencies to become civilized to good people.

Amin · August 3rd, 2022, 10:20 am

Friends I spent some time yesterday, rather unsuccessfully, trying to see if there could be some very simple way of extending the six moment Z-series representation to higher order in moments and coefficients of Z-series but very simple ideas did not seem to work.
I still have not tried ideas I had about some analytics that I wanted to share with friends. But starting today, I will devote two to three days to this project to see if I could systematically add a few moments and extend the Z-series construction program to to may be eighth or tenth order in moments, and seventh or ninth order in coefficients of time series.
So I hope to share the progress about this new program in a day or two. I hope it goes well.

Amin · August 4th, 2022, 12:55 pm

Friends, I am writing this post about gas in the Air-conditioners. It is summers and monsoon rainy season and it becomes hot and humid on many days. I had been trying to work in the drawing-room of our house where I would turn on the AC and then try to work. It worked well only for first few days and then later they would use some gas that would make me drowsy, sedated and lethargic. I still tried to continue to work in the same environment since it would be very difficult to work without air-conditioner but finally I gave up on that room and started working at other places. On many days I have very little productivity only because I cannot find a suitable environment to work anywhere.
Today I tried working in my mother's room after turning on air-conditioner. The first few hours were fine but they started releasing some gas there as well and the gas would make it very difficult for me to work there and then I simply had to turn off the AC and leave the room.
I do not use AC in my room since they have far more potent gases in AC of my room and it is totally impossible to survive in it for more than half to one hour and I start having extreme effect on my brain and my consciousness.

Amin · August 8th, 2022, 10:48 am

Friends I have been able to extend earlier six cumulant, up to fifth power Z-series representation to seven cumulants or sixth power Z-series. But I have not been able to add the eighth cumulant in calibration in a stable way. I will post the seven cumulant program with graphs later today but will have to think more about the strategy to tackle eighth cumulant that can possibly take up to a week.
I will be posting a few graphs and the new matlab program in a few hours. And also post some explanation of the new program tomorrow.

Amin · August 8th, 2022, 3:16 pm

Friends, I have tried plotting some graphs of densities with fitted cumulants and target standard cumulants along with the Z-series coefficients on the title. Not every combination of cumulants is possible. Here are the graphs. The blue density is that of standard normal.

.
.
.

Amin · August 8th, 2022, 4:46 pm

As a start, friends can try running the program with input cumulants as given in the above graphs in previous post.
,
,
,

function [b0,b] = FindZSeriesGivenCumulantsStandard07()


mu=0;
sigma=1.0;

 k1=0.0;
 k2=(1.0);
 k3=0;%1.0;
 k4=2.5;%-.50;
 k5=0.0;%-.5;%15;
 k6=150;%10;%100.0;
 k7=0;
% k8=1400;




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



[median] = FindMedianHermiteDensity06(mu1,mu,sigma,a0,a1,a2,a3,a4,a5);
a6=0;
a7=0;
a8=0;


%[median] = FindMedianHermiteDensity08Simple(mu1,mu,sigma,a0,a1,a2,a3,a4,a5,a6,a7,a8)

Xm=median

%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))+ ...
    a6* (((Xm-mu)/sigma).^6 - 15* ((Xm-mu)/sigma).^4 + 45*((Xm-mu)/sigma).^2-15)+ ...
    a7*(((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma))+ ...
    a8*(((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105));

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) + ...
    a6* (((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma))+...
    a7*(((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105)+ ...
    a8*(((Xm-mu)/sigma).^9 - 36* ((Xm-mu)/sigma).^7 + 378*((Xm-mu)/sigma).^5-1260*((Xm-mu)/sigma).^3+945*((Xm-mu)/sigma)));

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))+ ...
    a6* (((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105) + ...
    a7*(((Xm-mu)/sigma).^9 - 36* ((Xm-mu)/sigma).^7 + 378*((Xm-mu)/sigma).^5-1260*((Xm-mu)/sigma).^3+945*((Xm-mu)/sigma))+ ...
    a8/720/56*(((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));


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)+ ...
    a6* (((Xm-mu)/sigma).^9 - 36* ((Xm-mu)/sigma).^7 + 378*((Xm-mu)/sigma).^5-1260*((Xm-mu)/sigma).^3+945*((Xm-mu)/sigma)) + ...
    a7*(((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)+ ...
    a8*(((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)));


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))+ ...
    a6* (((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 ) + ...
    a7*(((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))+ ...
    a8*(((Xm-mu)/sigma).^12 - 66* ((Xm-mu)/sigma).^10 + 1485*((Xm-mu)/sigma).^8-13860*((Xm-mu)/sigma).^6+51975*((Xm-mu)/sigma).^4-62370.*((Xm-mu)/sigma).^2+10395));




dpxdX
d2pxdX2
d3pxdX3
d4pxdX4

%str=input('Look at derivatives results')

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


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 

%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

cc0=c0;
cc(1)=c1;
cc(2)=c2;
cc(3)=c3;
cc(4)=c4;
cc(5)=c5;

SeriesOrder=5;
NMoments=6;


[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
rmu1
rmu2
rmu3
rmu4
rmu5
rmu6


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



%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,6)

da(1,1)=cc0;
da(2:SeriesOrder+1,1)=cc(1:SeriesOrder);


NoOfCumulants=NMoments

ObjBest=0;
 for mm=1:NoOfCumulants
     if(C(mm)~=0)
     ObjBest=ObjBest+F(mm,1).^2/(abs(C(mm)));
     end
 end

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

  ObjNew=0;
  for mm=1:NoOfCumulants
      if(C(mm)~=0)
      ObjBest=ObjBest+F(mm,1).^2/(abs(C(mm)));
      end
  end
  
  if(ObjBest<ObjNew)
       
   else
       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.
  
  
  
  [F1,dF1] = CalculateCumulantsAndDerivativesFromMoments_0(C,b0,b,SeriesOrder,SeriesOrder,NoOfCumulants);
  %[F,dF,uu] = CalculateCumulantsAndDerivativesFromMoments_0A(C,b0,b,SeriesOrder,SeriesOrder,NoOfCumulants);
 ObjLast=0;
 for mm=1:NoOfCumulants
     if(C(mm)~=0)
     ObjLast=ObjLast+F(mm,1).^2/(abs(C(mm)));
     end
 end
%ObjLast
  
%nn

%F


[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.
rmu1
rmu2
rmu3
rmu4
rmu5
rmu6

str=input('Look at comparison of prior calibration');


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;


%For a rough calibration, I find the values of LambdaO, LambdaE and the
%coefficient a(6). %LambdaO is multiplier for odd coefficients. 
%LambdaE is multiplier for even coefficients
%I ask matlab to solve for seven equations (with constant coefficients)
%in least squares sense for three variables.LambdaO, LambdaE and the
%coefficient a(6). This is over-determined problem.
%Equation for each of seven moments is 
% dun_dLambadaO * DeltaLambadaO + dun_dLambadaE * DeltaLambadaE + dun_da6 *
% Deltaa6 = Defect_un
%Since first six moments are pre-calibrated, defect would almost always be
%in seventh moment equation.

a0=b0;
a(1:5)=b(1:5);
LambdaO=1;
LambdaE=1;
a_n=0;
NMoments=7;
SeriesOrder=5;
NZterms=6;
%Below I calculate defect and its derivatives with function below.
[F,dF] = CalculateCumulantsAndDerivativesFromMoments_Lambdas(C,a0,a,SeriesOrder,NZterms,NMoments,LambdaO,LambdaE,a_n)

% F
% F1
% dF
% dF1
% str=input('Look at F and df');

%Below I use a matlab specific command to solve the linear equations with
%given defect to come up with a starting guess for Newton method.
%x = lsqr(A,b)
x = lsqr(dF,F)

LambdaO=1.0;%LambdaO is multiplier for odd coefficients
LambdaE=1.0;%LambdaE is multiplier for even coefficients

LambdaO=LambdaO+x(1);
LambdaE=LambdaE+x(2);
a(6)=x(3);

a(1)=a(1)*LambdaO;
a(3)=a(3)*LambdaO;
a(5)=a(5)*LambdaO;

a0=a0*LambdaE;
a(2)=a(2)*LambdaE;
a(4)=a(4)*LambdaE;
%str=input('Look at X');

cc0=a0;
cc(1:6)=a(1:6);
%cc(6)=-.0002;

SeriesOrder=6;
NMoments=7;

%Below specify target cumulants for Newton iterations.
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,7)

da(1,1)=cc0;
da(2:SeriesOrder+1,1)=cc(1:SeriesOrder);


NoOfCumulants=NMoments

ObjBest=0;
 for mm=1:NoOfCumulants
     if(C(mm)~=0)
     ObjBest=ObjBest+F(mm,1).^2/sqrt(abs(C(mm)));
     end
 end

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

nn=0;
while((nn<30)&&((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);%.*Ratio;

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

[F,dF] = CalculateCumulantsAndDerivativesFromMoments_0(C,b0,b,SeriesOrder,7,NoOfCumulants);

  ObjNew=0;
  for mm=1:NoOfCumulants
      if(C(mm)~=0)
      ObjBest=ObjBest+F(mm,1).^2/sqrt(abs(C(mm)));
      end
  end
  
  if(ObjBest<ObjNew)
       
   else
       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.
rmu1
rmu2
rmu3
rmu4
rmu5
rmu6
rmu7
%rmu8
str=input('Look at new moments');

u1=Moments(1);
u2=Moments(2);
u3=Moments(3);
u4=Moments(4);
u5=Moments(5);
u6=Moments(6);
u7=Moments(7);

%Fitted cumulants below
Fk1=u1;
Fk2=u2-u1^2;
Fk3=u3-3*u2*u1+2*u1^3;
Fk4=u4-4*u3*u1-3*u2^2+12*u2*u1^2-6*u1^4;
Fk5=u5-5*u4*u1-10*u3*u2+20*u3*u1^2+30*u2^2*u1-60*u2*u1^3+24*u1^5;
Fk6=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;
Fk7= 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 ;





%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=91;  % 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 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('Fk1 = %.2f,k1=%.2f,Fk2=%.2f,k2=%.2f,Fk3=%.2f,k3 =%.2f,Fk4=%.2f,k4=%.2f,Fk5=%.2f,k5=%.2f,Fk6=%.2f,k6=%.3f,,Fk7=%.2f,k7=%.3f,b0=%.4f,b1=%.4f,b2=%.4f,b3=%.4f,b4=%.4f,b5=%.4f,b6=%.4f',Fk1,k1,Fk2,k2,Fk3,k3,Fk4,k4,Fk5,k5,Fk6,k6,Fk7,k7,b0,b(1),b(2),b(3),b(4),b(5),b(6)));
% %,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
%  
 legend({'Newton Density','Standard Gaussian Density'},'Location','northeast')
% 
% 
% b0
% b




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;
[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 [F,dF] = CalculateCumulantsAndDerivativesFromMoments_Lambdas(C,a0,a,SeriesOrder,NZterms,NMoments,LambdaO,LambdaE,a_n)
%[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;

NZterms=NZterms+1;
SeriesOrderNew=SeriesOrder+1;




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

EZ(1)=0;
EZ(2)=1;
for nn=3:NMoments*(SeriesOrder+1)+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


if(NMoments==7)
aLambda(1:SeriesOrder)=a;
a0Lambda=LambdaE*a0;
aLambda(2)=LambdaE*aLambda(2);
aLambda(4)=LambdaE*aLambda(4);


aLambda(1)=LambdaO*aLambda(1);
aLambda(3)=LambdaO*aLambda(3);
aLambda(5)=LambdaO*aLambda(5);

aLambda(6)=a_n;
end

if(NMoments==8)
aLambda(1:SeriesOrder)=a;
a0Lambda=LambdaE*a0;
aLambda(2)=LambdaE*aLambda(2);
aLambda(4)=LambdaE*aLambda(4);
aLambda(6)=LambdaE*aLambda(6);


aLambda(1)=LambdaO*aLambda(1);
aLambda(3)=LambdaO*aLambda(3);
aLambda(5)=LambdaO*aLambda(5);

aLambda(7)=a_n;
end


aLambda(SeriesOrderNew+1:NMoments*SeriesOrderNew+1)=0;
b0=a0Lambda;
b=aLambda;

for mm=1:NMoments
    if(mm>1)
        [b0,b] =SeriesProduct(a0Lambda,aLambda,b0,b,SeriesOrderNew*mm);
        b(SeriesOrderNew*mm+1:NMoments*SeriesOrderNew+1)=0;
    end
   % b0
   % b
%str=input('Look at numbers')    
    EXZ(mm+1,1)=b0;
    for pp2=1:SeriesOrderNew*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:SeriesOrderNew*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);

u5=EXZ(6,1);
u6=EXZ(7,1);
u7=EXZ(8,1);
if(NMoments==8)
u8=EXZ(9,1);
end
%u4=EXZ(5,1);

if(NMoments==7)
for mm=2:SeriesOrderNew+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);
end



%index 2 below corresponds to a(1)
ddu1(1)=a(1)*du1(2)+a(3)*du1(4)+a(5)*du1(6);
ddu2(1)=a(1)*du2(2)+a(3)*du2(4)+a(5)*du2(6);
ddu3(1)=a(1)*du3(2)+a(3)*du3(4)+a(5)*du3(6);
ddu4(1)=a(1)*du4(2)+a(3)*du4(4)+a(5)*du4(6);
ddu5(1)=a(1)*du5(2)+a(3)*du5(4)+a(5)*du5(6);
ddu6(1)=a(1)*du6(2)+a(3)*du6(4)+a(5)*du6(6);
ddu7(1)=a(1)*du7(2)+a(3)*du7(4)+a(5)*du7(6);


ddu1(2)=aa0*1+a(2)*1+a(4)*3;
ddu2(2)=     a(2)*du2(3)+a(4)*du2(5);
ddu3(2)=     a(2)*du3(3)+a(4)*du3(5);
ddu4(2)=     a(2)*du4(3)+a(4)*du4(5);
ddu5(2)=     a(2)*du5(3)+a(4)*du5(5);
ddu6(2)=     a(2)*du6(3)+a(4)*du6(5);
ddu7(2)=     a(2)*du7(3)+a(4)*du7(5);

ddu1(3)=15;
ddu2(3)=du2(7);
ddu3(3)=du3(7);
ddu4(3)=du4(7);
ddu5(3)=du5(7);
ddu6(3)=du6(7);
ddu7(3)=du7(7);

end



if(NMoments==8)
for mm=2:SeriesOrderNew+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



%index 2 below corresponds to a(1)
ddu1(1)=a(1)*du1(2)+a(3)*du1(4)+a(5)*du1(6);
ddu2(1)=a(1)*du2(2)+a(3)*du2(4)+a(5)*du2(6);
ddu3(1)=a(1)*du3(2)+a(3)*du3(4)+a(5)*du3(6);
ddu4(1)=a(1)*du4(2)+a(3)*du4(4)+a(5)*du4(6);
ddu5(1)=a(1)*du5(2)+a(3)*du5(4)+a(5)*du5(6);
ddu6(1)=a(1)*du6(2)+a(3)*du6(4)+a(5)*du6(6);
ddu7(1)=a(1)*du7(2)+a(3)*du7(4)+a(5)*du7(6);
ddu8(1)=a(1)*du8(2)+a(3)*du8(4)+a(5)*du8(6);


ddu1(2)=aa0*0+a(2)*1+a(4)*3 +a(6)*15;
ddu2(2)=     a(2)*du2(3)+a(4)*du2(5)+a(6)*du2(7);
ddu3(2)=     a(2)*du3(3)+a(4)*du3(5)+a(6)*du3(7);
ddu4(2)=     a(2)*du4(3)+a(4)*du4(5)+a(6)*du4(7);
ddu5(2)=     a(2)*du5(3)+a(4)*du5(5)+a(6)*du5(7);
ddu6(2)=     a(2)*du6(3)+a(4)*du6(5)+a(6)*du6(7);
ddu7(2)=     a(2)*du7(3)+a(4)*du7(5)+a(6)*du7(7);
ddu8(2)=     a(2)*du8(3)+a(4)*du8(5)+a(6)*du8(7);

ddu1(3)=0;
ddu2(3)=du2(8);
ddu3(3)=du3(8);
ddu4(3)=du4(8);
ddu5(3)=du5(8);
ddu6(3)=du6(8);
ddu7(3)=du7(8);
ddu8(3)=du8(8);

end




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

%k1=LambdaE*aa0+u1;
k1=aa0+u1;


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;

k8= + u8 - 8* u1* u7 -7 *( 4* u2* (-20* u3^2 + u6) ...
    -8 *u1^2* (180* u2^3 - 30* u3^2 - 45* u2* u4 + u6) ...
    + 16* u1* (45* u2^2* u3 - 5 *u3* u4 - 3* u2* u5) ...
    +48* u1^3 *(-40* u2* u3 + u5) ...
    +240* u1^4 *(15* u2^2 - u4) - 60* u2^2* u4 + 5* u4^2 + 8* u3 *u5 ...
    + 720* u1^8 - 2880* u1^6 *u2 + 90* u2^4 + 960* u1^5* u3);

%str=input('Look at comparison of k8');

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


dk1(1)=ddu1(1);
dk1(2)=ddu1(2);
dk1(3)=ddu1(3);



% 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
for mm=1:3
%dk1(mm)=du1(mm);
dk2(mm)=-2.* u1.*ddu1(mm)+ddu2(mm);
dk3(mm)=ddu3(mm)-3*ddu1(mm)*u2-3*u1*ddu2(mm)+2*3*u1^2*ddu1(mm);
dk4(mm)=ddu4(mm)-4*ddu3(mm)*u1-4*u3*ddu1(mm)-3*2*u2*ddu2(mm)+12*ddu2(mm)*u1^2+12*2*u2*u1*ddu1(mm)-6*4*u1^3*ddu1(mm);

if(NMoments>=5)
dk5(mm)=ddu5(mm)-5*ddu4(mm)*u1-5*u4*ddu1(mm)-10*ddu3(mm)*u2-10*u3*ddu2(mm)+20*ddu3(mm)*u1^2+20*2*u3*u1*ddu1(mm)+ ...
    30*2*u2*ddu2(mm)*u1+30*u2^2*ddu1(mm)-60*ddu2(mm)*u1^3-60*3*u2*u1^2*ddu1(mm)+24*5*u1^4*ddu1(mm);
end
if(NMoments>=6)
dk6(mm)=ddu6(mm)-6*ddu5(mm)*u1-6*u5*ddu1(mm)-15*ddu4(mm)*u2-15*u4*ddu2(mm)+30*ddu4(mm)*u1^2+30*2*u4*u1*ddu1(mm)-10*2*u3*ddu3(mm)+ ...
    120*ddu3(mm)*u2*u1-120*u3*ddu2(mm)*u1-120*u3*u2*ddu1(mm)-120*ddu3(mm)*u1^3-120*3*u3*u1^2*ddu1(mm)+30*3*u2^2*ddu2(mm) - ...
    270*2*u2*ddu2(mm)*u1^2+270*2*u2^2*u1*ddu1(mm)+360*ddu2(mm)*u1^4+360*4*u2*u1^3*ddu1(mm)-120*6*u1^5*ddu1(mm);
end
if(NMoments>=7)
    dk7(mm)= ddu7(mm)- 7 *ddu1(mm)* u6- 7 *u1* ddu6(mm) + 42*2* ddu1(mm)* u5+ 42* u1^2* ddu5(mm)  ...
        - 21* ddu2(mm) *u5- 21* u2 *ddu5(mm)- 210* 3*u1^2*ddu1(mm) *u4- 210* u1^3 *ddu4(mm) + ...
 210* ddu1(mm)*u2* u4+210* u1*ddu2(mm)* u4+210* u1*u2* ddu4(mm) - 35* ddu3(mm) *u4- 35* u3 *ddu4(mm) ...
 + 840 *4*u1^3*ddu1(mm)* u3 + 840 *u1^4* ddu3(mm)+ 210* 2*u2*ddu2(mm)* u3 + 210* u2^2* ddu3(mm)+ ...
 210*12* 3*u1^2*ddu1(mm)* u2^2 +210*12* u1^3* 2*u2*ddu2(mm)   - 1260* 2*u1*ddu1(mm) * u2* u3- 1260* u1^2 * ddu2(mm)* u3 ...
 - 1260* u1^2 * u2* ddu3(mm) - 630 *ddu1(mm)* u2^3- 630 *u1* 3*u2^2*ddu2(mm) ...
 +140*ddu1(mm)* u3^2+140*u1* 2*u3*ddu3(mm)  - 2520* 5*u1^4*ddu1(mm)* u2- 2520* u1^5* ddu2(mm) + 720* 7*u1^6*ddu1(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) - ddu1(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);


dk8(mm)= + ddu8(mm) - 8* ddu1(mm)* u7- 8* u1* ddu7(mm) -7 *( 4* ddu2(mm)* (-20* u3^2 + u6) +4* u2* (-20* 2*u3*ddu3(mm) + ddu6(mm)) ...
    -8 *2*u1*ddu1(mm)* (180* u2^3 - 30* u3^2 - 45* u2* u4 + u6)-8 *u1^2* (180*3* u2^2*ddu2(mm) - 30* 2*u3*ddu3(mm) - 45* ddu2(mm)* u4- 45* u2* ddu4(mm) + ddu6(mm)) ...
    + 16* ddu1(mm)*(45* u2^2* u3 - 5 *u3* u4 - 3* u2* u5) + 16* u1*(45* 2*u2*ddu2(mm)* u3+45* u2^2* ddu3(mm) - 5 *ddu3(mm)* u4- 5 *u3* ddu4(mm) - 3* ddu2(mm)* u5- 3* u2* ddu5(mm))  ...
    +48* 3*u1^2*ddu1(mm) *(-40* u2* u3 + u5)+48* u1^3 *(-40* ddu2(mm)* u3-40* u2* ddu3(mm) + ddu5(mm)) ...
    +240* 4*u1^3*ddu1(mm) *(15* u2^2 - u4)+240* u1^4 *(15* 2*u2*ddu2(mm) - ddu4(mm)) ...
    - 60* 2*u2*ddu2(mm)* u4 - 60* u2^2* ddu4(mm) + 5* 2*u4*ddu4(mm) + 8* ddu3(mm) *u5+ 8* u3 *ddu5(mm) ... 
    + 720* 8*u1^7*ddu1(mm) - 2880* 6*u1^5*ddu1(mm) *u2- 2880* u1^6 *ddu2(mm) + 90* 4*u2^3*ddu2(mm) + 960* 5*u1^4*ddu1(mm)* u3+ 960* u1^5* ddu3(mm)); 

%k8= + u8 - 8* u1* u7 -7 *( 4* u2* (-20* u3^2 + u6) ...
%    -8 *u1^2* (180* u2^3 - 30* u3^2 - 45* u2* u4 + u6) ...
%    + 16* u1* (45* u2^2* u3 - 5 *u3* u4 - 3* u2* u5) ...
%    +48* u1^3 *(-40* u2* u3 + u5) ...
%    +240* u1^4 *(15* u2^2 - u4) - 60* u2^2* u4 + 5* u4^2 + 8* u3 *u5 ...
%    + 720* u1^8 - 2880* u1^6 *u2 + 90* u2^4 + 960* u1^5* u3)



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
     for mm=1:3
    
        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

,
,
,

Amin · August 9th, 2022, 8:41 pm

Friends, I have mentioned several times that I want to start algorithmic trading in vanilla stocks. Several months have passed since I opened a trading account with Interactive Brokers. But unfortunately I could not fund my account due to foreign exchange restrictions in Pakistan. IB would not accept money in my account from a third party name (and therefore I cannot ask a foreign friend or acquaintance to wire money to my account) and as a Pakistani resident only way to fund my account is through a wire transfer. We have always found our country to be in foreign exchange crisis since our childhood. Therefore Pakistani banks are not allowed to make wire transfers from individuals to foreign businesses. I therefore looked for possibilities to open a bank account in a foreign country since I can make wire transfers to (individuals) my own account abroad. Most international banking requires large deposits. I recently found a relatively cheaper option with no balance requirements to open an international account with Caye Bank in small Central American country of Belize and I am putting together documents to apply for that bank account. The bank is supposedly good and safe.
But in past my experience with banks has been extremely poor. A very respected British bank stole my thousand pounds on insistence of CIA and mind control agencies. When I was in UK, a very reputed Pakistani bank also connived with mind control agencies to clear my check (for USD 2500) in past back dates despite that my check had already returned (and the guy I wrote the check had told me that check had returned but I found that check had cleared when I needed money in UK after several months) since my account was dormant and the check should never had cleared until I activated my bank account by depositing money in my account myself. Since when I was in UK, I spent a lot of money to stay in good hotels and move around in taxis as I was worried about my safety and also wanted to have good food. And therefore mind control agencies connived with both British and Pakistani banks to steal my money since they wanted to bleed me and make it difficult for me to keep spending for a longer time. I want to tell friends that I was working with a Japanese company from Pakistan before coming to UK on Tier one visa (in Sep 2010) and was making more than hundred thousand dollars per year despite being on very extremely high antipsychotics. And my saving gave mind control agencies a very tough time in UK as it became very difficult for them to drug my food despite drugging entire neighborhoods in London and other cities as I moved around a lot. And since my savings gave mind control agencies a very tough time, it was one reason mind control agencies never allowed me to make money after I returned to Pakistan from UK after a stint of seven months there.
I was trying to play with my non-funded account with IB but I cannot even download historical data without data subscriptions that require money in my IB account.
Good thing is that I would still have some spare time to think about improving the Newton method for densities and I want to try to improve the current program and take the Z-series representation of densities to eighth moment. I also plan to spend some time playing with gamma density and Laguerre polynomials if it could be possible to write a program similar to normal distribution based densities. I hope that I would be able to spend my spare time to do something interesting that friends would also like.

Amin · Yesterday, 12:50 pm

Friends, today I was thinking of some other possible and general way to find a good guess for Newton method in the context for moment/cumulant calibration. I think I have some new workable ideas. I am trying to work with a new method which is a combination of Taylor series and matrix manipulations in order to find a very good guess for Newton method. Though I am very optimistic, I still cannot be sure whether the new ideas would actually work in practice. I will be working with this new possible method in next few days and would let friends know about progress on the method in next 2-3 days. At this point I am working on it but there is still no guarantee that everything would fit together as I plan it to be. I will let friends know in a few days how things go and will post a new program if things go well.

Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation