Serving the Quantitative Finance Community

Amin
Topic Author
Posts: 1880
Joined: July 14th, 2002, 3:00 am

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

Friends, first the good news: I was able to modify my earlier program so that it can handle calculation of Z-Series with proper calibration of first eight cumulants. In my preliminary checks, the program is stable and works well.
The bad news is that the supposed new method I was thinking of did not work. One of the matrices I was working with turned out to not have full rank and could not be inverted. But I am still happy that I was able to modify the earlier program and the modified program can easily tackle eight cumulants. Another good thing is that  I calibrate for first eight cumulants in one go and without initially calibrating to smaller number of moments as I did a few days ago with seventh cumulant. I will be posting the new matlab program in a day.
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

Amin
Topic Author
Posts: 1880
Joined: July 14th, 2002, 3:00 am

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

Friends, I wrote the program for Z-series calibration to first eight cumulants. The program worked at a lot of places but was still unstable at some other places. But I have added a bisection calibration method that iteratively improves the initial seed. This bisection method is very good and sometimes even finds the result by itself without Newton method. The new program is far more stable but I am still checking it to see if I need to improve anything further. Earlier without bisection method, the program would work somewhere but would sometimes stop working when input cumulant parameters would be slightly perturbed. But this new version with analytic bisection guess to Newton is far more stable. Though if cumulants input to calibration are totally out of whack, the bisection method also does not converge. I want to check the new program more carefully and then will post the matlab program here possibly later today or 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

Amin
Topic Author
Posts: 1880
Joined: July 14th, 2002, 3:00 am

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

Here are graphs for five cases from the program that I am going to upload in next post. I have given cumulant values for five cases in my program and used the same program to run those cases one after the other.
This is a far better program than any of the earlier programs but it still requires some improvement and I will continue to work on it.
You have to be careful in specifying cumulants. When this new program does not work well, some cumulant is off, problematic and not in sync with other specified cumulants. It took me quite some time to become familiar with higher cumulants. I would earlier increase fifth cumulant too fast and the whole program would go off. When you play around do not increase medium cumulants (5th-7th cumulants) too fast.
Also please look at output of Newton guess bisection program's output and compare its output moments with specified moments and you would notice that they are quite close to getting matched. You can also play around with the objective function. I am also not an expert at Newton's method and therefore please feel free to make modifications to it.
.
.

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

Amin
Topic Author
Posts: 1880
Joined: July 14th, 2002, 3:00 am

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

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

.
.

Here are the matlab program files.
.
.
function [b0,b] = FindZSeriesGivenCumulantsStandard08A(k1,k2,k3,k4,k5,k6,k7,k8)

mu=0;
sigma=1.0;

%Case1
%  k1=0.0;
%  k2=(1.0);
%  k3=-.45;%1.0;
%  k4=1.750;%-.50;
%  k5=-3;%-.5;%15;
%  k6=8;%10;%100.0;
%  k7=-10;
%  k8=200;

%Case 2
%  k1=0.0;
%  k2=(1.0);
%  k3=-.40;%1.0;
%  k4=1.5;%-.50;
%  k5=-3.80;%-.5;%15;
%  k6=11.0;%10;%100.0;
%  k7=-13.0;
%  k8=1500;

%Case 3
%  k1=0.0;
%  k2=(1.0);
%  k3=-.5;%1.0;
%  k4=1.150;%-.50;
%  k5=-2;%-.5;%15;
%  k6=12;%10;%100.0;
%  k7=-11;
%  k8=1000;

%Case 4
%  k1=0.0;
%  k2=(1.0);
%  k3=-.9;%1.0;
%  k4=1.950;%-.50;
%  k5=-5;%-.5;%15;
%  k6=19.5;%10;%100.0;
%  k7=-23;
%  k8=3850;

%Case 5
%   k1=0.0;
%  k2=(1.0);
%  k3=-.45;%1.0;
%  k4=.850;%-.50;
%  k5=-2.15;%-.5;%15;
%  k6=5.5;%10;%100.0;
%  k7=-10;
%  k8=1500;

%Case 6
k1=0.0;
k2=(1.0);
k3=-1.3;%1.0;
k4=2.50;%-.50;
k5=-3.0;%-.5;%15;
k6=8.0;%10;%100.0;
k7=-40;
k8=2500;

[mu,sigma,a0,a1,a2,a3,a4,a5,mu1,rmu1,rmu2,rmu3,rmu4,rmu5,rmu6] = GenerateDensityFromCumulantsX06Simple(mu,sigma,k1,k2,k3,k4,k5,k6)
rmu7=k1*rmu6+6*k2*rmu5+15*k3*rmu4+20*k4*rmu3+15*k5*rmu2+6*k6*rmu1+k7;
rmu8=k1*rmu7+7*k2*rmu6+21*k3*rmu5+35*k4*rmu4+35*k5*rmu3+21*k6*rmu2+7*k7*rmu1+k8;

rmu(1)=rmu1;
rmu(2)=rmu2;
rmu(3)=rmu3;
rmu(4)=rmu4;
rmu(5)=rmu5;
rmu(6)=rmu6;
rmu(7)=rmu7;
rmu(8)=rmu8;

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

Xm=median;

%str=input('Look at medians');
%Below calculate the density and its derivatives at median.

px=exp(-0.5* ((Xm-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((Xm-mu)/sigma)+  a2* (((Xm-mu)/sigma).^2-1) + ...
a3* ( ((Xm-mu)/sigma).^3 - 3*((Xm-mu)/sigma) )+ ...
a4* (((Xm-mu)/sigma).^4 - 6* ((Xm-mu)/sigma).^2 + 3)+ ...
a5* (((Xm-mu)/sigma).^5 - 10* ((Xm-mu)/sigma).^3 + 15*((Xm-mu)/sigma)));

dpxdX=exp(-0.5* ((Xm-mu)/sigma).^2)/sigma^2/sqrt(2*pi).* -(a0* ((Xm-mu)/sigma)+  a1* (((Xm-mu)/sigma).^2-1) + ...
a2* ( ((Xm-mu)/sigma).^3 - 3*((Xm-mu)/sigma) )+ ...
a3* (((Xm-mu)/sigma).^4 - 6* ((Xm-mu)/sigma).^2 + 3)+ ...
a4* (((Xm-mu)/sigma).^5 - 10* ((Xm-mu)/sigma).^3 + 15*((Xm-mu)/sigma))+ ...
a5* (((Xm-mu)/sigma).^6 - 15* ((Xm-mu)/sigma).^4 + 45*((Xm-mu)/sigma).^2-15));

d2pxdX2=exp(-0.5* ((Xm-mu)/sigma).^2)/sigma^3/sqrt(2*pi).* (a0* (((Xm-mu)/sigma).^2-1) + ...
a1* ( ((Xm-mu)/sigma).^3 - 3*((Xm-mu)/sigma) )+ ...
a2* (((Xm-mu)/sigma).^4 - 6* ((Xm-mu)/sigma).^2 + 3)+ ...
a3* (((Xm-mu)/sigma).^5 - 10* ((Xm-mu)/sigma).^3 + 15*((Xm-mu)/sigma))+ ...
a4* (((Xm-mu)/sigma).^6 - 15* ((Xm-mu)/sigma).^4 + 45*((Xm-mu)/sigma).^2-15) + ...
a5* (((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma)));

d3pxdX3=exp(-0.5* ((Xm-mu)/sigma).^2)/sigma^4/sqrt(2*pi).* -(a0* (((Xm-mu)/sigma).^3 - 3*((Xm-mu)/sigma) )+ ...
a1* (((Xm-mu)/sigma).^4 - 6* ((Xm-mu)/sigma).^2 + 3)+ ...
a2* (((Xm-mu)/sigma).^5 - 10* ((Xm-mu)/sigma).^3 + 15*((Xm-mu)/sigma))+ ...
a3* (((Xm-mu)/sigma).^6 - 15* ((Xm-mu)/sigma).^4 + 45*((Xm-mu)/sigma).^2-15) + ...
a4* (((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma))+ ...
a5* (((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105));

d4pxdX4=exp(-0.5* ((Xm-mu)/sigma).^2)/sigma^5/sqrt(2*pi).* ( a0* (((Xm-mu)/sigma).^4 - 6* ((Xm-mu)/sigma).^2 + 3)+ ...
a1* (((Xm-mu)/sigma).^5 - 10* ((Xm-mu)/sigma).^3 + 15*((Xm-mu)/sigma))+ ...
a2* (((Xm-mu)/sigma).^6 - 15* ((Xm-mu)/sigma).^4 + 45*((Xm-mu)/sigma).^2-15) + ...
a3* (((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma))+ ...
a4* (((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105)+ ...
a5* (((Xm-mu)/sigma).^9 - 36* ((Xm-mu)/sigma).^7 + 378*((Xm-mu)/sigma).^5-1260*((Xm-mu)/sigma).^3+945*((Xm-mu)/sigma)));

d5pxdX5=exp(-0.5* ((Xm-mu)/sigma).^2)/sigma^6/sqrt(2*pi).* -( a0* (((Xm-mu)/sigma).^5 - 10* ((Xm-mu)/sigma).^3 + 15*((Xm-mu)/sigma))+ ...
a1* (((Xm-mu)/sigma).^6 - 15* ((Xm-mu)/sigma).^4 + 45*((Xm-mu)/sigma).^2-15) + ...
a2* (((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma))+ ...
a3* (((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105)+ ...
a4* (((Xm-mu)/sigma).^9 - 36* ((Xm-mu)/sigma).^7 + 378*((Xm-mu)/sigma).^5-1260*((Xm-mu)/sigma).^3+945*((Xm-mu)/sigma))+ ...
a5* (((Xm-mu)/sigma).^10 - 45* ((Xm-mu)/sigma).^8 + 630*((Xm-mu)/sigma).^6-3150*((Xm-mu)/sigma).^4+4725*((Xm-mu)/sigma).^2- 945 ));

d6pxdX6=exp(-0.5* ((Xm-mu)/sigma).^2)/sigma^7/sqrt(2*pi).* (a0* (((Xm-mu)/sigma).^6 - 15* ((Xm-mu)/sigma).^4 + 45*((Xm-mu)/sigma).^2-15) + ...
a1* (((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma))+ ...
a2* (((Xm-mu)/sigma).^8 - 28* ((Xm-mu)/sigma).^6 + 210*((Xm-mu)/sigma).^4-420*((Xm-mu)/sigma).^2+105)+ ...
a3* (((Xm-mu)/sigma).^9 - 36* ((Xm-mu)/sigma).^7 + 378*((Xm-mu)/sigma).^5-1260*((Xm-mu)/sigma).^3+945*((Xm-mu)/sigma))+ ...
a4* (((Xm-mu)/sigma).^10 - 45* ((Xm-mu)/sigma).^8 + 630*((Xm-mu)/sigma).^6-3150*((Xm-mu)/sigma).^4+4725*((Xm-mu)/sigma).^2- 945 ) + ...
a5* (((Xm-mu)/sigma).^11 - 55* ((Xm-mu)/sigma).^9 + 990*((Xm-mu)/sigma).^7-6930*((Xm-mu)/sigma).^5+17325*((Xm-mu)/sigma).^3-10395*((Xm-mu)/sigma) ));

%Below calculate standard normal density and its derivatives at median.
Zm=0;
pz=exp(-0.5* (Zm.^2))/sqrt(2*pi);
dpzdZ=-Zm*exp(-0.5* (Zm.^2))/sqrt(2*pi);
d2pzdZ2=(Zm^2-1)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
d3pzdZ3=-(Zm^3-3*Zm)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
d4pzdZ4=(Zm^4-6*Zm^2+3)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
d5pzdZ5=-(Zm^5-10*Zm^3+15*Zm)*exp(-0.5* (Zm.^2))/sqrt(2*pi);
d6pzdZ6=(Zm^6-15*Zm^4+45*Zm.^2-15)*exp(-0.5* (Zm.^2))/sqrt(2*pi);

dXdZ=pz/px;

d2XdZ2=(dpzdZ-(dpxdX *dXdZ^2))/px

d3XdZ3=(d2pzdZ2- (2* dpxdX* dXdZ *d2XdZ2+dXdZ *(dXdZ^2* d2pxdX2+dpxdX *d2XdZ2)))/px

d4XdZ4=(d3pzdZ3-(3* d2XdZ2* (dXdZ^2* d2pxdX2+dpxdX *d2XdZ2)+3* dpxdX *dXdZ *d3XdZ3+dXdZ* (3 *dXdZ* d2pxdX2* d2XdZ2+dXdZ^3* d3pxdX3+dpxdX* d3XdZ3)))/px

d5XdZ5=(d4pzdZ4- (6 *(dXdZ^2 *d2pxdX2+dpxdX *d2XdZ2) *d3XdZ3+4 *d2XdZ2* (3* dXdZ *d2pxdX2 *d2XdZ2+dXdZ^3* d3pxdX3+dpxdX* d3XdZ3)+4* dpxdX* dXdZ* d4XdZ4+dXdZ *(3* d2pxdX2* d2XdZ2^2+6* dXdZ^2* d2XdZ2* d3pxdX3+4 *dXdZ* d2pxdX2* d3XdZ3+dXdZ^4* d4pxdX4+dpxdX* d4XdZ4)))/px

d6XdZ6=(d5pzdZ5- (10* d3XdZ3 *(3 *dXdZ* d2pxdX2* d2XdZ2+dXdZ^3* d3pxdX3+dpxdX *d3XdZ3)+10* (dXdZ^2* d2pxdX2+dpxdX *d2XdZ2)* d4XdZ4+5* d2XdZ2* (3* d2pxdX2* d2XdZ2^2+6* dXdZ^2* d2XdZ2* d3pxdX3+4 *dXdZ* d2pxdX2* d3XdZ3+dXdZ^4* d4pxdX4+dpxdX* d4XdZ4)+5 *dpxdX* dXdZ* d5XdZ5+dXdZ *(15* dXdZ* d2XdZ2^2* d3pxdX3+10* d2pxdX2* d2XdZ2* d3XdZ3+10* dXdZ^2 *d3pxdX3* d3XdZ3+10* dXdZ^3* d2XdZ2* d4pxdX4+5* dXdZ* d2pxdX2* d4XdZ4+dXdZ^5* d5pxdX5+dpxdX* d5XdZ5)))/px

d7XdZ7=(d6pzdZ6- (20 *(3* dXdZ* d2pxdX2* d2XdZ2+dXdZ^3 *d3pxdX3+dpxdX *d3XdZ3)* d4XdZ4+15* d3XdZ3* (3 *d2pxdX2 *d2XdZ2^2+6* dXdZ^2* d2XdZ2* d3pxdX3+4 *dXdZ* d2pxdX2* d3XdZ3+dXdZ^4* d4pxdX4+dpxdX* d4XdZ4)+15* (dXdZ^2* d2pxdX2+dpxdX* d2XdZ2)* d5XdZ5+6* d2XdZ2* (15* dXdZ* d2XdZ2^2* d3pxdX3+10* d2pxdX2* d2XdZ2 *d3XdZ3+10 *dXdZ^2 *d3pxdX3 *d3XdZ3+10 *dXdZ^3 *d2XdZ2* d4pxdX4+5 *dXdZ *d2pxdX2 *d4XdZ4+dXdZ^5 *d5pxdX5+dpxdX* d5XdZ5)+6* dpxdX* dXdZ *d6XdZ6+dXdZ* (15 *d2XdZ2^3 *d3pxdX3+60 *dXdZ *d2XdZ2 *d3pxdX3* d3XdZ3+10 *d2pxdX2 *d3XdZ3^2+45 *dXdZ^2* d2XdZ2^2 *d4pxdX4+20* dXdZ^3 *d3XdZ3* d4pxdX4+15 *d2pxdX2* d2XdZ2 *d4XdZ4+15 *dXdZ^2* d3pxdX3 *d4XdZ4+15 *dXdZ^4* d2XdZ2* d5pxdX5+6 *dXdZ *d2pxdX2* d5XdZ5+dXdZ^6* d6pxdX6+dpxdX *d6XdZ6)))/px

%Below calculate iniotial guess coefficients of Z-series to be input to
%Newton method.

c0=median;
c1=dXdZ;
c2=1/2*1*d2XdZ2;
c3=1/6*1*d3XdZ3;
c4=1/24*1*d4XdZ4;
c5=1/120*1*d5XdZ5;
c6=1/720*1*d6XdZ6;
c7=1/720/7*d7XdZ7;

%earlier program worked only because odd coefficients that are related to
%even moments were reasonablyclose to true values. Even coefficients
%responsible for odd moments were usually grossly off. Below I try to make
%sure that sign of even coefficients matches the sign of odd moments. This
%is not completely definitive but usually right.
cc0=c0;
cc(1)=c1;
if(sign(k3)~=sign(c2))
cc(2)=-1*c2;
else
cc(2)=c2;
end
cc(3)=c3;
if(sign(k5)~=sign(c4))
cc(4)=-1*c4;
else
cc(4)=c4;
end
cc(5)=c5;
if(sign(k7)~=sign(c6))
cc(6)=-1*c6;
else
cc(6)=c6;
end

cc(7)=c7;

cc0=-(cc(2)+cc(4)*3+cc(6)*15);

SeriesOrder=7;
NMoments=8;

[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
rmu

str=input('Look at moments--before')

%Below is the function that finds a guess for Newton method using heuristic
%bisection.

[cc0,cc] = ImproveNewtonGuess02(cc0,cc,rmu,25);

%str=input('Look at initial match of moments for the Newton guess');
%Below specify target cumulants.
C(1:NMoments+1)=0;
C(1)=k1;
C(2)=k2;
C(3)=k3;
C(4)=k4;
C(5)=k5;
C(6)=k6;
C(7)=k7;
C(8)=k8;

%Below calculate F, the defect vector from target cumulants. And calculate
%dF which is the matrix consisting of derivatives of defect vector from
%target cumulants.

[F,dF] = CalculateCumulantsAndDerivativesFromMoments_0(C,cc0,cc,SeriesOrder,NMoments,8)

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

NoOfCumulants=NMoments
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments);
ObjBest=0;

%ObjBest=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
%   abs((rmu(5)-Moments(5))^(1/2.5))+abs((rmu(6)-Moments(6))^(1/3.0))+abs((rmu(7)-Moments(7))^(1/3.5))+ ...
%   abs((rmu(8)-Moments(8))^(1/4.0));

ObjBest=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
abs((rmu(5)-Moments(5))^(1/2.0))+abs((rmu(6)-Moments(6))^(1/2.0))+abs((rmu(7)-Moments(7))^(1/2.0))+ ...
abs((rmu(8)-Moments(8))^(1/2.0));

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

nn=0;
while((nn<100)&&((abs(F(1,1))>.000000000001) || (abs(F(2,1))>.000000000001) || (abs(F(3,1))>.000000000001) || (abs(F(4,1))>.00000000001)|| (abs(F(5,1))>.00000000001)|| (abs(F(6,1))>.00000000001)  ))

nn=nn+1;
%Below Newton matrix equation to improve the Z-series coefficients guess at previous step.

da=da-dF\F;

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

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

[IsValidFlag] = CheckIsValidDensity(b0,b)

[Moments] = CalculateMomentsOfZSeries(b0,b,SeriesOrder,NMoments);
ObjNew=0;
%ObjNew=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
%   abs((rmu(5)-Moments(5))^(1/2.5))+abs((rmu(6)-Moments(6))^(1/3.0))+abs((rmu(7)-Moments(7))^(1/3.5))+ ...
%   abs((rmu(8)-Moments(8))^(1/4.0));

ObjNew=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
abs((rmu(5)-Moments(5))^(1/2.0))+abs((rmu(6)-Moments(6))^(1/2.0))+abs((rmu(7)-Moments(7))^(1/2.0))+ ...
abs((rmu(8)-Moments(8))^(1/2.0));

if((ObjBest>ObjNew) &&( IsValidFlag))

ObjBest=ObjNew;
b0Best=b0;
bBest(1:SeriesOrder)=b(1:SeriesOrder);
end

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

end
b0=b0Best;%Best;
b(1:SeriesOrder)=bBest;%Best(1:SeriesOrder);

%Above are the best chosen Z-series coefficients based on the objective function.
%But my objective function can be improved.

[Moments] = CalculateMomentsOfZSeries(b0,b,SeriesOrder,NMoments)
%Above are raw moments calculated from the resulting Z-series from Newton
%Method

%Below are raw moments that were input to calibration.
rmu
str=input('Look at comparison of moments calibration');

%Below we do calculations to construct the density from our newly derived
%Z-series.

dNn=.2/2;   % Normal density subdivisions width. would change with number of subdivisions
Nn=90;  % No of normal density subdivisions

NnMid=4.0;

Z(1:Nn)=(((1:Nn)-6.5/1)*dNn-NnMid);
%Z
%str=input('Look at Z');
%ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
%ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);
%ZProb(2:Nn-1)=normcdf(.5*Z(2:Nn-1)+.5*Z(3:Nn),0,1)-normcdf(.5*Z(2:Nn-1)+.5*Z(1:Nn-2),0,1);

w(1:Nn)=b0;
for nn=1:SeriesOrder
w(1:Nn)=w(1:Nn)+b(nn)*Z(1:Nn).^nn;
end
wnStart=1;
Dfw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1

Dfw(nn) = (w(nn + 1) - w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
end

pw(1:Nn)=0;
for nn = wnStart+1:Nn-1

pw(nn) = (normpdf(Z(nn),0, 1))/abs(Dfw(nn));

end

%
% %Below is the construction of density of our originial analytic density for
% %plotting purposes
%
dMm=.02/10;
Mm=1250*5;
X(1:Mm)=-5+dMm+dMm*(1:Mm);

%
% %Below is the newly constructed density so that its output central moments
% %match the central moments input by the user.
% pdf=exp(-0.5* ((X-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X-mu)/sigma)+  a2* (((X-mu)/sigma).^2-1) + ...
%     a3* ( ((X-mu)/sigma).^3 - 3*((X-mu)/sigma) )+ ...
%     a4* (((X-mu)/sigma).^4 - 6* ((X-mu)/sigma).^2 + 3)+ ...
%     a5* (((X-mu)/sigma).^5 - 10* ((X-mu)/sigma).^3 + 15*((X-mu)/sigma))+ ...
%     a6* (((X-mu)/sigma).^6 - 15* ((X-mu)/sigma).^4 + 45*((X-mu)/sigma).^2-15))+ ...
%     a7*(((Xm-mu)/sigma).^7 - 21* ((Xm-mu)/sigma).^5 + 105*((Xm-mu)/sigma).^3-105*((Xm-mu)/sigma));
%
% %Below is standard Gaussian PDF.
pdfG=exp(-0.5* ((X-mu)/sigma).^2)/sigma/sqrt(2*pi);
%
plot(w(wnStart+1:Nn-1),pw(wnStart+1:Nn-1),'r',X(1:Mm),pdfG(1:Mm),'b');

title(sprintf('Fm1 = %.4f,m1=%.4f,Fm2=%.2f,m2=%.2f,Fm3=%.2f,m3 =%.2f,Fm4=%.2f,m4=%.2f,Fm5=%.2f,m5=%.2f,Fm6=%.1f,m6=%.1f,Fm7=%.1f,m7=%.1f,Fm8=%.1f,m8=%.1f',Moments(1),rmu1,Moments(2),rmu2,Moments(3),rmu3,Moments(4),rmu4,Moments(5),rmu5,Moments(6),rmu6,Moments(7),rmu7,Moments(8),rmu8));
%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

legend({'Newton Density','Analytical Density','Standard Gaussian Density'},'Location','northeast')

cc0
b0
cc
b

end


.
.
.
.
function [cc0Best,ccBest] = ImproveNewtonGuess02(cc0,cc,rmu,Iterations)
%This function finds a good guess for Newton method by heuristic bisection.

SeriesOrder=7;
NMoments=8;
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments);

%ObjectiveBest=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
%    abs((rmu(5)-Moments(5))^(1/2.5))+abs((rmu(6)-Moments(6))^(1/3.0))+abs((rmu(7)-Moments(7))^(1/3.5))+ ...
%    abs((rmu(8)-Moments(8))^(1/4.0));

ObjectiveBest=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
abs((rmu(5)-Moments(5))^(1/2.0))+abs((rmu(6)-Moments(6))^(1/2.0))+abs((rmu(7)-Moments(7))^(1/2.0))+ ...
abs((rmu(8)-Moments(8))^(1/2.0));

ccBest=cc;
cc0Best=cc0;
r1=.5;%bisection weight one
r2=.5;%bisection weight two

for nn=1:Iterations

if(nn>1)
%r1=r1+.025;
%r2=r2-.025;
end

[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
Ratio2=(r1*rmu(2)+r2*Moments(2))/Moments(2);
cc(1)=cc(1)*sqrt(Ratio2);

[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)

Ratio3=(r1*rmu(3)+r2*Moments(3))/Moments(3);
cc(2)=cc(2)*(Ratio3);

cc0=-(cc(2)+3*cc(4)+15*cc(6));

[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)

Ratio4=(r1*rmu(4)+r2*Moments(4))/Moments(4);
cc(3)=cc(3)*(Ratio4);

[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)

Ratio5=(r1*rmu(5)+r2*Moments(5))/Moments(5);
cc(4)=cc(4)*(Ratio5);

cc0=-(cc(2)+3*cc(4)+15*cc(6));
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
Ratio6=(r1*rmu(6)+r2*Moments(6))/Moments(6);
cc(5)=cc(5)*(Ratio6);

[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
Ratio7=(r1*rmu(7)+r2*Moments(7))/Moments(7);
cc(6)=cc(6)*(Ratio7);

[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
Ratio8=(r1*rmu(8)+r2*Moments(8))/Moments(8);
cc(7)=cc(7)*(Ratio8);

cc0=-(cc(2)+3*cc(4)+15*cc(6));
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)

[IsValidFlag] = CheckIsValidDensity(cc0,cc)

%ObjectiveNew=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
%    abs((rmu(5)-Moments(5))^(1/2.5))+abs((rmu(6)-Moments(6))^(1/3.0))+abs((rmu(7)-Moments(7))^(1/3.5))+ ...
%    abs((rmu(8)-Moments(8))^(1/4.0));

ObjectiveNew=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
abs((rmu(5)-Moments(5))^(1/2.0))+abs((rmu(6)-Moments(6))^(1/2.0))+abs((rmu(7)-Moments(7))^(1/2.0))+ ...
abs((rmu(8)-Moments(8))^(1/2.0));

if((ObjectiveNew<ObjectiveBest) &&( IsValidFlag))
ObjectiveBest=ObjectiveNew;

ccBest=cc;
cc0Best=cc0;
end

end

[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments)
rmu
str=input('Look at moments--In function---After');

end


.
.
.
.
function [IsValidFlag] = CheckIsValidDensity(cc0,cc)
%Roughly checks if the density is valid. It is not a definitive check but
%just a rough one but mostly works.

Z=2.0;
Xp1=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;

Z=3.0;
Xp2=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;

Z=4.65;
Xp3=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;

Z=-2.0;
Xn1=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;

Z=-3.0;
Xn2=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;

Z=-4.65;
Xn3=cc0 + cc(1) * Z+ cc(2) * Z^2+ cc(3) * Z^3+ cc(4) * Z^4+ cc(5) * Z^5+ cc(6) * Z^6+ cc(7) * Z^7;

IsValidFlag=0;
if( ((Xp2>Xp1) && (Xp3>Xp2)) &&(Xp1>Xn1) && ((Xn2<Xn1) && (Xn3<Xn2)))
IsValidFlag=1;
end

end


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

%In this program, you specify the first six cumulants through input parametersand the program creates a density based on these
%cumulantsents.

%The algorithm takes a base gaussian density that has it mean mu and
%sq-root of variance given by sigma.
%This mean and sigma can be different from first and sq-root of second
%central moment(variance) of the target density that we specified as mu1 and mu2 earlier.

%Below is the construction of grid to calculate pdf, numerically integrate the pdf
%to find cdf. And to numerically calculate the six moments on the constructed density
%so they can be matched with input moments.

dNn=.02/10;
Nn=1250*10;
X(1:Nn)=-10+dNn+dNn*(1:Nn);

%Below, we calculate first six raw-moments from the input cumulants.
%These raw moments would be used to calculate the coefficients of hermite
%polynomials so that output central moments calculated on new density match
%the input central moments.

mu1=k1;
rmu1=k1;
rmu2=k1*rmu1+k2;
rmu3=k1*rmu2+2*k2*rmu1+k3;
rmu4=k1*rmu3+3*k2*rmu2+3*k3*rmu1+k4;
rmu5=k1*rmu4+4*k2*rmu3+6*k3*rmu2+4*k4*rmu1+k5;
rmu6=k1*rmu5+5*k2*rmu4+10*k3*rmu3+10*k4*rmu2+5*k5*rmu1+k6;
%rmu7=k1*rmu6+6*k2*rmu5+15*k3*rmu4+20*k4*rmu3+15*k5*rmu2+6*k6*rmu1+k7;
%rmu8=k1*rmu7+7*k2*rmu6+21*k3*rmu5+35*k4*rmu4+35*k5*rmu3+21*k6*rmu2+7*k7*rmu1+k8;

%Below is base gaussian pdf.
pdfG=exp(-0.5* ((X-mu)/sigma).^2)/sigma/sqrt(2*pi);

%Below I calculate the coefficients of polynomials making the density so
%that the resulting density has the same as specified raw moments.
%This follows my description of method 4 in the post 1491 given at
%https://forum.wilmott.com/viewtopic.php?f=4&t=99702&p=871093#p871047

MOrder=6+6
SeriesOrder=6;
SeriesOrder=6;
[EZ] = CalculateGaussianMoments(MOrder)

%str=input('Look at Z  Expectations')
CMatrix(1,1:SeriesOrder)=EZ(1:SeriesOrder);
CMatrix(2,1:SeriesOrder)=EZ(2:SeriesOrder+1);
CMatrix(3,1:SeriesOrder)=EZ(3:SeriesOrder+2);
CMatrix(4,1:SeriesOrder)=EZ(4:SeriesOrder+3);
CMatrix(5,1:SeriesOrder)=EZ(5:SeriesOrder+4);
CMatrix(6,1:SeriesOrder)=EZ(6:SeriesOrder+5);
%CMatrix(7,1:SeriesOrder)=EZ(7:SeriesOrder+6);
%CMatrix(8,1:SeriesOrder)=EZ(8:SeriesOrder+7);

Mu(1,1)=1.0;
Mu(2,1)=rmu1;
Mu(3,1)=rmu2;
Mu(4,1)=rmu3;
Mu(5,1)=rmu4;
Mu(6,1)=rmu5;
%Mu(7,1)=rmu6;
%Mu(8,1)=rmu7;

Coeff=CMatrix\Mu;

Coeff

%str=input('Look at Mu results');
%Below I convert from simple polynomials to hermite polynomials so that
%resulting density remains identical with the initial density.

[CoeffaH] = ConvertZCoeffsToHCoeffsAB(Coeff,SeriesOrder)

a0=CoeffaH(1);
a1=CoeffaH(2);
a2=CoeffaH(3);
a3=CoeffaH(4);
a4=CoeffaH(5);
a5=CoeffaH(6);
a6=0;%CoeffaH(7);
a7=0;%CoeffaH(8);
a8=0;

CoeffaH(7)=0;
CoeffaH(8)=0;

%str=input('Look at Mu results---2');

%Below is the newly constructed density so that its output raw moments
%match the raw moments input by the user.
%Sometimes this density is negative at places but we still use it to derive
%an initial guess for construction of Z-series Newton density.

pdf=exp(-0.5* ((X-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X-mu)/sigma)+  a2* (((X-mu)/sigma).^2-1) + ...
a3* ( ((X-mu)/sigma).^3 - 3*((X-mu)/sigma) )+ ...
a4* (((X-mu)/sigma).^4 - 6* ((X-mu)/sigma).^2 + 3)+ ...
a5* (((X-mu)/sigma).^5 - 10* ((X-mu)/sigma).^3 + 15*((X-mu)/sigma))+ ...
a6* (((X-mu)/sigma).^6 - 15* ((X-mu)/sigma).^4 + 45*((X-mu)/sigma).^2-15)+ ...
a7* (((X-mu)/sigma).^7 - 21* ((X-mu)/sigma).^5 + 105*((X-mu)/sigma).^3-105*((X-mu)/sigma))+ ...
a8* (((X-mu)/sigma).^8 - 28* ((X-mu)/sigma).^6 + 210*((X-mu)/sigma).^4-420*((X-mu)/sigma).^2+105));

%Below cdf is calculated by doing numerical integration of the density. I
%have replaced integration with summations.
cdf0(1:Nn)=0;
for nn=2:Nn
cdf0(nn)=cdf0(nn-1)+exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2/2* (((X(nn)-mu)/sigma).^2-1) + ...
a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15) + ...
a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
end

%Below is numerically calculated value of first raw moment. N in mu1N
%stands for numerical.

rmu1N=0;
for nn=2:Nn
rmu1N=rmu1N+X(nn) .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2/2* (((X(nn)-mu)/sigma).^2-1) + ...
a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15)+ ...
a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
end

%Below is numerically calculated value of second raw moment.

rmu2N=0;
for nn=2:Nn
rmu2N=rmu2N+(X(nn)-0*rmu1N).^2 .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2* (((X(nn)-mu)/sigma).^2-1) + ...
a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15)+ ...
a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
end

%Below is numerically calculated value of third raw moment.
rmu3N=0;
for nn=2:Nn
rmu3N=rmu3N+(X(nn)-0*rmu1N).^3 .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2* (((X(nn)-mu)/sigma).^2-1) + ...
a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15)+ ...
a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
end

%Below is numerically calculated value of fourth raw moment.
rmu4N=0;
for nn=2:Nn
rmu4N=rmu4N+(X(nn)-0*rmu1N).^4 .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2* (((X(nn)-mu)/sigma).^2-1) + ...
a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15)+ ...
a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
end

%Below is numerically calculated value of fifth raw moment.
rmu5N=0;
for nn=2:Nn
rmu5N=rmu5N+(X(nn)-0*rmu1N).^5 .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2* (((X(nn)-mu)/sigma).^2-1) + ...
a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15) + ...
a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
end

% %Below is numerically calculated value of sixth raw moment.
% rmu6N=0;
% for nn=2:Nn
% rmu6N=rmu6N+(X(nn)-0*rmu1N).^6 .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2* (((X(nn)-mu)/sigma).^2-1) + ...
%     a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
%     a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
%     a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
%     a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15) + ...
%     a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
%     a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
% end
%
% rmu7N=0;
% for nn=2:Nn
% rmu7N=rmu7N+(X(nn)-0*rmu1N).^7 .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2* (((X(nn)-mu)/sigma).^2-1) + ...
%     a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
%     a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
%     a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
%     a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15) + ...
%     a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
%     a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
% end
%
% rmu8N=0;
% for nn=2:Nn
% rmu8N=rmu8N+(X(nn)-0*rmu1N).^8 .*exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X(nn)-mu)/sigma)+ a2* (((X(nn)-mu)/sigma).^2-1) + ...
%     a3* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
%     a4* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
%     a5* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
%     a6* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15) + ...
%     a7* (((X(nn)-mu)/sigma).^7 - 21* ((X(nn)-mu)/sigma).^5 + 105*((X(nn)-mu)/sigma).^3-105*((X(nn)-mu)/sigma))+ ...
%     a8* (((X(nn)-mu)/sigma).^8 - 28* ((X(nn)-mu)/sigma).^6 + 210*((X(nn)-mu)/sigma).^4-420*((X(nn)-mu)/sigma).^2+105))*dNn;
% end
%Below is analytically calculated cdf.

cdf(1:Nn)=normcdf(((X(1:Nn)-mu)/sigma))*a0+exp(-0.5* ((X(1:Nn)-mu)/sigma).^2)/sqrt(2*pi).* (- a1 - a2* ((X(1:Nn)-mu)/sigma) - ...
a3* ( ((X(1:Nn)-mu)/sigma).^2 - 1 )- ...
a4* (((X(1:Nn)-mu)/sigma).^3 - 3* ((X(1:Nn)-mu)/sigma)) - ...
a5* (((X(1:Nn)-mu)/sigma).^4 - 6* ((X(1:Nn)-mu)/sigma).^2 + 3) - ...
a6* (((X(1:Nn)-mu)/sigma).^5 - 10* ((X(1:Nn)-mu)/sigma).^3 + 15*((X(1:Nn)-mu)/sigma)));

pdf02=exp(-0.5* ((X-mu)/sigma).^2)/sigma/sqrt(2*pi).* (CoeffaH(1)+ CoeffaH(2)* ((X-mu)/sigma)+  CoeffaH(3)* (((X-mu)/sigma).^2-1) + ...
CoeffaH(4)* ( ((X-mu)/sigma).^3 - 3*((X-mu)/sigma) )+ ...
CoeffaH(5)* (((X-mu)/sigma).^4 - 6* ((X-mu)/sigma).^2 + 3)+ ...
CoeffaH(6)* (((X-mu)/sigma).^5 - 10* ((X-mu)/sigma).^3 + 15*((X-mu)/sigma))+ ...
CoeffaH(7)* (((X-mu)/sigma).^6 - 15* ((X-mu)/sigma).^4 + 45*((X-mu)/sigma).^2-15)+ ...
CoeffaH(8)* (((X-mu)/sigma).^7 - 21* ((X-mu)/sigma).^5 + 105*((X-mu)/sigma).^3-105*((X-mu)/sigma)));%+ ...
%a8/720/56* (((X-mu)/sigma).^8 - 28* ((X-mu)/sigma).^6 + 210*((X-mu)/sigma).^4-420*((X-mu)/sigma).^2+105));

plot(X(1:Nn),cdf0(1:Nn),'g',X(1:Nn),cdf(1:Nn),'r',X(1:Nn),pdf(1:Nn),'b',X(1:Nn),pdfG(1:Nn),'m',X(1:Nn),pdf02(1:Nn),'k');
legend({'Numerically Calculated CDF','Analytically Calculated CDF','Newly Constructed Moment Matched PDF','Base Gaussian PDF'},'Location','northeast')

%Below program checks if density becomes slightly negative.
IsNegative=0;
for nn=2:Nn
if(pdf(nn)<0)
IsNegative=1;
end
end

%cdf
%Below compare the input raw moments rmu's with moments calculated on
%newly constructed density denoted by rmu#N
rmu1N
rmu1
rmu2N
rmu2
rmu3N
rmu3
rmu4N
rmu4
rmu5N
rmu5
%rmu6N
%rmu6
%rmu7N
%rmu7
%rmu8N
%rmu8

%str=input('Look at match of moments');

%Below are coefficients on hermite polynomials found out so that output
%moments in newly constructed density match input central moment values.

a1
a2
a3
a4
a5
%a6
%a7
%a8

IsNegative

%str=input('Look if the density is negative');

end


.
.
.
function [median] = FindMedianHermiteDensity06(mu1,mu,sigma,a0,a1,a2,a3,a4,a5)
%Uses Newton-Raphson to calculate the median.

X=mu1;

cdf=normcdf(((X-mu)/sigma))*a0+exp(-0.5* ((X-mu)/sigma).^2)/sqrt(2*pi).* (- a1 - a2* ((X-mu)/sigma) - ...
a3* ( ((X-mu)/sigma).^2 - 1 )- ...
a4* (((X-mu)/sigma).^3 - 3* ((X-mu)/sigma)) - ...
a5* (((X-mu)/sigma).^4 - 6* ((X-mu)/sigma).^2 + 3));% - ...
%a6/720* (((X-mu)/sigma).^5 - 10* ((X-mu)/sigma).^3 + 15*((X-mu)/sigma)));

nn=0;

while (abs(cdf-.5)> .00000000001)
nn=nn+1;
pdf=exp(-0.5* ((X-mu)/sigma).^2)/sigma/sqrt(2*pi).* (a0+ a1* ((X-mu)/sigma)+  a2* (((X-mu)/sigma).^2-1) + ...
a3* ( ((X-mu)/sigma).^3 - 3*((X-mu)/sigma) )+ ...
a4* (((X-mu)/sigma).^4 - 6* ((X-mu)/sigma).^2 + 3)+ ...
a5* (((X-mu)/sigma).^5 - 10* ((X-mu)/sigma).^3 + 15*((X-mu)/sigma)));%+ ...
%a6/720* (((X-mu)/sigma).^6 - 15* ((X-mu)/sigma).^4 + 45*((X-mu)/sigma).^2-15));

X=X-(cdf-.5)/pdf;
cdf=normcdf(((X-mu)/sigma))*a0+exp(-0.5* ((X-mu)/sigma).^2)/sqrt(2*pi).* (- a1 - a2* ((X-mu)/sigma) - ...
a3* ( ((X-mu)/sigma).^2 - 1 )- ...
a4* (((X-mu)/sigma).^3 - 3* ((X-mu)/sigma)) - ...
a5* (((X-mu)/sigma).^4 - 6* ((X-mu)/sigma).^2 + 3));% - ...
%a6/720* (((X-mu)/sigma).^5 - 10* ((X-mu)/sigma).^3 + 15*((X-mu)/sigma)));

cdf
end

nn
cdf
mu1
median=X

% dNn=.002/100;
% Nn=1250*1000;
% X(1:Nn)=-10+dNn+dNn*(1:Nn);
%
% Mm=round((median+10-dNn)/dNn);
%
%
%
% cdf0(1:Nn)=0;
% for nn=2:Mm
% cdf0(nn)=cdf0(nn-1)+exp(-0.5* ((X(nn)-mu)/sigma).^2)/sigma/sqrt(2*pi).* (1+ a1* ((X(nn)-mu)/sigma)+ a2/2* (((X(nn)-mu)/sigma).^2-1) + ...
%     a3/6* ( ((X(nn)-mu)/sigma).^3 - 3*((X(nn)-mu)/sigma) )+ ...
%     a4/24* (((X(nn)-mu)/sigma).^4 - 6* ((X(nn)-mu)/sigma).^2 + 3) + ...
%     a5/120* (((X(nn)-mu)/sigma).^5 - 10* ((X(nn)-mu)/sigma).^3 + 15*((X(nn)-mu)/sigma))+ ...
%     a6/720* (((X(nn)-mu)/sigma).^6 - 15* ((X(nn)-mu)/sigma).^4 + 45*((X(nn)-mu)/sigma).^2-15))*dNn;
% end
%
%
% cdf0(Mm)
%
% str=input('Look at median cdf');

end


.
.
.
function [F,dF] = CalculateCumulantsAndDerivativesFromMoments_0(C,a0,a,SeriesOrder,NZterms,NMoments)
%[F,dF] = CalculateCumulantsAndDerivativesFromMoments(C,a0,a,SeriesOrder,SeriesOrder,NoOfCumulants);

%a(1)=a1;
%a(2)=a2;
%a(3)=a3;
%a(4)=a4;
%a(5)=a5;
%a(6)=a6;
%a(7)=a7;

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

EZ(1)=0;
EZ(2)=1;
for nn=3:NMoments*SeriesOrder+NZterms+2
if rem(nn,2)==1
EZ(nn)=0;
else
EZ(nn)=EZ(nn-2)*(nn-1);
EZ(nn);
end
end
EZ;

EXZ(1,1)=1;
for pp1=1:NZterms
EXZ(1,pp1+1)=EZ(pp1);
end

a(SeriesOrder+1:NMoments*SeriesOrder+1)=0;
b0=a0;
b=a;

for mm=1:NMoments
if(mm>1)
[b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*mm);
b(SeriesOrder*mm+1:NMoments*SeriesOrder+1)=0;
end
% b0
% b
%str=input('Look at numbers')
EXZ(mm+1,1)=b0;
for pp2=1:SeriesOrder*mm

EXZ(mm+1,1)=EXZ(mm+1,1)+b(pp2).*EZ(pp2);
end
for pp1=1:NZterms
EXZ(mm+1,pp1+1)=b0.*EZ(pp1);
for pp2=1:SeriesOrder*mm

EXZ(mm+1,pp1+1)=EXZ(mm+1,pp1+1)+b(pp2).*EZ(pp2+pp1);
end
end
end

u1=EXZ(2,1);
u2=EXZ(3,1);
u3=EXZ(4,1);
u4=EXZ(5,1);

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

k2=u2-u1^2;
k3=u3-3*u2*u1+2*u1^3;
k4=u4-4*u3*u1-3*u2^2+12*u2*u1^2-6*u1^4;

% du1(1)=1;%----2
% du2(1)=2*EXZ(2,1);
% du3(1)=3*EXZ(3,1);
% du4(1)=4*EXZ(4,1);

du1(1)=0;%----2
du2(1)=0;
du3(1)=0;
du4(1)=0;

for mm=2:SeriesOrder+1
du1(mm)=EXZ(1,mm);
du2(mm)=2*EXZ(2,mm);
du3(mm)=3*EXZ(3,mm);
du4(mm)=4*EXZ(4,mm);

end

if(NMoments>=5)
u5=EXZ(6,1);
k5=u5-5*u4*u1-10*u3*u2+20*u3*u1^2+30*u2^2*u1-60*u2*u1^3+24*u1^5;
du5(1)=0;

for mm=2:SeriesOrder+1
du5(mm)=5*EXZ(5,mm);
end

end

if(NMoments>=6)
u6=EXZ(7,1);
k6=u6-6*u5*u1-15*u4*u2+30*u4*u1^2-10*u3^2+120*u3*u2*u1-120*u3*u1^3+30*u2^3 - ...
270*u2^2*u1^2+360*u2*u1^4-120*u1^6;
du6(1)=0;

for mm=2:SeriesOrder+1
du6(mm)=6*EXZ(6,mm);
end

end

if(NMoments>=7)
u7=EXZ(8,1);
%k7=u7-u1*u6-6*k2*u5-15*k3*u4-20*k4*u3-15*k5*u2-6*k6*u1
k7= u7- 7 *u1* u6 + 42* u1^2* u5 - 21* u2 *u5- 210* u1^3 *u4+ ...
210* u1*u2* u4 - 35* u3 *u4 + 840 *u1^4* u3 + 210* u2^2* u3 + ...
210*12* u1^3* u2^2   - 1260* u1^2 * u2* u3  - 630 *u1* u2^3 + 140*u1* u3^2  - 2520* u1^5* u2 + 720* u1^7 ;
%str=input('Look at k7 and k71')

du7(1)=0;

for mm=2:SeriesOrder+1
du7(mm)=7*EXZ(7,mm);
end

end
if(NMoments>=8)

u8=EXZ(9,1);
k8 = u8 - u1* u7 - 7 *k2* u6 - 21* k3* u5 - 35* k4* u4 - ...
35* k5* u3 - 21* k6* u2 - 7 *k7* u1;
du8(1)=0;

for mm=2:SeriesOrder+1
du8(mm)=8*EXZ(8,mm);
end

end

% for mm=2:SeriesOrder+1
% du1(mm)=EXZ(1,mm);
% du2(mm)=2*EXZ(2,mm);
% du3(mm)=3*EXZ(3,mm);
% du4(mm)=4*EXZ(4,mm);
% du5(mm)=5*EXZ(5,mm);
% du6(mm)=6*EXZ(6,mm);
% du7(mm)=7*EXZ(7,mm);
% du8(mm)=8*EXZ(8,mm);
% end

dk1(1)=1.0;
dk1(2)=0.0;
dk1(3)=1.0;
dk1(4)=0.0;
if(SeriesOrder>=4)
dk1(5)=3.0;
end
if(SeriesOrder>=5)
dk1(6)=0.0;
end
if(SeriesOrder>=6)
dk1(7)=15.0;
end
if(SeriesOrder>=7)
dk1(8)=0.0;
end
if(SeriesOrder>=8)
dk1(9)=105.0;
end

% dk2(1)=0;
% dk3(1)=0;
% dk4(1)=0;
% dk5(1)=0;
% dk6(1)=0;
% dk7(1)=0;
% dk8(1)=0;

for mm=2:SeriesOrder+1
%dk1(mm)=du1(mm);
dk2(mm)=-2.* u1.*du1(mm)+du2(mm);
dk3(mm)=du3(mm)-3*du1(mm)*u2-3*u1*du2(mm)+2*3*u1^2*du1(mm);
dk4(mm)=du4(mm)-4*du3(mm)*u1-4*u3*du1(mm)-3*2*u2*du2(mm)+12*du2(mm)*u1^2+12*2*u2*u1*du1(mm)-6*4*u1^3*du1(mm);

if(NMoments>=5)
dk5(mm)=du5(mm)-5*du4(mm)*u1-5*u4*du1(mm)-10*du3(mm)*u2-10*u3*du2(mm)+20*du3(mm)*u1^2+20*2*u3*u1*du1(mm)+ ...
30*2*u2*du2(mm)*u1+30*u2^2*du1(mm)-60*du2(mm)*u1^3-60*3*u2*u1^2*du1(mm)+24*5*u1^4*du1(mm);
end
if(NMoments>=6)
dk6(mm)=du6(mm)-6*du5(mm)*u1-6*u5*du1(mm)-15*du4(mm)*u2-15*u4*du2(mm)+30*du4(mm)*u1^2+30*2*u4*u1*du1(mm)-10*2*u3*du3(mm)+ ...
120*du3(mm)*u2*u1-120*u3*du2(mm)*u1-120*u3*u2*du1(mm)-120*du3(mm)*u1^3-120*3*u3*u1^2*du1(mm)+30*3*u2^2*du2(mm) - ...
270*2*u2*du2(mm)*u1^2+270*2*u2^2*u1*du1(mm)+360*du2(mm)*u1^4+360*4*u2*u1^3*du1(mm)-120*6*u1^5*du1(mm);
end
if(NMoments>=7)
dk7(mm)= du7(mm)- 7 *du1(mm)* u6- 7 *u1* du6(mm) + 42*2* du1(mm)* u5+ 42* u1^2* du5(mm)  ...
- 21* du2(mm) *u5- 21* u2 *du5(mm)- 210* 3*u1^2*du1(mm) *u4- 210* u1^3 *du4(mm) + ...
210* du1(mm)*u2* u4+210* u1*du2(mm)* u4+210* u1*u2* du4(mm) - 35* du3(mm) *u4- 35* u3 *du4(mm) ...
+ 840 *4*u1^3*du1(mm)* u3 + 840 *u1^4* du3(mm)+ 210* 2*u2*du2(mm)* u3 + 210* u2^2* du3(mm)+ ...
210*12* 3*u1^2*du1(mm)* u2^2 +210*12* u1^3* 2*u2*du2(mm)   - 1260* 2*u1*du1(mm) * u2* u3- 1260* u1^2 * du2(mm)* u3 ...
- 1260* u1^2 * u2* du3(mm) - 630 *du1(mm)* u2^3- 630 *u1* 3*u2^2*du2(mm) ...
+140*du1(mm)* u3^2+140*u1* 2*u3*du3(mm)  - 2520* 5*u1^4*du1(mm)* u2- 2520* u1^5* du2(mm) + 720* 7*u1^6*du1(mm);

%dk7(mm)= du7(mm)- 7 *du1(mm)* u6 + 42* u1^2* u5 - 21* u2 *u5- 210* u1^3 *u4+ ...
% 210* u1*u2* u4 - 35* u3 *u4 + 840 *u1^4* u3 + 210* u2^2* u3 + ...
% 210*12* u1^3* u2^2   - 1260* u1^2 * u2* u3  - 630 *u1* u2^3 + 140*u1* u3^2  - 2520* u1^5* u2 + 720* u1^7

%dk7(mm)= du7(mm) - du1(mm)* u6- u1* du6(mm) - 6* dk2(mm)* u5 - 6 *k2 *du5(mm) - 15* dk3(mm)* u4- 15* k3* du4(mm) ...
%    - 20* dk4(mm)* u3- 20* k4* du3(mm) -15* dk5(mm)* u2-15* k5* du2(mm) - 6* dk6(mm)* u1- 6 *k6* du1(mm);

%str=input('Look at comparison of derivatives of 7th cumulant');

end
if(NMoments>=8)

dk8(mm)= du8(mm) - du1(mm)* u7- u1* du7(mm) - 7* dk2(mm)* u6- 7* k2* du6(mm) - 21* dk3(mm)* u5- 21* k3* du5(mm) ...
-35* dk4(mm)* u4-35* k4* du4(mm) - 35* dk5(mm)* u3- 35* k5* du3(mm) - 21* dk6(mm)* u2- 21* k6 *du2(mm) ...
- 7* dk7(mm)* u1- 7 *k7* du1(mm);
end
end

F(1,1)=k1-C(1);
F(2,1)=k2-C(2);
F(3,1)=k3-C(3);
F(4,1)=k4-C(4);
if(NMoments>=5)
F(5,1)=k5-C(5);
end
if(NMoments>=6)
F(6,1)=k6-C(6);
end
if(NMoments>=7)
F(7,1)=k7-C(7);
end
if(NMoments>=8)
F(8,1)=k8-C(8);
end

for mm=1:SeriesOrder+1
dF(1,mm)=dk1(mm);
dF(2,mm)=dk2(mm);
dF(3,mm)=dk3(mm);
dF(4,mm)=dk4(mm);
if(NMoments>=5)
dF(5,mm)=dk5(mm);
end
if(NMoments>=6)
dF(6,mm)=dk6(mm);
end
if(NMoments>=7)
dF(7,mm)=dk7(mm);
end
if(NMoments>=8)
dF(8,mm)=dk8(mm);
end
end

end


.
.
.
function [Moments] = CalculateMomentsOfZSeries(a0,a,SeriesOrder,NMoments)

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

EZ(1)=0;
EZ(2)=1;
LogEZ(2)=0;
for nn=3:NMoments*SeriesOrder
if rem(nn,2)==1
EZ(nn)=0;
else
EZ(nn)=EZ(nn-2)*(nn-1);
LogEZ(nn)=log(EZ(nn-2))+log(nn-1);
%EZ(nn);
end
end
%EZ
%EZ(1:30)
%str=input('Look at numbers')
a(SeriesOrder+1:SeriesOrder*NMoments+1)=0;
b0=a0;
b=a;

for mm=1:NMoments
if(mm>1)
%[b0,b] =SeriesProductLogarithmic(a0,a,b0,b,SeriesOrder*mm);

[b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*mm);

%b0-bb0

%b-bb

%str=input('Look at two products')

b(SeriesOrder*mm+1:SeriesOrder*NMoments+1)=0;
end
Logb=log(abs(b));
Signb=sign(b);
EXZ(mm,1)=b0;
for pp2=1:SeriesOrder*mm
if rem(pp2,2)==0
%EXZ(mm,1)=EXZ(mm,1)+exp(Logb(pp2)+LogEZ(pp2)).*Signb(pp2);
EXZ(mm,1)=EXZ(mm,1)+b(pp2).*EZ(pp2);
%           b(pp2).*EZ(pp2)
%           exp(Logb(pp2)+LogEZ(pp2)).*Signb(pp2)
%           mm
%           pp2
%str=input('Look at moment values')
end
end
end

for nn=1:NMoments
Moments(nn)=EXZ(nn,1);
end

end


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

Amin
Topic Author
Posts: 1880
Joined: July 14th, 2002, 3:00 am

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

Friends, In another day or two, I want to post a matlab program that takes moments from non-negative positive skewed densities (like densities of lognormal and CEV SDEs and stoch vol asset SDEs) and convert those moments of positive skewed densities into moments of normal or bessel type densities. So once we have specified moments of our skewed and positive densities, we would choose a  suitable transformation and convert moments of non-negative positively skewed densities into moments of bessel type densities. Later we would find a Z-series representation that matches the moments of bessel type density and then transform the Bessel Z-series into Z-series of original non-negative skewed density.
Though it might be interesting, I do not think there is any special need for dealing directly with gamma based densities. We can always transform the moments of non-negative skewed distributions into moments of bessel type densities and then do all the Z-series calculations and then transform the  newly calculated  Z-series into  Z-series of distribution in original coordinates.
After preparing all this infrastructure, we can do away with SDEs and possibly deal directly with the Z-series densities and define their evolution in a way that matches the evolution of the options, derivatives and other contingent claims. Not only can we define variances, we can specify higher moments of evolution and all of this while staying calibrated to implied volatility surface.
We can easily calculate all greeks. There is usually some difficulty in finding vega (sensitivity to implied volatility) in SDE models, but while doing standardization, we have seen that Z-series scales directly with variance (implied variance in case of derivatives) so imagine calibrating a volatility surface with Z-series and then slightly change the variance standardization coefficient and rescale the density and calculate the sensitivity of option prices to change in implied variance which can easily be converted to sensitivity in option prices to implied volatility which is our vega. All we have to do is to perturb the standardization of the fitted implied distribution slightly and re-calculate the option prices.
If we have calculated Z-series representation of an SDE based distribution, we can also use it in SDE framework to calculate vega and other greeks with respect to implied volatility just by rescaling the Z-series as I already mentioned.
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

Amin
Topic Author
Posts: 1880
Joined: July 14th, 2002, 3:00 am

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

Friends, in my experiments I noticed that this new function to improve the Newton guess is slightly better than the previous function. Please play around with this new function for Newton seed and replace the previous function if you find that this new one is better.
I could not work yesterday but hopefully my program for finding moments of transformed densities from moments of densities in original coordinates would be ready by mid week.
I am posting the new function to improve the initial Newton guess below. Replace the original function  "ImproveNewtonGuess02" in last matlab program with this function  "ImproveNewtonGuess08"
.
.
.
function [cc0Best,ccBest] = ImproveNewtonGuess08(cc0,cc,rmu,Iterations)
%This function finds a good guess for Newton method by heuristic bisection.

SeriesOrder=7;
NMoments=8;
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments);

%ObjectiveBest=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
%    abs((rmu(5)-Moments(5))^(1/2.5))+abs((rmu(6)-Moments(6))^(1/3.0))+abs((rmu(7)-Moments(7))^(1/3.5))+ ...
%    abs((rmu(8)-Moments(8))^(1/4.0));

ObjectiveBest=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
abs((rmu(5)-Moments(5))^(1/2.0))+abs((rmu(6)-Moments(6))^(1/2.0))+abs((rmu(7)-Moments(7))^(1/2.0))+ ...
abs((rmu(8)-Moments(8))^(1/2.0));

ccBest=cc;
cc0Best=cc0;
r1=.5;%bisection weight one
r2=.5;%bisection weight two

for nn=1:Iterations

[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments);
Ratio2=(r1*rmu(2)+r2*Moments(2))/Moments(2);
cc(1)=cc(1)*sqrt(Ratio2);

[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments);

Ratio3=(r1*rmu(3)+r2*Moments(3))/Moments(3);
cc(2)=cc(2)*(Ratio3);

cc0=-(cc(2)+3*cc(4)+15*cc(6));

[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments);

M4Target=rmu(4)-(-3*Moments(2).^2);
M4=Moments(4)-(-3*Moments(2).^2);

Ratio4=(r1*M4Target+r2*M4)/M4;
cc(3)=cc(3)*(Ratio4);

[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments);

%k5=u5-5*u4*u1-10*u3*u2+20*u3*u1^2+30*u2^2*u1-60*u2*u1^3+24*u1^5;
M5Target=rmu(5)-(- 10* Moments(3) .*Moments(2) );
M5=Moments(5)-(- 10* Moments(3) .*Moments(2) );

Ratio5=(.5*M5Target+.5*M5)/M5;
%Ratio5=(r1*rmu(5)+r2*Moments(5))/Moments(5);
cc(4)=cc(4)*(Ratio5);

cc0=-(cc(2)+3*cc(4)+15*cc(6));
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments);
%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;
%u6-15*u4*u2-10*u3^2+30*u2^3
M6Target=rmu(6)-(- 15* Moments(4) .*Moments(2)- 10* Moments(3).^2 + 30* Moments(2).^3 );
M6=Moments(6)-(- 15* Moments(4) .*Moments(2)- 10* Moments(3).^2 + 30* Moments(2).^3 );

Ratio6=(.5*M6Target+.5*M6)/M6;
%Ratio6=(r1*rmu(6)+r2*Moments(6))/Moments(6);
cc(5)=cc(5)*(Ratio6);

%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 ;
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments);
M7Target=rmu(7)-(- 21* Moments(2) .*Moments(5)- 35* Moments(3) *Moments(4) + 210* Moments(2).^2* Moments(3) );
M7=Moments(7)-(- 21* Moments(2) .*Moments(5)- 35* Moments(3) *Moments(4) + 210* Moments(2).^2* Moments(3) );

Ratio7=(.5*M7Target+.5*M7)/M7;
%Ratio7=(.5*rmu(7)+.5*Moments(7))/Moments(7);
%if(abs(Moments(7))>abs(rmu(7)))
cc(6)=cc(6)*(Ratio7);
%end
cc0=-(cc(2)+3*cc(4)+15*cc(6));
%cc0=-(cc(2)+3*cc(4)+15*cc(6));
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments);

Ratio8=(r1*rmu(8)+r2*Moments(8))/Moments(8);
cc(7)=cc(7)*(Ratio8);

cc0=-(cc(2)+3*cc(4)+15*cc(6));
[Moments] = CalculateMomentsOfZSeries(cc0,cc,SeriesOrder,NMoments);

[IsValidFlag] = CheckIsValidDensity(cc0,cc)

%ObjectiveNew=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
%    abs((rmu(5)-Moments(5))^(1/2.5))+abs((rmu(6)-Moments(6))^(1/3.0))+abs((rmu(7)-Moments(7))^(1/3.5))+ ...
%    abs((rmu(8)-Moments(8))^(1/4.0));

ObjectiveNew=100000*(abs(rmu(1)-Moments(1)))+abs(rmu(2)-Moments(2))+abs((rmu(3)-Moments(3))^(1/1.5))+abs((rmu(4)-Moments(4))^(1/2.0)) + ...
abs((rmu(5)-Moments(5))^(1/2.0))+abs((rmu(6)-Moments(6))^(1/2.0))+abs((rmu(7)-Moments(7))^(1/2.0))+ ...
abs((rmu(8)-Moments(8))^(1/2.0));

if((ObjectiveNew<ObjectiveBest) &&( IsValidFlag))
ObjectiveBest=ObjectiveNew;

ccBest=cc;
cc0Best=cc0;
end

end

[Moments] = CalculateMomentsOfZSeries(cc0Best,ccBest,SeriesOrder,NMoments)
rmu
str=input('Look at moments--In function---After');

end


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

Amin
Topic Author
Posts: 1880
Joined: July 14th, 2002, 3:00 am

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

Friends, there are two ways to transform the moments from original coordinates to transformed coordinates.
One way is to construct the density in original coordinates, transform the density and calculate moments. This is very cheaply done (as opposed to what friends might think) if original density is Bessel/normal based density and transformed density is CEV or lognormal density.
However converting moments from CEV or lognormal based densities to bessel densities might not be very simple since CEV or lognormal based density has to be constructed first.
Second way is to use a "Taylor series for the transformation function" around the mean in original coordinates and take expectations. Expansion polynomials become central moments of respective order after taking expectations. This works well for the mean (of the transformation) and first few moments but higher moments do poorly. I did this with eight terms/moments in the Taylor expansion.
But when we are converting from Bessel/normal based coordinates, even though our Z-series might have only eight elements, we can calculate higher central moments than first eight moments and they can be added to expectations to improve accuracy. But I am afraid there are some accuracy issues with powers of standard normal when we are raising standard normal to 80th power but I remember that I earlier wrote a log-based Z-series moment calculation function in which powers of standard normal are calculated in logarithmic coordinates and since these high powers have extremely small coefficients which are again calculated in logarithmic form and these are added (to find logarithm of the product) and usually the resulting logarithm is a reasonably small number in original coordinates which can be retrieved by taking exponential of the product. I will try seeing if it improves accuracy.
Anyway, I will post the new matlab program tomorrow with a lot of basic routines to play around with transformation of series and their moments and I hope friends will find the basic routines useful. I will be posting the new program 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

Amin
Topic Author
Posts: 1880
Joined: July 14th, 2002, 3:00 am

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

When we were working with SDEs in context of Z-series, we learnt that in order to get probability density of functions of SDE variable, we have to Taylor expand the function around the median of Z-series of SDE variable and then use the coefficients of Z-series in Taylor derivatives calculations and it worked very well.
However, in order to find the expectations of functions of the SDE/stochastic variable, we need to Taylor expand the function of stochastic variable around mean of the stochastic variable and then replace the expansion polynomial powers with central moments of that order and find expectation of the functions of the stochastic variable. I am sure many friends have used something like this before but this is very convenient in our framework.
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

Amin
Topic Author
Posts: 1880
Joined: July 14th, 2002, 3:00 am

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

But I am afraid there are some accuracy issues with powers of standard normal when we are raising standard normal to 80th power but I remember that I earlier wrote a log-based Z-series moment calculation function in which powers of standard normal are calculated in logarithmic coordinates and since these high powers have extremely small coefficients which are again calculated in logarithmic form and these are added (to find logarithm of the product) and usually the resulting logarithm is a reasonably small number in original coordinates which can be retrieved by taking exponential of the product. I will try seeing if it improves accuracy.
.
In a day or two, I will also upload this program that friends can use when they want higher accuracy for calculation of moments of Z-series.
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

Amin
Topic Author
Posts: 1880
Joined: July 14th, 2002, 3:00 am

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

Friends, my work is going alright and I hope to be posting new matlab programs later today or tomorrow. But I am writing this post to mention the new persecution and torture activity from mind control agents. For past few days, I have been working in the drawing room of my house with AC on. Usually there would be very small amount of gas in the air-conditioners but it would be generally tolerable and I could continue to work. But today, just a bit ago, they increased gas remarkably and it became difficult to breathe properly and I had to turn off the AC and open the room doors. With my breathing, the gas would settle on inside of my mouth and I would feel severe constriction of muscles inside my mouth around the jaws.
All day yesterday, they continued to charge my back very severely and there would be extreme itch in my back and I would continue to request the mind control agents to decrease the intensity of EM waves but to no avail.
They also asked my mother yesterday to add drugs to food and  she drugged the cooked food in our house and I could not work later yesterday night after having food.
Tomorrow I will get an antipsychotic injection and it would become difficult for me to work for several days. Though I do feel slightly better (as compared to earlier) after the doctor increased the period of injections from two weeks to three weeks but I still do not feel perfectly good and want to request friends to please ask the mind control agencies to end my antipsychotics.
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

Amin
Topic Author
Posts: 1880
Joined: July 14th, 2002, 3:00 am

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

Friends, sorry for this delay but it will take me another two or three days to complete the work. I want to complete writing a bit more than twenty small matlab routines to find for example, CDF of the Z-series densities, transformation of the densities and finding densities of functions of stochastic variables when density of original variable has been found in Z-series, finding expectations of functions of the stochastic variable from its moments, finding option prices analytically for a Z-series density or for transformation of the Z-series density, finding some greeks of the option prices, sensitivities of option prices to coefficients of Z-series etc. I want to write small independent matlab routines that friends can easily transfer to their programming languages of choice.
Similarly they add mind control chemicals to tea if I can find tea of some less known brand at a store, they approach the store later and drug the tea there. There is no question of getting good coffee and I have given up trying to get good coffee and I only try to get good tea of some lesser known brand that still has not started to add mind control drugs to their product at manufacturing source.
CIA pays every Pakistan army chief several tens of millions of dollars every year. This is the only reason mind control has continued and has increased in Pakistan for more than past twenty years. Army persuades or forces manufacturers of coffee, tea, beverages, medicines, soaps, lotions, toothpastes and a large number of other products to add mind control drugs to their products at manufacturing source. Importers are forced to add mind control drugs to their imported products before selling it to local distributors. All of this is on top of the accepted practice in which army soldiers or agents go from shop to shop in entire markets to drug food and beverages they want to drug in order to target some individual.
At this point, I want to tell friends even though Pakistan army does all the work to drug the food and lobby for mind control, it is American agents who use the mind control equipment and retard intelligent and ambitious people. Generals of Pakistani army work very closely with American mind control agents since tens of millions of dollars are paid to them in return of their cooperation. I really believe that Pakistan army generals have received bribes worth 5-8 billion dollars in past twenty years but have dented our country of more than a hundred and fifty billion dollars by retarding several hundreds of intelligent enterprising young brains.
I want to tell Americans that most educated Pakistanis understand that Americans are mostly very nice people who want development in our poor countries and do not want to see our people suffer in poverty  but American army(which is absolutely not representative of American people) is run by neocons, racist conservatives, religious conservatives and other religious vested interests and they actively thwart the attempts of our nations to succeed economically.  When people say bad things about America and shout anti-American slogans, they really do not realize that it is mostly American army run by vested interests that controls diplomacy and does all machinations in our country. Otherwise most of  American citizens do not want anything bad to happen to our country and would be willing to help us if they somehow can.
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

Amin
Topic Author
Posts: 1880
Joined: July 14th, 2002, 3:00 am

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

Friends, a few hours ago, I received antipsychotic fluanxol depot shot in my arm. Since I had sensed that things were getting heated for past few days, I was very worried about getting a good injection. I knew if they would anticipate where I would buy the injection, they would try to make me buy a special injection that has mind control drugs in it. I am sure all the fluanxol injections in the market have special mind control drugs in them but these injection remain on the shelves in pharmacies for more than six months to a year. And mind control agents mostly want to add drugs specifically for the current condition and mental state of target's neurotransmitters at the present time. So I was very worried that they would try to sell me drugged injections today since mind control agencies were becoming bitter to control me.
I left the house on my car to buy the injections. I tried to buy them at a FDFP Fazal Din Pharma Plus in model town but they did not have the injections. I went to ferozpur road from model town and drove towards DHA/Walton side. In DHA, I stopped at a large Mehmood Pharmacy in front of packages mall. I knew they used to have injections and I had earlier bought Fluanxol injections from there.
Here I want to tell friends that mind control agents can project foreign thought and voice to skull on people remotely and Pakistan army has approached many pharmacies in posh areas of the city and has told them to abide by the instructions when thought is projected on them. They are already told in advance to follow the instructions when they would hear a foreign sound in their brain.
So I was at this Mehmood pharmacy where I asked the salesperson to sell me the injections. The guy waited for a bit and told me they did not have the injections. I was slightly suspicious and asked another salesman right next and he cut the first salesman saying, "He is talking about depot injections. We have them." I asked them to look for the injections but just in another twenty seconds the other salesman (without any special effort at finding the injections) also said that they did not have the injections. I was very suspicious but could not do anything and had to leave. But I did realize that mind control agents might be trying to remotely ask them to not sell the injections and this meant that they wanted me to buy injections from their choice of pharmacy.
I was a bit upset at this time and I drove into Cantt. and went to this round market. I forgot the name of this market but I recall they used to have a petrol pump and then a large pizza hut in other side of the round traffic circle. A road also comes directly to this market from sherpao bridge from Gulberg side. I entered CSH pharmacy there and bought the injection. But the salesman  gave me a very fancy ice bag made of plastic that was possibly worth more than the injections I had bought. He told me that he had run out of other ordinary ice bags and had to give me this special ice bag. I was not particularly suspicious but still I kept thinking why would he give me this special ice bag.
But I was very sure that this last injection was good and I dropped the previous injection bought from CSH. The injection is cheap in Pakistan and costs merely 2.5 dollars worth of Pakistani rupees.
So I was able to get the good injection but the bad part of the story is that things are getting heated and mind control agencies are trying to somehow control me again. I am very afraid that they would try to inject me at night while sleeping if I ever sleep outside of my locked room. I am worried since sometimes I sleep for a little (an hour or half) when working at night in the drawing room as I would basically fall asleep while working. And our drawing room cannot be properly locked from inside. I am very afraid that mind control agents are looking for some opportunity to somehow retard me again for a few weeks.
I also want to make this request to Pakistani friends and also foreign embassies in Pakistan to keep a vigil since I am very afraid that Pakistan army would start drugging all of the fluanxol injections all over 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

Amin
Topic Author
Posts: 1880
Joined: July 14th, 2002, 3:00 am

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

Friends,  my  small program to improve the Newton guess worked reasonably well but I was convinced that there could be better methods to improve the initial Newton guess. I continued to work on another method to improve the Newton guess last night and have had substantial success. I have ideas about at least  two different new methods to improve Newton guess and I want to work with these ideas for a few days. I really believe that a better method for initial Newton guess would really help the use of Z-series everywhere and it is more important than other routine things. These methods are more intuitive and I am sure that friends would easily be able to make changes and improve them.
It will take another two to three days for the effects of injection to catch up with me and I want to finish this work before that.
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

Amin
Topic Author
Posts: 1880
Joined: July 14th, 2002, 3:00 am

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

Friends, I have been able to work with a new method to improve the initial Newton guess. It works better than the previous method and is totally intuitive. Everything in the method is perfectly explainable as opposed to what we had in the previous method for Newton guess. So I am sure that experts would be able to take a good basic framework and make their own modifications to further improve the method. I am working to check the method a bit more and would post the new routine later today at night or possibly tomorrow with a very detailed explanation of the method.
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

Amin
Topic Author
Posts: 1880
Joined: July 14th, 2002, 3:00 am

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

Friends, even though the prototype is working well, I am trying to make my program for finding initial guess for Newton method a good program. And the final program may take a few days. At this point I think I would be able to post a good program on weekend. I really hope to post a good program. I also think there may not be any need for initial construction of polynomial densities after this new program. I will keep friends updated.
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