Friends, I want to break the good news that I have been able to fit the moments of lognormal and CEV type (in SDE jargon) random variables(in original coordinates) to coefficients of a Z-series. The fit to moments is exact when they are derived from a real distribution. The program works very well.
One caveat is that when volatility is very high, even though fit to first prescribed eight moments is exact in my algorithm, there are times when actual SDE distribution is different from the Z-series distribution even though both fit the same first eight moments. This is because when volatility is very high and tails of SDE distribution are very large, it is difficult to pin it down with first eight moments.
For the kind of volatilities we see for stocks (<40% per year), this algorithm will work with precision out to more than five years. For very high volatilities, you will have to prescribe more than eight moments to pin down the distribution and the algorithm will work well. Volatility (<40% per year) is not a limitation of the algorithm, it is a limitation to get exact fit to distribution with first eight moments (as I have in current matlab implementation of the algorithm).
I am very sure the algorithm would work perfectly well (with precision) for fat-tailed distributions of the kind we come across in asset management.
The new program is not based on inverse series, rather it is a modification of the program that I used for solution of SDEs a few weeks earlier. I made the same algorithm more precise and it works very well for CEV and lognormal type SDE distributions when volatility is not extremely high. When it does not work, it still fits the prescribed moments perfectly well but, in those cases, we need more than eight moments to pin down the distribution.
The algorithm is general and can easily be used to find Z-series coefficients to any number of moments and fit to first eight moments is not a limitation of the algorithm (but of the current program).
The algorithm is based on series of one dimensional Newton iterations and fits the specified distribution with precision whenever specified moments are taken from a real distribution (I got them from CEV and Lognromal SDEs in original coordinates).
I will post the new program later tonight or possibly tomorrow.
Serving the Quantitative Finance Community
Re: Breakthrough in the theory of stochastic differential equations and their simulation
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, I have written the Z-series probability distribution calculation from moments programs in the form as independent functions. I will be posting them tomorrow after making minor changes and writing comments.
I have also written functions that calculate the probability density function, cumulative density function and value of Z once you input the value of X.
[$]X=c0 + c1 Z+ c2 Z^2+ c3 Z^3+ c4 Z^4+ c5 Z^5+ c6 Z^6+ c7 Z^7[$]
Since in our framework given above, the value of X is known after value of Z is specified. But mostly we want to know pdf and cdf directly from value of X and we do not know what is the corresponding value of Z. I have done exactly that in a new function that inverts the above Z-series and finds Z given a value of X. As we know inverse series is only approximate, this initial (but very close) guess is fed to Newton method to calculate precise value of Z that solves for X. After finding appropriate Z, we can easily calculate probability density function and cumulative density function of specific input value of X associated with calculated Z.
In case, you derived your Z-series from some analytic distribution formula using derivatives, the above formulas will calculate the same precise values of cdf and pdf that you would find from your closed form formulas of the distribution.
I will be posting my programs tomorrow.
I have also written functions that calculate the probability density function, cumulative density function and value of Z once you input the value of X.
[$]X=c0 + c1 Z+ c2 Z^2+ c3 Z^3+ c4 Z^4+ c5 Z^5+ c6 Z^6+ c7 Z^7[$]
Since in our framework given above, the value of X is known after value of Z is specified. But mostly we want to know pdf and cdf directly from value of X and we do not know what is the corresponding value of Z. I have done exactly that in a new function that inverts the above Z-series and finds Z given a value of X. As we know inverse series is only approximate, this initial (but very close) guess is fed to Newton method to calculate precise value of Z that solves for X. After finding appropriate Z, we can easily calculate probability density function and cumulative density function of specific input value of X associated with calculated Z.
In case, you derived your Z-series from some analytic distribution formula using derivatives, the above formulas will calculate the same precise values of cdf and pdf that you would find from your closed form formulas of the distribution.
I will be posting my programs tomorrow.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends here are the matlab functions required for Z-series density calibration to first eight raw moments.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
When you run the program, you will see a comparison of input moments (rMu) and fitted moments.
rMu =
1.0e+05 *
Columns 1 through 7
0.000012499863483 0.000021197093700 0.000050274347857 0.000171368317209 0.000856094970826 0.006320436615909 0.068698344472325
Column 8
1.084082326087001
FittedMoments =
1.0e+05 *
Columns 1 through 7
0.000012499863483 0.000021197093700 0.000050274347857 0.000171368317209 0.000856094970826 0.006320436615909 0.068698344472325
Column 8
1.084082326087000
Look at comparison of fitted moments with input moments.
And you will see the following graph of the density.
.
.
Code: Select all
function [c0,c,FittedMoments] = CalculateZSeriesDensityFromRawMoments02(rMu)
%The above function takes an array of first eight raw moments as input.
%And returns Z-series Coeffs upto 7th power of Z as in
%Y=c0+c(1)*Z+c(2)*Z^2+c(3)*Z^3+c(4)*Z^4+c(5)*Z^5+c(6)*Z^6+c(7)*Z^7
%Fitted Moments are raw moments of the Z-series density and usually are exactly
%the same as input raw moments rMu(1:8).i.e FittedMoments(1:8)==rMu(1:8)
%with high precision.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%If you do not have valid moments from a density available, please
%uncomment the following input moments to run the program.
rMu(1)=0.000012499863483*10^5;
rMu(2)=0.000021197093700*10^5;
rMu(3)=0.000050274347857*10^5;
rMu(4)=0.000171368317209*10^5;
rMu(5)=0.000856094970826*10^5;
rMu(6)=0.006320436615909*10^5;
rMu(7)=0.068698344472325*10^5;
rMu(8)=1.084082326087001*10^5;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
SeriesOrder=7; %Z-series till 7th power. hardcoded here cannot chagne.
mOrder=8; %Number of raw moments in calibration. hardcoded here cannot chagne.
NMoments=mOrder;
%Find central Moments from raw moments
[Mu1,cMu] = ConvertRawMomentsToCentralMoments(rMu,mOrder);
%Find standardized Moments from central moments
sMu(1)=0;
sMu(2)=1;
sMu(3)=cMu(3)/cMu(2).^1.5;
sMu(4)=cMu(4)/cMu(2).^2.0;
sMu(5)=cMu(5)/cMu(2).^2.5;
sMu(6)=cMu(6)/cMu(2).^3.0;
sMu(7)=cMu(7)/cMu(2).^3.5;
sMu(8)=cMu(8)/cMu(2).^4.0;
iter=5; %Number of iterations in the calibration.
%Please experiment and play around.
bGuess(2:7)=0; %We have assigned initial values of coefficients
bGuess(1)=1; %other than c(1) to zero. We are starting with no prior information.
[c0,c] = PreSmoothingGuessAdvancedFromGuessBestNew(sMu,bGuess,iter);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Comment this block if you do not want multidimensional Newton pass
%Even though preSmoothing guess is very close to optimal, slight
%inaccuracies can change the shape of probabioity density so this
%Multidimensional Newton is usually required.
[F,dF] = CalculateMomentsAndDerivatives_0(sMu,c0,c,SeriesOrder,SeriesOrder,NMoments);
da(1,1)=c0;
da(2:SeriesOrder+1,1)=c(1:SeriesOrder);
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%Replace with your own more intelligent objective function if you like.
ObjBest=100000*(abs(sMu(1)-Moments(1)))+abs(sMu(2)-Moments(2))+abs((sMu(3)-Moments(3))^(1/1.5))+abs((sMu(4)-Moments(4))^(1/2.0)) + ...
abs((sMu(5)-Moments(5))^(1/2.0))+abs((sMu(6)-Moments(6))^(1/2.0))+abs((sMu(7)-Moments(7))^(1/2.0))+ ...
abs((sMu(8)-Moments(8))^(1/2.0));
b0Best=c0;
bBest(1:SeriesOrder)=c(1:SeriesOrder);
nn=0;
while((nn<20)&&((abs(F(1,1))>.000000000001) || (abs(F(2,1))>.000000000001) || (abs(F(3,1))>.000000000001) || (abs(F(4,1))>.00000000001)|| (abs(F(5,1))>.00000000001)|| (abs(F(6,1))>.00000000001)|| (abs(F(7,1))>.00000000001)|| (abs(F(8,1))>.00000000001) ))
nn=nn+1;
%Below Newton matrix equation to improve the Z-series coefficients guess at previous step.
da=da-dF\F;
%b0=median
b0=da(1,1);
b(1:SeriesOrder)=da(2:SeriesOrder+1,1);
%[F,dF] = CalculateCumulantsAndDerivativesFromMoments_0(C,b0,b,SeriesOrder,SeriesOrder,NoOfCumulants);
[F,dF] = CalculateMomentsAndDerivatives_0(sMu,b0,b,SeriesOrder,SeriesOrder,NMoments);
[IsValidFlag] = CheckIsValidDensity(b0,b);
[Moments] = CalculateMomentsOfZSeries(b0,b,SeriesOrder,NMoments);
ObjNew=100000*(abs(sMu(1)-Moments(1)))+abs(sMu(2)-Moments(2))+abs((sMu(3)-Moments(3))^(1/1.5))+abs((sMu(4)-Moments(4))^(1/2.0)) + ...
abs((sMu(5)-Moments(5))^(1/2.0))+abs((sMu(6)-Moments(6))^(1/2.0))+abs((sMu(7)-Moments(7))^(1/2.0))+ ...
abs((sMu(8)-Moments(8))^(1/2.0));
if((ObjBest>ObjNew) &&( IsValidFlag))
ObjBest=ObjNew;
b0Best=b0;
bBest(1:SeriesOrder)=b(1:SeriesOrder);
end
da(1,1)=b0;
da(2:SeriesOrder+1,1)=b(1:SeriesOrder);
end
c0=b0Best;%Best;
c(1:SeriesOrder)=bBest;%Best(1:SeriesOrder);
%Above are the best chosen Z-series coefficients based on the objective function.
%But my objective function can be improved.
%Multidimensional Newton root search ends. Comment till here if you like to
%exclude this
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Make sure you do not comment the lines below
c0=c0*sqrt(cMu(2)); %convert from standardized to central moments
c=c*sqrt(cMu(2)); %convert from standardized to central moments
c0=c0+Mu1; %convert from central to raw moments by adding mean.
[FittedMoments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Please Uncomment this if you want to see a comparison of input moments and
%fitted moments
rMu
FittedMoments
str=input('Look at comparison of fitted moments with input moments');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%COMMENT lines below if you do not want to see the density of new Z-series random variable with
%fitted moments.
%First calculate normal random variable on a grid below
dNn=.1/2; % Normal density subdivisions width. would change with number of subdivisions
Nn=45*4; % No of normal density subdivisions
NnMid=((1+Nn)/2)*dNn;
Z(1:Nn)=(((1:Nn)*dNn)-NnMid);
%Z
%str=input('Look at Z');
%Now calculate Y as a function of normal random variable.
Y(1:Nn)=c0;
for nn=1:7
Y(1:Nn)=Y(1:Nn)+c(nn)*Z(1:Nn).^nn;
end
%Now take change of densities derivative of Y with respect to normal
DfY(1:Nn)=0;
for nn=2:Nn-1
DfY(nn) = (Y(nn + 1) - Y(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
DfY(Nn)=DfY(Nn-1);
DfY(1)=DfY(2);
%Now calculate the density of Y from density of normal random variable
%using change of probability derivative.
pY(1:Nn)=0;
for nn = 1:Nn
pY(nn) = (normpdf(Z(nn),0, 1))/abs(DfY(nn));
end
plot(Y(1:Nn),pY(1:Nn),'b');
title('Density of the Z-Series Variable with Fitted Moments');
%str=input('Look at the density of Z-series random variable')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end.
.
Code: Select all
function [c0,c] = PreSmoothingGuessAdvancedFromGuessBestNew(cmu,cin,iter)
%Mul=1.0;
SeriesOrder=7;
NMoments=8;
EZ(1)=0;
EZ(2)=1;
for nn=3:NMoments*SeriesOrder+SeriesOrder+2
if rem(nn,2)==1
EZ(nn)=0;
else
EZ(nn)=EZ(nn-2)*(nn-1);
EZ(nn);
end
end
c0=-(cin(2)+3*cin(4)+15*cin(6)); %This is for zero mean condition.
c=cin;
MaxIter=10;
for nn=1:iter %increase the iterations over coefficients if neeeded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
tolerance=.0000001*cmu(3);
mm=0;
while ((abs(M3-cmu(3)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(2)=c(2)-(M3-cmu(3))/dc2;
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,3,2,EZ);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
tolerance=.0000001*cmu(4);
mm=0;
while ((abs(M4-cmu(4)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(3)=c(3)-(M4-cmu(4))/dc3;
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,4,2,EZ);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
tolerance=.0000001*cmu(5);
mm=0;
while ((abs(M5-cmu(5)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(4)=c(4)-(M5-cmu(5))/dc4;
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,5,2,EZ);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
tolerance=.0000001*cmu(6);
mm=0;
while ((abs(M6-cmu(6)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(5)=c(5)-(M6-cmu(6))/dc5;
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,6,2,EZ);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M7,dc6] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,7,EZ);
tolerance=.0000001*cmu(7);
mm=0;
while ((abs(M7-cmu(7)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(6)=c(6)-(M7-cmu(7))/dc6;
[M7,dc6] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,7,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,7,2,EZ);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M8,dc7] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,8,EZ);
tolerance=.0000001*cmu(8);
mm=0;
while ((abs(M8-cmu(8)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(7)=c(7)-(M8-cmu(8))/dc7;
[M8,dc7] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,8,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,8,2,EZ);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
.
.
Code: Select all
function [EnMoment,dEnMomentdCoeff] = CalculateParticularMomentAndDerivativeOfItsCoeff(a0,a,SeriesOrder,nMoment,EZ)
% if(SeriesOrder>=8)
% a0=-(a(2)+3*a(4)+15*a(6)+105*a(8));% ---1
% end
% if(SeriesOrder<8)
% a0=-(a(2)+3*a(4)+15*a(6));% ---1
% end
%
% EZ(1)=0;
% EZ(2)=1;
% for nn=3:NMoments*SeriesOrder+NZterms+2
% if rem(nn,2)==1
% EZ(nn)=0;
% else
% EZ(nn)=EZ(nn-2)*(nn-1);
% EZ(nn);
% end
% end
% EZ;
%EXZ(1,1)=1;
%for pp1=1:NZterms
% EXZ(1,pp1+1)=EZ(pp1);
%end
a(SeriesOrder+1:nMoment*SeriesOrder+1+SeriesOrder)=0;
b0=a0;
b=a;
for mm=2:nMoment
[b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*mm+SeriesOrder+1);
b(SeriesOrder*mm+1:nMoment*SeriesOrder+SeriesOrder+1)=0;
if(mm==nMoment-1)
dEnMomentdCoeff=b0.*EZ(nMoment-1);
for pp2=1:SeriesOrder*nMoment-1
dEnMomentdCoeff=dEnMomentdCoeff+b(pp2).*EZ(pp2+nMoment-1);
end
dEnMomentdCoeff=dEnMomentdCoeff*nMoment;
end
if(mm==nMoment)
EnMoment=b0;
for pp2=1:SeriesOrder*nMoment
EnMoment=EnMoment+b(pp2).*EZ(pp2);
end
end
end
end
.
.
Code: Select all
function [c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,nMoment,iter,EZ)
MaxIter=7;
SeriesOrder=7;
c0=-(c(2)+3*c(4)+15*c(6));
if(nMoment==3)
for nn=1:iter %increase the iterations over coefficients if neeeded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
tolerance=.0000001*cmu(3);
mm=0;
while ((abs(M3-cmu(3)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(2)=c(2)-(M3-cmu(3))/dc2;
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
if(nMoment==4)
for nn=1:iter %increase the iterations over coefficients if neeeded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
tolerance=.0000001*cmu(3);
mm=0;
while ((abs(M3-cmu(3)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(2)=c(2)-(M3-cmu(3))/dc2;
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
tolerance=.0000001*cmu(4);
mm=0;
while ((abs(M4-cmu(4)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(3)=c(3)-(M4-cmu(4))/dc3;
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(nMoment==5)
for nn=1:iter %increase the iterations over coefficients if neeeded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
tolerance=.0000001*cmu(3);
mm=0;
while ((abs(M3-cmu(3)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(2)=c(2)-(M3-cmu(3))/dc2;
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
tolerance=.0000001*cmu(4);
mm=0;
while ((abs(M4-cmu(4)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(3)=c(3)-(M4-cmu(4))/dc3;
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
tolerance=.0000001*cmu(5);
mm=0;
while ((abs(M5-cmu(5)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(4)=c(4)-(M5-cmu(5))/dc4;
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(nMoment==6)
for nn=1:iter %increase the iterations over coefficients if neeeded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
tolerance=.0000001*cmu(5);
mm=0;
while ((abs(M3-cmu(3)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(2)=c(2)-(M3-cmu(3))/dc2;
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
tolerance=.0000001*cmu(5);
mm=0;
while ((abs(M4-cmu(4)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(3)=c(3)-(M4-cmu(4))/dc3;
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
tolerance=.0000001*cmu(5);
mm=0;
while ((abs(M5-cmu(5)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(4)=c(4)-(M5-cmu(5))/dc4;
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
tolerance=.0000001*cmu(6);
mm=0;
while ((abs(M6-cmu(6)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(5)=c(5)-(M6-cmu(6))/dc5;
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(nMoment==7)
for nn=1:iter %increase the iterations over coefficients if neeeded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
tolerance=.0000001*cmu(3);
mm=0;
while ((abs(M3-cmu(3)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(2)=c(2)-(M3-cmu(3))/dc2;
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
tolerance=.0000001*cmu(4);
mm=0;
while ((abs(M4-cmu(4)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(3)=c(3)-(M4-cmu(4))/dc3;
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
tolerance=.0000001*cmu(5);
mm=0;
while ((abs(M5-cmu(5)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(4)=c(4)-(M5-cmu(5))/dc4;
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
tolerance=.0000001*cmu(6);
mm=0;
while ((abs(M6-cmu(6)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(5)=c(5)-(M6-cmu(6))/dc5;
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M7,dc6] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,7,EZ);
tolerance=.0000001*cmu(7);
mm=0;
while ((abs(M7-cmu(7)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(6)=c(6)-(M7-cmu(7))/dc6;
[M7,dc6] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,7,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if(nMoment==8)
for nn=1:iter %increase the iterations over coefficients if neeeded
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
tolerance=.0000001*cmu(3);
mm=0;
while ((abs(M3-cmu(3)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(2)=c(2)-(M3-cmu(3))/dc2;
[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
tolerance=.0000001*cmu(4);
mm=0;
while ((abs(M4-cmu(4)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(3)=c(3)-(M4-cmu(4))/dc3;
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
tolerance=.0000001*cmu(5);
mm=0;
while ((abs(M5-cmu(5)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(4)=c(4)-(M5-cmu(5))/dc4;
[M5,dc4] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,5,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
tolerance=.0000001*cmu(6);
mm=0;
while ((abs(M6-cmu(6)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(5)=c(5)-(M6-cmu(6))/dc5;
[M6,dc5] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,6,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M7,dc6] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,7,EZ);
tolerance=.0000001*cmu(7);
mm=0;
while ((abs(M7-cmu(7)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(6)=c(6)-(M7-cmu(7))/dc6;
[M7,dc6] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,7,EZ);
end
c0=-(c(2)+3*c(4)+15*c(6));
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[M8,dc7] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,8,EZ);
tolerance=.0000001*cmu(8);
mm=0;
while ((abs(M8-cmu(8)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(7)=c(7)-(M8-cmu(8))/dc7;
[M8,dc7] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,8,EZ);
end
[SecondMoment] = CalculateSecondMomentC7(c0,c);
c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end
end
end
.
.
Code: Select all
function [SecondMoment] = CalculateSecondMomentC7(a0,a)
SecondMoment=a0^2+a(1)^2 +2* a0* a(2) +3*(a(2).^2+2* a(1).* a(3) +2* a0* a(4))+ ...
15*(a(3).^2 +2 *a(2).* a(4) +2 *a(1).* a(5)+2* a0* a(6))+105*(a(4).^2 +2* a(3).* a(5)+2* a(2).* a(6) +2* a(1).* a(7) )+ ...
945*(a(5).^2 +2* a(4).* a(6) +2 *a(3).* a(7) )+(10395)*(a(6).^2+2* a(5).* a(7) )+(135135)*a(7).^2;
end
.
.
Code: Select all
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
.
.
Code: Select all
function [F,dF] = CalculateMomentsAndDerivatives_0(Ms,a0,a,SeriesOrder,NZterms,NMoments)
%[F,dF] = CalculateCumulantsAndDerivativesFromMoments(C,a0,a,SeriesOrder,SeriesOrder,NoOfCumulants);
%a(1)=a1;
%a(2)=a2;
%a(3)=a3;
%a(4)=a4;
%a(5)=a5;
%a(6)=a6;
%a(7)=a7;
%aa0=a0;
if(NMoments>8)
a0=-(a(2)+3*a(4)+15*a(6)+105*a(8));% ---1
end
if(NMoments<=8)
a0=-(a(2)+3*a(4)+15*a(6));% ---1
end
EZ(1)=0;
EZ(2)=1;
for nn=3:NMoments*SeriesOrder+NZterms+2
if rem(nn,2)==1
EZ(nn)=0;
else
EZ(nn)=EZ(nn-2)*(nn-1);
EZ(nn);
end
end
EZ;
EXZ(1,1)=1;
for pp1=1:NZterms
EXZ(1,pp1+1)=EZ(pp1);
end
a(SeriesOrder+1:NMoments*SeriesOrder+1)=0;
b0=a0;
b=a;
for mm=1:NMoments
if(mm>1)
[b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*mm);
b(SeriesOrder*mm+1:NMoments*SeriesOrder+1)=0;
end
% b0
% b
%str=input('Look at numbers')
EXZ(mm+1,1)=b0;
for pp2=1:SeriesOrder*mm
EXZ(mm+1,1)=EXZ(mm+1,1)+b(pp2).*EZ(pp2);
end
for pp1=1:NZterms
EXZ(mm+1,pp1+1)=b0.*EZ(pp1);
for pp2=1:SeriesOrder*mm
EXZ(mm+1,pp1+1)=EXZ(mm+1,pp1+1)+b(pp2).*EZ(pp2+pp1);
end
end
end
%u1=EXZ(2,1);
u2=EXZ(3,1);
u3=EXZ(4,1);
u4=EXZ(5,1);
u1=a0+a(2);
if(SeriesOrder>=4)
u1=a0+a(2)+3*a(4);
end
if(SeriesOrder>=6)
u1=u1+15*a(6);
end
if(SeriesOrder>=8)
u1=u1+105*a(8);
end
%k2=u2-u1^2;
%k3=u3-3*u2*u1+2*u1^3;
%k4=u4-4*u3*u1-3*u2^2+12*u2*u1^2-6*u1^4;
du1(1)=1;%----2
% du1(2)=0;%----2
% du1(3)=1;%----2
% du1(4)=0;%----2
% du1(5)=1;%----2
% du1(6)=0;%----2
% du1(7)=15;%----2
% du1(8)=0;%----2
% du1(9)=105;%----2
% du1(10)=0;%----2
du2(1)=2*EXZ(2,1);
du3(1)=3*EXZ(3,1);
du4(1)=4*EXZ(4,1);
%du1(1)=1;%----2
%du2(1)=0;
%du3(1)=0;
%du4(1)=0;
for mm=2:SeriesOrder+1
du1(mm)=EXZ(1,mm);
du2(mm)=2*EXZ(2,mm);
du3(mm)=3*EXZ(3,mm);
du4(mm)=4*EXZ(4,mm);
end
if(NMoments>=5)
u5=EXZ(6,1);
%du5(1)=5*EXZ(5,1);
for mm=1:SeriesOrder+1
du5(mm)=5*EXZ(5,mm);
end
end
if(NMoments>=6)
u6=EXZ(7,1);
%du6(1)=0;
for mm=1:SeriesOrder+1
du6(mm)=6*EXZ(6,mm);
end
end
if(NMoments>=7)
u7=EXZ(8,1);
%str=input('Look at k7 and k71')
%du7(1)=0;
for mm=1:SeriesOrder+1
du7(mm)=7*EXZ(7,mm);
end
end
if(NMoments>=8)
u8=EXZ(9,1);
for mm=1:SeriesOrder+1
du8(mm)=8*EXZ(8,mm);
end
end
if(NMoments>=9)
u9=EXZ(10,1);
for mm=1:SeriesOrder+1
du9(mm)=9*EXZ(9,mm);
end
end
if(NMoments>=10)
u10=EXZ(11,1);
for mm=1:SeriesOrder+1
du10(mm)=10*EXZ(10,mm);
end
end
F(1,1)=u1-Ms(1);
F(2,1)=u2-Ms(2);
F(3,1)=u3-Ms(3);
F(4,1)=u4-Ms(4);
if(NMoments>=5)
F(5,1)=u5-Ms(5);
end
if(NMoments>=6)
F(6,1)=u6-Ms(6);
end
if(NMoments>=7)
F(7,1)=u7-Ms(7);
end
if(NMoments>=8)
F(8,1)=u8-Ms(8);
end
if(NMoments>=9)
F(9,1)=u9-Ms(9);
end
if(NMoments>=10)
F(10,1)=u10-Ms(10);
end
for mm=1:SeriesOrder+1
dF(1,mm)=du1(mm);
dF(2,mm)=du2(mm);
dF(3,mm)=du3(mm);
dF(4,mm)=du4(mm);
if(NMoments>=5)
dF(5,mm)=du5(mm);
end
if(NMoments>=6)
dF(6,mm)=du6(mm);
end
if(NMoments>=7)
dF(7,mm)=du7(mm);
end
if(NMoments>=8)
dF(8,mm)=du8(mm);
end
if(NMoments>=9)
dF(9,mm)=du9(mm);
end
if(NMoments>=10)
dF(10,mm)=du10(mm);
end
end
end
.
.
Code: Select all
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
.
.
Code: Select all
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
.
.
Code: Select all
function [IsValidFlag] = CheckIsValidDensity(cc0,cc)
%Roughly checks if the density is valid. It is not a definitive check but
%just a rough one but mostly works.
Z=2.0;
Xp1=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;
Z=3.0;
Xp2=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;
Z=4.65;
Xp3=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;
Z=-2.0;
Xn1=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;
Z=-3.0;
Xn2=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;
Z=-4.65;
Xn3=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;
IsValidFlag=0;
if( ((Xp2>Xp1) && (Xp3>Xp2)) &&(Xp1>Xn1) && ((Xn2<Xn1) && (Xn3<Xn2)))
IsValidFlag=1;
end
end
.
.
When you run the program, you will see a comparison of input moments (rMu) and fitted moments.
rMu =
1.0e+05 *
Columns 1 through 7
0.000012499863483 0.000021197093700 0.000050274347857 0.000171368317209 0.000856094970826 0.006320436615909 0.068698344472325
Column 8
1.084082326087001
FittedMoments =
1.0e+05 *
Columns 1 through 7
0.000012499863483 0.000021197093700 0.000050274347857 0.000171368317209 0.000856094970826 0.006320436615909 0.068698344472325
Column 8
1.084082326087000
Look at comparison of fitted moments with input moments.
And you will see the following graph of the density.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, I want to describe my numerical strategy/algorithm to fit the Z-series coefficients to prescribed moments.
Following are the main points
1. The whole procedure is iterative.
2. In our standardized setup, the first moment is always zero. First moment of a Z-series density is given as [$]c_0 + c_2 + 3 c_4 + 15 c_6[$] which gives us the zero mean condition on [$]c_0[$] as [$]c_0 = -( c_2 + 3 c_4 + 15 c_6) [$], so whenever [$]c_2[$] or [$]c_4[$] or [$]c_6[$] changes we again adjust [$]c_0[$] to make sure that zero mean condition is satisfied. You will see the above equation again and again in the iterative algorithm.
3. In our standardized setup, second standard moment is always one. Whenever any coefficient changes, it changes the effective second moment and we have to divide all coefficients with square root of (new) second moment which ensures that second moment goes to one after the division with (altered new second moment after changing any coefficient).
By using step 2 and 3 above after changing any coefficient, we ensure that first and second moments remain standardized to zero and one.
4. The most important clue I got when I was playing with moments was that coefficient c_1 is associated with second moment, c_2 is associated with third moment, c_3 is associated with fourth moment and so on. If you calibrate moments to an even (around median) density with zero skew, you will notice that there will only be non-zero values for c_1, c_3 , c_5, c_7 and c_2, c_4, c_6 will be exactly zero. So over time I gained some intuition that each of the moment has a major coefficient of one smaller degree associated with it. Though other coefficients affect each moment, the principal effect after accounting for effect from smaller coefficients is from the major coefficient.
5. Main algorithm is very simple. I start with c_0=0 and c_1 =1, and then I add third moment to calibration. I use Newton-Raphson to alter c_2 until third moment is perfectly calibrated. Then I go back and make sure using step 2 by altering c_0 that first moment goes to zero and ,step 3, divide by square-root of new second moment to ensure that second moment is standardized again. Obviously when I alter c_0 and divide by square-root of first moment, the 3rd moment that we calibrated with Newton-Raphson will slightly go out of calibration with third moment so we repeat the whole process with Newton-Raphson again and then standardize again.
Once first three moments are calibrated, I alter c_3 according to Newton-Raphson so that fourth moment is perfectly calibrated. But perfect calibration here means that first three moments would be slightly out again so calibration is repeated on them to bring them in line again.
Once first four moments have been calibrated, I change c_4 using Newton-Raphson so that fifth moment is perfectly calibrated. And then I repeat on all previous moments so that they remain in calibration.
So we continue till eight moments using the same procedure.
After a few iterations of the algorithm, the coefficients are firmly in the catchment of multi-dimensional Newton method that finds exact calibration to moments in 2-3 iterations after that.
The algorithm works because once we calibrate a moment with Newton-Raphson using associated major coefficient, the newly created mismatch in smaller moments is very little and decreases fast with every iteration of the algorithm.
Though I tried the algorithm on only first eight moments, you can safely use the algorithm for something like 12 moments or higher if you wish.
Once you have the appropriate intuition, it is a very simple algorithm.
Following are the main points
1. The whole procedure is iterative.
2. In our standardized setup, the first moment is always zero. First moment of a Z-series density is given as [$]c_0 + c_2 + 3 c_4 + 15 c_6[$] which gives us the zero mean condition on [$]c_0[$] as [$]c_0 = -( c_2 + 3 c_4 + 15 c_6) [$], so whenever [$]c_2[$] or [$]c_4[$] or [$]c_6[$] changes we again adjust [$]c_0[$] to make sure that zero mean condition is satisfied. You will see the above equation again and again in the iterative algorithm.
3. In our standardized setup, second standard moment is always one. Whenever any coefficient changes, it changes the effective second moment and we have to divide all coefficients with square root of (new) second moment which ensures that second moment goes to one after the division with (altered new second moment after changing any coefficient).
By using step 2 and 3 above after changing any coefficient, we ensure that first and second moments remain standardized to zero and one.
4. The most important clue I got when I was playing with moments was that coefficient c_1 is associated with second moment, c_2 is associated with third moment, c_3 is associated with fourth moment and so on. If you calibrate moments to an even (around median) density with zero skew, you will notice that there will only be non-zero values for c_1, c_3 , c_5, c_7 and c_2, c_4, c_6 will be exactly zero. So over time I gained some intuition that each of the moment has a major coefficient of one smaller degree associated with it. Though other coefficients affect each moment, the principal effect after accounting for effect from smaller coefficients is from the major coefficient.
5. Main algorithm is very simple. I start with c_0=0 and c_1 =1, and then I add third moment to calibration. I use Newton-Raphson to alter c_2 until third moment is perfectly calibrated. Then I go back and make sure using step 2 by altering c_0 that first moment goes to zero and ,step 3, divide by square-root of new second moment to ensure that second moment is standardized again. Obviously when I alter c_0 and divide by square-root of first moment, the 3rd moment that we calibrated with Newton-Raphson will slightly go out of calibration with third moment so we repeat the whole process with Newton-Raphson again and then standardize again.
Once first three moments are calibrated, I alter c_3 according to Newton-Raphson so that fourth moment is perfectly calibrated. But perfect calibration here means that first three moments would be slightly out again so calibration is repeated on them to bring them in line again.
Once first four moments have been calibrated, I change c_4 using Newton-Raphson so that fifth moment is perfectly calibrated. And then I repeat on all previous moments so that they remain in calibration.
So we continue till eight moments using the same procedure.
After a few iterations of the algorithm, the coefficients are firmly in the catchment of multi-dimensional Newton method that finds exact calibration to moments in 2-3 iterations after that.
The algorithm works because once we calibrate a moment with Newton-Raphson using associated major coefficient, the newly created mismatch in smaller moments is very little and decreases fast with every iteration of the algorithm.
Though I tried the algorithm on only first eight moments, you can safely use the algorithm for something like 12 moments or higher if you wish.
Once you have the appropriate intuition, it is a very simple algorithm.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, in 2-3 days, I will post 5-7 matlab programs for dealing with Z-series densities. Some of the programs I have posted before.
These programs will include
1. Calculation of probability density function, cumulative density function and associated value of standard normal.
2. Conversion from Z-series to Hermite series and back.
3. Calculation of variances using hermite polynomial series.
4. Addition of independent densities using squared coefficients of hermite polynomials. You do not need convolutions or transforms for this anymore. This can make large changes for non-linear filtering and will make non-linear filtering relatively simpler.
5. Calculation of Z-series representation of functions of the original Z-series random variable.
These programs will include
1. Calculation of probability density function, cumulative density function and associated value of standard normal.
2. Conversion from Z-series to Hermite series and back.
3. Calculation of variances using hermite polynomial series.
4. Addition of independent densities using squared coefficients of hermite polynomials. You do not need convolutions or transforms for this anymore. This can make large changes for non-linear filtering and will make non-linear filtering relatively simpler.
5. Calculation of Z-series representation of functions of the original Z-series random variable.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Happy new year 2023 to all friends. I wish that this new year brings happiness, prosperity and success to all of us. And I hope to present better research on my thread for friends this year than I have been able to do in past few years.
Friends, my research work on Z-series densities has become quite mature. We can almost do anything with Z-series densities that we can do with other densities and this presentation is probably most general and it encompasses most other densities. What is very interesting that we can find the Z-series density once we are given first few moments of the random variable using my new algorithm. I also presented a method based on derivatives of densities using which any known density can be converted into Z-series and presented an example program that finds Z-series density representation from any lognormal density.
Here I presented the program to convert arbitrary lognormal density into Z-series representation which can easily be generalized to any other density with continuous derivatives to find its Z-series representation: https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1635#p873551 and
https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1635#p873552
Since all of this is very interesting development for stochastics, probability theory and statistics, I want to write a paper on Z-series densities, general formulas related to Z-series densities, and their derivation from known distributions, and also their derivation from first eight moments of the random variable. I hope that first version of the paper would be ready in 7-10 days. This is totally new research but open on internet so I want to request friends to make sure nobody else could present the ideas as genuinely their own in some journal quickly even though most probably my fears are completely baseless. Of course all friends are free to write their research based on top of ideas presented in this thread and get it published but I hope they will acknowledge the basic work I have done.
Anyway, I am going to write a formal paper and present it here in 10-15 days.
Friends, my research work on Z-series densities has become quite mature. We can almost do anything with Z-series densities that we can do with other densities and this presentation is probably most general and it encompasses most other densities. What is very interesting that we can find the Z-series density once we are given first few moments of the random variable using my new algorithm. I also presented a method based on derivatives of densities using which any known density can be converted into Z-series and presented an example program that finds Z-series density representation from any lognormal density.
Here I presented the program to convert arbitrary lognormal density into Z-series representation which can easily be generalized to any other density with continuous derivatives to find its Z-series representation: https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1635#p873551 and
https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1635#p873552
Since all of this is very interesting development for stochastics, probability theory and statistics, I want to write a paper on Z-series densities, general formulas related to Z-series densities, and their derivation from known distributions, and also their derivation from first eight moments of the random variable. I hope that first version of the paper would be ready in 7-10 days. This is totally new research but open on internet so I want to request friends to make sure nobody else could present the ideas as genuinely their own in some journal quickly even though most probably my fears are completely baseless. Of course all friends are free to write their research based on top of ideas presented in this thread and get it published but I hope they will acknowledge the basic work I have done.
Anyway, I am going to write a formal paper and present it here in 10-15 days.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends a bit after writing the above post I went out for a walk as my only new year resolution is to lose weight. I had walked one and a half kilometer and I was walking on an almost empty service road. There were a lot of people and also quite some traffic on the main road. Two people in black jackets who had covered their faces with black woven face-masks (not the covid medical ones but professional face-masks) approached me from behind on a motorcycle. They stopped me and started asking why I was walking on the service road. Though they were not in uniform, they claimed they were police. They started threatening me that I should not walk on the road and I must walk only on a track in the park otherwise they would arrest me in the future. They continued to harass me for quite a bit until several people around noticed and approached us and then both men left on the motorcycle.
I have walked on the same path several thousand times since I have been living in the present house in 2005. Some local police staff even knows me and greets me when they sometime come across me in very early morning since I have been doing that for more than a decade and they well recognize me.
I want to tell friends that people behind my persecution are becoming impatient and restless with all of this when things are not going their way and they might ask Pakistan army to use force on me. Pakistan army lifts people who they do not like on a daily basis, tortures them for a few days in confinement and then releases them on random streets and then army feigns ignorance about all of this and it is said that unknown people lifted the victim. Everybody knows that army lifted the victims but nobody mentions army and the word of "unknown people" is used about perpetrators of such incident. Police and courts are never able to help any such victims.
When I started telling people about mind control nobody believed me and only later most of them realized the truth about what I was trying to tell them. I want to tell friends that I would have been killed/shot in New York 25 years ago by some planted mugger (with complete planning) or would have been given poison when my talent was discovered and they failed to retard me since it was impossible for many people to stand or accept a Muslim with special neurotransmitters. I remained alive only because a large number of people in the university had already come to know about it and if I were killed, most of them would have known who was behind it. I bitterly wanted to remain in United States and work there but once I realized that most people wanted to retard me with an iron fist, I knew I had no future there and I decided to return to Pakistan and I know I made the right decision.
Later I remained alive since my family was assured when they cooperated that I would not be killed or made "special" and again a large number of people in my broader family, army and several influential people in Pakistan got to know about it. But the way people behind my persecution are becoming impatient again, I am seriously afraid they might do what they never did 25 years ago in New York. I am very afraid they might give money to some mugger to shoot me or create some other way to end my life with full cooperation and planning of Pakistan army. I want to tell friends very clearly that if I had premature seemingly accidental death, it would never be accidental, it would be a completely pre-planned murder.
I have walked on the same path several thousand times since I have been living in the present house in 2005. Some local police staff even knows me and greets me when they sometime come across me in very early morning since I have been doing that for more than a decade and they well recognize me.
I want to tell friends that people behind my persecution are becoming impatient and restless with all of this when things are not going their way and they might ask Pakistan army to use force on me. Pakistan army lifts people who they do not like on a daily basis, tortures them for a few days in confinement and then releases them on random streets and then army feigns ignorance about all of this and it is said that unknown people lifted the victim. Everybody knows that army lifted the victims but nobody mentions army and the word of "unknown people" is used about perpetrators of such incident. Police and courts are never able to help any such victims.
When I started telling people about mind control nobody believed me and only later most of them realized the truth about what I was trying to tell them. I want to tell friends that I would have been killed/shot in New York 25 years ago by some planted mugger (with complete planning) or would have been given poison when my talent was discovered and they failed to retard me since it was impossible for many people to stand or accept a Muslim with special neurotransmitters. I remained alive only because a large number of people in the university had already come to know about it and if I were killed, most of them would have known who was behind it. I bitterly wanted to remain in United States and work there but once I realized that most people wanted to retard me with an iron fist, I knew I had no future there and I decided to return to Pakistan and I know I made the right decision.
Later I remained alive since my family was assured when they cooperated that I would not be killed or made "special" and again a large number of people in my broader family, army and several influential people in Pakistan got to know about it. But the way people behind my persecution are becoming impatient again, I am seriously afraid they might do what they never did 25 years ago in New York. I am very afraid they might give money to some mugger to shoot me or create some other way to end my life with full cooperation and planning of Pakistan army. I want to tell friends very clearly that if I had premature seemingly accidental death, it would never be accidental, it would be a completely pre-planned murder.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Though Z-series densities can be extremely helpful in inference and in so many other ways, they can also be extremely helpful for simulation of random variables of difficult densities where only imperfect and difficult methods for simulation exist. All you have to do is to be able to find derivatives of density of random variable and use them to find equivalent Z-series random variable representation like I have done for lognormal density as in post 1636, and 1637. And then simulation of Z-series random variables only requires simulation of standard gaussian random variable which is extremely well-understood.
I want to add several examples of conversion of random variables of different densities to random variables of Z-series in my new upcoming paper.
I want to add several examples of conversion of random variables of different densities to random variables of Z-series in my new upcoming paper.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
As friends are aware that an antipsychotic injection is forced on me due to machinations of mind control agency and people backing it. American army gives money to Pakistan army generals (several hundred million dollars every year) who among other things ask Pakistani psychiatrists to declare me schizophrenic despite that I am perfectly mentally sane and in better mental shape than most of the psychiatrists themselves. psychiatrist following the orders of Pak army will strictly detain me for several months if I do not take the injection. Most psychiatrists have a private army of about a dozen poorly educated and low paid young men and there is little oversight about actions of the psychiatrists.
Anyway an antipsychotic injection is due in another two days. Most Fluanxol injections by lundbeck distributed in Lahore city already have mind control drugs added at manufacturing source. And the effect of mind control chemicals added to injections is so severe that it is very difficult to work for first ten days after the injection and I remain in an extremely poor state of mind with far decreased consciousness.
I remained better for past three weeks since when I went to Kot Addu three weeks ago, the first thing I did in the city was to buy antipsychotic injection before they would get the idea of drugging it and apparently the supply of injections in Southern Punjab probably did not have mind control chemicals in it.
If it were not for these injections, mind control would not become effective on me. Mind control agencies have to ask psychiatrists to declare me mentally ill so that I have to take some medicine/injections and then every medicine or injections is distributed in the entire market with mind control chemicals to control me.
Despite that injections distributed in Lahore city are already drugged with mind control chemicals, I am very afraid that mind control agencies might try to add even more mind control chemicals in special injections to control me. I really want to request mainland European embassies to please keep a vigil on any such activity of Pakistan army and try to stop them from any such thing. I really want to request good European people to urge their embassies to force the military and influential people in Pakistan to end my mind control. I would have asked good American people to do the same but American embassy is not controlled by American people, it is controlled by backers of industrial-military complex where people are dying to retard me forever.
Now several people who do not like me are waiting for my getting injected.
Anyway an antipsychotic injection is due in another two days. Most Fluanxol injections by lundbeck distributed in Lahore city already have mind control drugs added at manufacturing source. And the effect of mind control chemicals added to injections is so severe that it is very difficult to work for first ten days after the injection and I remain in an extremely poor state of mind with far decreased consciousness.
I remained better for past three weeks since when I went to Kot Addu three weeks ago, the first thing I did in the city was to buy antipsychotic injection before they would get the idea of drugging it and apparently the supply of injections in Southern Punjab probably did not have mind control chemicals in it.
If it were not for these injections, mind control would not become effective on me. Mind control agencies have to ask psychiatrists to declare me mentally ill so that I have to take some medicine/injections and then every medicine or injections is distributed in the entire market with mind control chemicals to control me.
Despite that injections distributed in Lahore city are already drugged with mind control chemicals, I am very afraid that mind control agencies might try to add even more mind control chemicals in special injections to control me. I really want to request mainland European embassies to please keep a vigil on any such activity of Pakistan army and try to stop them from any such thing. I really want to request good European people to urge their embassies to force the military and influential people in Pakistan to end my mind control. I would have asked good American people to do the same but American embassy is not controlled by American people, it is controlled by backers of industrial-military complex where people are dying to retard me forever.
Now several people who do not like me are waiting for my getting injected.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
A very annoying practice of US army is that among many other high-handed practices, they stop emails and linkedin messages of people they do not like. In past ten years there were numerous instances when they stopped my emails. Just one incident I want to mention here is that I sent emails to each and every member of upper and lower houses of UK when I was there in 2010. I tried to tell them about my persecution and how brazenly London city was getting drugged. Aside from Britvic, there are four or five large departmental chains in UK and all of them were manufacturing beverages, juices and bottled water that was thoroughly drugged with mind control chemicals. And every single one of the emails I sent to members of both houses of parliament I sent was stopped. I went to Parliament house in London in person but I was not allowed to enter inside despite that it was my legal right. They simply asked me that I should come after a week as it was very cold. They were allowing everyone else to go inside. And they were becoming belligerent when I protested and I never wanted a fight so I had to return.
There have been many many other instances when I tried to approach people in other countries and my emails or linkedin messages were stopped.
There have been many many other instances when I tried to approach people in other countries and my emails or linkedin messages were stopped.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, I started writing my paper this evening and wrote first three pages. I hope to quickly complete the paper in next 5-7 days. Here is the image to first page while you can download the initial three pages from download file link below.
.
.
.
Here is the file containing initial three pages.
.
.
.
Here is the file containing initial three pages.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, in order to show how well the Z-series density construction from input moments, I have created a few examples.
I used three parameters of generalized gamma distribution to calculate its first eight raw moments. I fed these moments as input to my Z-series density construction program posted above in post 1653 above on the same page. Then I compare the Z-series density from fitted moments to original generalized gamma density from the same three parameters that were used to calculate the moments. And then graph both Z-series density and generalized gamma density on top of each other.
Here are a few graphs from the above procedure.((I will also post matlab program to create these comparison graphs)
The following graphs show three parameters generalized Gamma distribution(as in wikipedia; https://en.wikipedia.org/wiki/Generaliz ... stribution)
as input and we calculates its eight raw moments. These first eight raw moments are input to our Z-series density moment fitting program (post 1683) above and you can see that match from our Z-Series fitted density to generalized gamma distribution is mostly exact.
Values of Parameters of the original Generalized gamma distribution a, d, and p are given on title of the graph.
.
.
.
.
Matlab program to generate above graphs is given below.
.
.
.
.
I will like to mention here that sometimes density fitted with Z-series random variables is not a good fit. Mostly it can have two reasons.
1. Some more moments are needed to perfectly fit the Z-series density and first eight moments might not be enough to pin the density.
2. Sometime fitted density of Z-series variable goes into origin(zero) and then folds back onto itself and becomes positive again. In such cases, we would have to add both branches of the density to get one single estimate of the density. I have still not done that in my experimental programs and would do that in next few days. Please keep checking the thread. After adding both branches of the density we would get a perfect estimate of the true density.
I really want to request friends to not confuse this new method with Gram-Charlier or edgeworth expansions. It fits perfectly to input moments and gives a shape of distribution that is very close to original shape when the true shape is known as in above examples.
I used three parameters of generalized gamma distribution to calculate its first eight raw moments. I fed these moments as input to my Z-series density construction program posted above in post 1653 above on the same page. Then I compare the Z-series density from fitted moments to original generalized gamma density from the same three parameters that were used to calculate the moments. And then graph both Z-series density and generalized gamma density on top of each other.
Here are a few graphs from the above procedure.((I will also post matlab program to create these comparison graphs)
The following graphs show three parameters generalized Gamma distribution(as in wikipedia; https://en.wikipedia.org/wiki/Generaliz ... stribution)
as input and we calculates its eight raw moments. These first eight raw moments are input to our Z-series density moment fitting program (post 1683) above and you can see that match from our Z-Series fitted density to generalized gamma distribution is mostly exact.
Values of Parameters of the original Generalized gamma distribution a, d, and p are given on title of the graph.
.
.
.
.
Matlab program to generate above graphs is given below.
.
Code: Select all
function [] = FitMomentsOfGeneralizedGammaRV(a,d,p)
%The following program takes three parameters of
%generalized Gamma distribution(as in wikipedia; https://en.wikipedia.org/wiki/Generalized_gamma_distribution)
%as input and calculates its raw moments. These first eight raw moments are
%input to our Z-series density moment fitting program (post 1683) at web
%address below
%https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1650#p873803
%and you can see that
%match to generalized gamma distribution is mostly exact.
%Below calculate first eight raw moments of generalized gamma distribution.
rMu(1:8)=0;
for nn=1:8
rMu(nn)=a^nn*gamma((d+nn)/p)/gamma(d/p);
end
%Feed those raw moments to our Z-series random variable estimation program
[c0,c,FittedMoments] = CalculateZSeriesDensityFromRawMoments02(rMu);
%Below construct generalized gamma density on Gaussian grid.
%First calculate normal random variable on a grid below
dNn=.1/2; % Normal density subdivisions width. would change with number of subdivisions
Nn=45*4; % No of normal density subdivisions
NnMid=((1+Nn)/2)*dNn;
Z(1:Nn)=(((1:Nn)*dNn)-NnMid);
%Z
%str=input('Look at Z');
%Now calculate Y as a function of normal random variable.
Y(1:Nn)=c0;
for nn=1:7
Y(1:Nn)=Y(1:Nn)+c(nn)*Z(1:Nn).^nn;
end
%Now take change of densities derivative of Y with respect to normal
DfY(1:Nn)=0;
for nn=2:Nn-1
DfY(nn) = (Y(nn + 1) - Y(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
DfY(Nn)=DfY(Nn-1);
DfY(1)=DfY(2);
%Now calculate the density of Y from density of normal random variable
%using change of probability derivative.
pY(1:Nn)=0;
for nn = 1:Nn
pY(nn) = (normpdf(Z(nn),0, 1))/abs(DfY(nn));
end
%Below calculate generalized gamma distribution from defintion of its
%probability density given the parameters a, d and p.
X(1)=0;
X(5000+1)=Y(Nn);
Delta=(Y(Nn))/5000;
X(2:5001)=X(1)+(1:5000)*Delta;
Xpdf = p/a^d/gamma(d/p).*X.^(d-1).*exp(-(X./a).^p);
plot(Y(1:Nn),pY(1:Nn),'b',X(1:5001),Xpdf(1:5001),'g');
title(sprintf('Comparison of Z-Series RV Density with Fitted Moments and Original Genralized Gamma Density, a = %.3f,d=%.3f,p=%.3f', a,d,p));
%title('Comparison of Z-Series RV Density with Fitted Moments and Original Genralized Gamma Density');
legend({'Z-series Random Variable Density','Generalized Gamma Density'},'Location','northeast')
str=input('Look at comparison between Z-series Fitted Density and original Generalized Gamma Density');
end.
.
I will like to mention here that sometimes density fitted with Z-series random variables is not a good fit. Mostly it can have two reasons.
1. Some more moments are needed to perfectly fit the Z-series density and first eight moments might not be enough to pin the density.
2. Sometime fitted density of Z-series variable goes into origin(zero) and then folds back onto itself and becomes positive again. In such cases, we would have to add both branches of the density to get one single estimate of the density. I have still not done that in my experimental programs and would do that in next few days. Please keep checking the thread. After adding both branches of the density we would get a perfect estimate of the true density.
I really want to request friends to not confuse this new method with Gram-Charlier or edgeworth expansions. It fits perfectly to input moments and gives a shape of distribution that is very close to original shape when the true shape is known as in above examples.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, Pakistan army is drugging Lahore city at a very fast pace again on behest of CIA. Earlier they were confined to neighborhoods close to Johar Town where I live but not they are drugging food and beverages in even remote areas. This has been going on for past 5-6 days. I really thought that people will take notice of it and stop the army but they are ruthlessly drugging the whole city. They go into markets and drug the water, food and beverages on the shelves in entire markets. People who are aware of tactics of CIA know what I am talking about.
I want to request good American friends to pressurize their government to stop ruthless drugging of Lahore city on behest of CIA to corner me. I also believe that CIA, and Pak army have connived with my family if I could be targeted with special drugs, then I might lose my senses and my father would restart my treatment saying that I lost my mental balance and he had to act to help me in the extreme case. CIA and Pak Army have are machinating again to detain me with some psychiatrist and restart injections and antipsychotic medicine filled with mind control drugs.
I really want to request good American friends to intervene now before it is too late.
I also want to request Europeans to ask their embassies to keep a vigil and continue to monitor the tactics of CIA as they drug Lahore city. Please ask them to stay in touch with the matter using their sources to know the reality of mind control as it is also used on a lot of talented Europeans. This way European nations would also have true knowledge about the nature of mind control.
I have been able to hold for past few days since I have been very careful about food. I still got hit with bad food several times but was able to recover and remain perfectly sane since there is still some good food available on stores if one is careful.
I again want to request American friends to help me and force CIA to stop drugging Lahore city.
I want to request good American friends to pressurize their government to stop ruthless drugging of Lahore city on behest of CIA to corner me. I also believe that CIA, and Pak army have connived with my family if I could be targeted with special drugs, then I might lose my senses and my father would restart my treatment saying that I lost my mental balance and he had to act to help me in the extreme case. CIA and Pak Army have are machinating again to detain me with some psychiatrist and restart injections and antipsychotic medicine filled with mind control drugs.
I really want to request good American friends to intervene now before it is too late.
I also want to request Europeans to ask their embassies to keep a vigil and continue to monitor the tactics of CIA as they drug Lahore city. Please ask them to stay in touch with the matter using their sources to know the reality of mind control as it is also used on a lot of talented Europeans. This way European nations would also have true knowledge about the nature of mind control.
I have been able to hold for past few days since I have been very careful about food. I still got hit with bad food several times but was able to recover and remain perfectly sane since there is still some good food available on stores if one is careful.
I again want to request American friends to help me and force CIA to stop drugging Lahore city.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, my work is going very slowly and on most of the past days I could not work. Also because I remained down after taking drugged food. But here is the next incremental version of the Z-series paper. I will really try to complete the first version by the end of coming weekend if things go well. Here is the link to current incremental version in which I have described how to find Z-series of function transformation of a random variable if original Z-series of the random variable is known. I have also described how to convert from Z-series coordinates to hermite polynomial coordinates. I have also described how to add two independent Z-series random variables through squared addition of hermite polynomial representation coefficients.
I want to mention that function transformations of the Z-series random variables are done in original Z-series coordinates while analytics like addition of independent variables and variances is done in Hermite polynomial coordinates.
Here is the new incremental version of the paper. [attachment=0]Z-SeriesIntroductionA1.pdf[/attachment]
I will really try to complete the first version by the end of coming weekend if things go well and I am not hit hard by mind control drugs being spread all over the city by military secret services.
I want to mention that function transformations of the Z-series random variables are done in original Z-series coordinates while analytics like addition of independent variables and variances is done in Hermite polynomial coordinates.
Here is the new incremental version of the paper. [attachment=0]Z-SeriesIntroductionA1.pdf[/attachment]
I will really try to complete the first version by the end of coming weekend if things go well and I am not hit hard by mind control drugs being spread all over the city by military secret services.
- Attachments
-
- Z-SeriesIntroductionA1.pdf
- (176.62 KiB) Downloaded 35 times
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, another new idea that I want to include in my new paper is about non-linear regressions with Z-Series Hermite polynomials. The standard linear regression that we currently use is as if we are regressing on data represented by first hermite polynomial. But with Z-Series, we can extend the idea of linear regression on first hermite to non-linear regressions on multiple hermite polynomials.
Suppose we want to do a nonlinear regression on W from two explanatory variables X and Y
I will give brief details to outline the new method in a few points.
First we would have to find Z-series representation of all the variables involved including regressand and regressors. Usually this would be done by calculating moments of all the data and fitting a Z-series density on each of the regressand and regressors.
For example in our case we would have.
First we find Z-Series of regressand
[$]W=a_0 + a_1 \, Z_1 + a_2 \, Z_1^2 + a_3 \, Z_1^3 + a_4 \, Z_1^4 + ... +a_7 \, Z_1^7[$]
Then we find Z-series of regressors
[$]X=b_0 + b_1 \, Z_2 + b_2 \, Z_2^2 + b_3 \, Z_2^3 + b_4 \, Z_2^4 + ... +b_7 \, Z_2^7[$]
[$]Y=c_0 + c_1 \, Z=3 + c_2 \, Z_3^2 + c_3 \, Z_3^3 + c_4 \, Z_3^4 + ... +c_7 \, Z_3^7[$]
Then we convert from Z-Series representation to hermite representation as I have done many times on this forum and also mentioned in last incomplete version of my paper above as
[$]W=ah_0 + ah_1 \, H_1(Z_1) + ah_2 \, H_2(Z_1) + ah_3 \, H_3(Z_1) + ah_4 \, H_4(Z_1) + ... +ah_7 \, H_7(Z_1)[$]
and
[$]X=bh_0 + bh_1 \, H_1(Z_2) + bh_2 \, H_2(Z_2) + bh_3 \, H_3(Z_2) + bh_4 \, H_4(Z_2) + ... +bh_7 \, H_7(Z_2)[$]
[$]Y=ch_0 + ch_1 \, H_1(Z_3) + ch_2 \, H_2(Z_3) + ch_3 \, H_3(Z_3) + ch_4 \, H_4(Z_3) + ... +ch_7 \, H_7(Z_3)[$]
The purpose of converting to hermites is that we can convert our problem of regression into several sub-problems since hermite polynomials are orthogonal to each other.
We want to regress each hermite polynomial term of the regressand on same order hermite polynomial terms of the regressors
Now we regress 1st hermite term of regressand on first hermite terms of regressors as
Here [$]k_{11}[$] and [$]k_{12}[$] are regression coefficients that we will find through ordinary least squares algorithm.
[$]ah_1 \, H_1(Z_1) =k_{11} (bh_1 \, H_1(Z_2) )+k_{12} (ch_1 \, H_1(Z_3))[$]
and for second hermite polynomial regress as
Here [$]k_{21}[$] and [$]k_{22}[$] are regression coefficients that we will find through ordinary least squares algorithm.
[$]ah_2 \, H_2(Z_1) =k_{21} (bh_2 \, H_2(Z_2) )+k_{22} (ch_2 \, H_2(Z_3))[$]
and for third polynomial ,we will regress as
[$]ah_3 \, H_3(Z_1) =k_{31} (bh_3 \, H_3(Z_2) )+k_{32} (ch_3 \, H_3(Z_3))[$]
after doing these regressions to possibly seventh order, and finding [$]k_{n1}[$] and [$]k_{n2}[$] through least squares algorithm, we could write the non-linear regression relationship solution for W as
[$]W=k_0+\displaystyle\sum_{n=1}^{7} \, k_{n1} \, bh_n \, H_n(Z_2) +\displaystyle\sum_{n=1}^{7} \, k_{n2} \, ch_n \, H_n(Z_3) [$]
Friends, I am sure that this new method and its variants would be able to find all sorts of non-linear relationships that we could not find with linear regressions. And it would also make new AI a lot more efficient.
I urge friends to play with this method. I will also include it in the applications of Z-series in my working research paper.
Suppose we want to do a nonlinear regression on W from two explanatory variables X and Y
I will give brief details to outline the new method in a few points.
First we would have to find Z-series representation of all the variables involved including regressand and regressors. Usually this would be done by calculating moments of all the data and fitting a Z-series density on each of the regressand and regressors.
For example in our case we would have.
First we find Z-Series of regressand
[$]W=a_0 + a_1 \, Z_1 + a_2 \, Z_1^2 + a_3 \, Z_1^3 + a_4 \, Z_1^4 + ... +a_7 \, Z_1^7[$]
Then we find Z-series of regressors
[$]X=b_0 + b_1 \, Z_2 + b_2 \, Z_2^2 + b_3 \, Z_2^3 + b_4 \, Z_2^4 + ... +b_7 \, Z_2^7[$]
[$]Y=c_0 + c_1 \, Z=3 + c_2 \, Z_3^2 + c_3 \, Z_3^3 + c_4 \, Z_3^4 + ... +c_7 \, Z_3^7[$]
Then we convert from Z-Series representation to hermite representation as I have done many times on this forum and also mentioned in last incomplete version of my paper above as
[$]W=ah_0 + ah_1 \, H_1(Z_1) + ah_2 \, H_2(Z_1) + ah_3 \, H_3(Z_1) + ah_4 \, H_4(Z_1) + ... +ah_7 \, H_7(Z_1)[$]
and
[$]X=bh_0 + bh_1 \, H_1(Z_2) + bh_2 \, H_2(Z_2) + bh_3 \, H_3(Z_2) + bh_4 \, H_4(Z_2) + ... +bh_7 \, H_7(Z_2)[$]
[$]Y=ch_0 + ch_1 \, H_1(Z_3) + ch_2 \, H_2(Z_3) + ch_3 \, H_3(Z_3) + ch_4 \, H_4(Z_3) + ... +ch_7 \, H_7(Z_3)[$]
The purpose of converting to hermites is that we can convert our problem of regression into several sub-problems since hermite polynomials are orthogonal to each other.
We want to regress each hermite polynomial term of the regressand on same order hermite polynomial terms of the regressors
Now we regress 1st hermite term of regressand on first hermite terms of regressors as
Here [$]k_{11}[$] and [$]k_{12}[$] are regression coefficients that we will find through ordinary least squares algorithm.
[$]ah_1 \, H_1(Z_1) =k_{11} (bh_1 \, H_1(Z_2) )+k_{12} (ch_1 \, H_1(Z_3))[$]
and for second hermite polynomial regress as
Here [$]k_{21}[$] and [$]k_{22}[$] are regression coefficients that we will find through ordinary least squares algorithm.
[$]ah_2 \, H_2(Z_1) =k_{21} (bh_2 \, H_2(Z_2) )+k_{22} (ch_2 \, H_2(Z_3))[$]
and for third polynomial ,we will regress as
[$]ah_3 \, H_3(Z_1) =k_{31} (bh_3 \, H_3(Z_2) )+k_{32} (ch_3 \, H_3(Z_3))[$]
after doing these regressions to possibly seventh order, and finding [$]k_{n1}[$] and [$]k_{n2}[$] through least squares algorithm, we could write the non-linear regression relationship solution for W as
[$]W=k_0+\displaystyle\sum_{n=1}^{7} \, k_{n1} \, bh_n \, H_n(Z_2) +\displaystyle\sum_{n=1}^{7} \, k_{n2} \, ch_n \, H_n(Z_3) [$]
Friends, I am sure that this new method and its variants would be able to find all sorts of non-linear relationships that we could not find with linear regressions. And it would also make new AI a lot more efficient.
I urge friends to play with this method. I will also include it in the applications of Z-series in my working research paper.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal