Serving the Quantitative Finance Community

 
User avatar
Amin
Topic Author
Posts: 2012
Joined: July 14th, 2002, 3:00 am

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

November 28th, 2022, 3:02 pm

Cuch, very sorry for delay in my reply. Below is the Mathematica code to solve any 1st order or nth order ODE. I am writing the code in text format below and explain it again in a next code window.
.
Clear[t,y,Zi,Wi,Yi,Y,y0,i,i1,p,n,ZAns];
n=1;
Array[y,n+1,0];
y[1]:= y[0]^2;
Array[y0,n,0];
y0[0]=1;
p=8;
Y:=y0[0];
For[k=1,k<n,k++,Y=Y+y0[k] *t^k/k!];
Array[Zi,n+p-1,0];
Array[Wi,n+p-1,0];
Array[Yi,n+p-1,0];
Array[ti,n+p-1,0];
Zi[n]=(y[n]/.t-> ti[n]);
For [i=n,i<p+n-1,i++,(Zi[i+1]=0;For [k=0,k<n,k++,(Zi[i+1]=Zi[i+1]+D[Zi[i],y[k]]*(y[k+1]/.t-> ti[i+1]));]);];
For[i=n,i<= p+n-1,i++,(Wi[i]=Zi[i];For[k=0,k<n,k++,(Wi[i]=(Wi[i]/.y[k]-> y0[k]))];)];
For[i=n,i<= p+n-1,i++,(Yi[i]=Wi[i];For[i1=i,i1>=1,i1--,(Yi[i1-1]=Integrate[Yi[i1],{ti[i1],0,ti[i1-1]}])];Y=Y+(Yi[0]/.ti[0]-> t);)]
ZAns=Collect[Y,t,Simplify]//PolynomialForm[#,TraditionalOrder -> False]&
.
In Below I explain the above code with some comments.
Clear[t,y,Zi,Wi,Yi,Y,y0,i,i1,p,n,ZAns];   This line is handy to clear variables especially when you want to solve for more than one ODE.
n=1;                                      n is order of the ODE. For 1st order ODE n=1; 
Array[y,n+1,0];                           Define a workhorse variable.
y[1]:= y[0]^2;                            Specify the ODE. y[1]= dy/dt. y[0]=y[t] and similarly for higher order ODE's y[2]=d2y/dt2
Array[y0,n,0];
y0[0]=1;                                  Specify initial conditions. y0[0]=y(0) similalry for 2nd order ODE, you will also have to specify y0[1]=dy(0)/dt and so on.
p=8;                                       Number of powers in expansion of series.
Y:=y0[0];                                  Allocate first initial condition.
For[k=1,k<n,k++,Y=Y+y0[k] *t^k/k!];       Allocation of higher order initial conditions.
Array[Zi,n+p-1,0];
Array[Wi,n+p-1,0];
Array[Yi,n+p-1,0];
Array[ti,n+p-1,0];
Zi[n]=(y[n]/.t-> ti[n]);
For [i=n,i<p+n-1,i++,(Zi[i+1]=0;For [k=0,k<n,k++,(Zi[i+1]=Zi[i+1]+D[Zi[i],y[k]]*(y[k+1]/.t-> ti[i+1]));]);];
For[i=n,i<= p+n-1,i++,(Wi[i]=Zi[i];For[k=0,k<n,k++,(Wi[i]=(Wi[i]/.y[k]-> y0[k]))];)];
For[i=n,i<= p+n-1,i++,(Yi[i]=Wi[i];For[i1=i,i1>=1,i1--,(Yi[i1-1]=Integrate[Yi[i1],{ti[i1],0,ti[i1-1]}])];Y=Y+(Yi[0]/.ti[0]-> t);)];     This is integration loop which iteratively integrates starting from zero. replace zero with initial time when it is different from zero.
ZAns=Collect[Y,t,Simplify]//PolynomialForm[#,TraditionalOrder -> False]&
.
I think your ODE(2) can be very easily solved. ODE(1) seems singular at t=0 but can possibly be solved by changing starting value to some positive time and specifying initial conditions there.
Please feel free to use the mathematica ODE solver above and come back with your comments.
Sorry, I have been pretty busy with work on my trading algorithms and want to apply them on NASDAQ stocks so little time left for other things. Will come back soon to tackle the discussion.
In my paper, I have given a mathematica routine that can solve arbitrary systems of nth order ODEs (at least in theory). I have also given an example explaining the algorithm where I used the algorithm to solve a system of three 2nd order non-linear ODEs. 
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
 
User avatar
Cuchulainn
Posts: 18230
Joined: July 16th, 2004, 7:38 am
Location: Llanfairpwllgwyngyllgogerychwyrndrobwllllantysiliogogogoch

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

November 29th, 2022, 11:23 am

Thanks, Amin
Code is important but it is not the driver. It just produces a number. 
We are not on the same wavelength, unfortunately. My two examples were concrete insights into the bespoke qualitative behaviour of ODEs.

„A good technical writer, trying not to be obvious about it, but says everything twice: formally and informally. Or maybe three times.“ —  Donald Ervin Knuth
They sought it with thimbles, they sought it with care; They pursued it with forks and hope; They threatened its life with a railway-share; They charmed it with smiles and soap.
 
User avatar
Cuchulainn
Posts: 18230
Joined: July 16th, 2004, 7:38 am
Location: Llanfairpwllgwyngyllgogerychwyrndrobwllllantysiliogogogoch

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

November 29th, 2022, 11:34 am

For completeness

[$]dx/dt = x^2, x(0) = 1[$] general solution [$]x = 1/(1-t)[$] (B)

. The solution can be found by inspection or Picard iteration serie 
. The solution does blow up afer t > 1
. Solving (B) by ODE solvers, including yours ... a prioiri quantitatve results on stability and accurscy

These are the compelling issues in numerical analysis, like gravity in physics.
They sought it with thimbles, they sought it with care; They pursued it with forks and hope; They threatened its life with a railway-share; They charmed it with smiles and soap.
 
User avatar
Amin
Topic Author
Posts: 2012
Joined: July 14th, 2002, 3:00 am

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

November 30th, 2022, 3:21 pm

Sorry Cuch, a lot of other things to do. I will come back to your examples soon when I will have some more time. I am sure they are very simple.

Friends, I have decided to write 3-5 formal papers about the material in this thread that would ideally be submitted to Wilmott magazine and I will soon be approaching Wilmott people for their requirements and guidelines. 
The research subject of first paper would be finding the dynamics of variance/volatility SDE so that volatility surface is closely fitted across time given that asset SDE in two dimensional correlated SDE system is given either by lognormal or CEV dynamics.
All through this paper, we will remain mainly in a single dimensional SDE framework other than just a few small excursions into two dimensions.
For my first paper, I am cleaning up my work on one dimensional SDEs with addition of cumulants. And also on time integral of the mean-reverting and other SDEs. I want to solve for distribution of integral of volatility/variance SDE by backing it out from option prices using lognormal or CEV formula. And then I will find the dynamics of volatility SDE so that it perfectly fits the implied distribution of integral of variance. While backing out integral of volatility, I will also have to do accounting for correlation and contribution of correlated part of SV SDE but this should be simple given the insights we have learnt about variances and their addition within Z-series framework. All of this will be done in a one dimensional SDE framework. We have already done 70% of the work in this thread. I will be changing a few things like I do not want to go into cumulants at all (even though they should be used a lot for inference) and will be using new alternative methods. I will also be making some changes to calculation of time integrals of SDEs. Only major part that remains to be done is bootstrapping of integral of variance SDE and correlated part of the SDE from the cross-section of option prices. Since we want to remain in the one dimensional domain, I really hope that computational time for option chains out to 4-5 years will be in seconds and not in minutes. 
And my time frame to complete the unfinished part is a few weeks. I will be posting the code on Wilmott as usual. 
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
 
User avatar
Cuchulainn
Posts: 18230
Joined: July 16th, 2004, 7:38 am
Location: Llanfairpwllgwyngyllgogerychwyrndrobwllllantysiliogogogoch

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

December 1st, 2022, 10:37 am

Sorry Cuch, a lot of other things to do. I will come back to your examples soon when I will have some more time. I am sure they are very simple.
Take your time. Don't want be a party pooper, but the mathematical and numerical  issues are far from simple.
They sought it with thimbles, they sought it with care; They pursued it with forks and hope; They threatened its life with a railway-share; They charmed it with smiles and soap.
 
User avatar
Cuchulainn
Posts: 18230
Joined: July 16th, 2004, 7:38 am
Location: Llanfairpwllgwyngyllgogerychwyrndrobwllllantysiliogogogoch

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

December 1st, 2022, 11:27 am

Cuch, very sorry for delay in my reply. Below is the Mathematica code to solve any 1st order or nth order ODE. I am writing the code in text format below and explain it again in a next code window.
.
Clear[t,y,Zi,Wi,Yi,Y,y0,i,i1,p,n,ZAns];
n=1;
Array[y,n+1,0];
y[1]:= y[0]^2;
Array[y0,n,0];
y0[0]=1;
p=8;
Y:=y0[0];
For[k=1,k<n,k++,Y=Y+y0[k] *t^k/k!];
Array[Zi,n+p-1,0];
Array[Wi,n+p-1,0];
Array[Yi,n+p-1,0];
Array[ti,n+p-1,0];
Zi[n]=(y[n]/.t-> ti[n]);
For [i=n,i<p+n-1,i++,(Zi[i+1]=0;For [k=0,k<n,k++,(Zi[i+1]=Zi[i+1]+D[Zi[i],y[k]]*(y[k+1]/.t-> ti[i+1]));]);];
For[i=n,i<= p+n-1,i++,(Wi[i]=Zi[i];For[k=0,k<n,k++,(Wi[i]=(Wi[i]/.y[k]-> y0[k]))];)];
For[i=n,i<= p+n-1,i++,(Yi[i]=Wi[i];For[i1=i,i1>=1,i1--,(Yi[i1-1]=Integrate[Yi[i1],{ti[i1],0,ti[i1-1]}])];Y=Y+(Yi[0]/.ti[0]-> t);)]
ZAns=Collect[Y,t,Simplify]//PolynomialForm[#,TraditionalOrder -> False]&
.
In Below I explain the above code with some comments.
Clear[t,y,Zi,Wi,Yi,Y,y0,i,i1,p,n,ZAns];   This line is handy to clear variables especially when you want to solve for more than one ODE.
n=1;                                      n is order of the ODE. For 1st order ODE n=1; 
Array[y,n+1,0];                           Define a workhorse variable.
y[1]:= y[0]^2;                            Specify the ODE. y[1]= dy/dt. y[0]=y[t] and similarly for higher order ODE's y[2]=d2y/dt2
Array[y0,n,0];
y0[0]=1;                                  Specify initial conditions. y0[0]=y(0) similalry for 2nd order ODE, you will also have to specify y0[1]=dy(0)/dt and so on.
p=8;                                       Number of powers in expansion of series.
Y:=y0[0];                                  Allocate first initial condition.
For[k=1,k<n,k++,Y=Y+y0[k] *t^k/k!];       Allocation of higher order initial conditions.
Array[Zi,n+p-1,0];
Array[Wi,n+p-1,0];
Array[Yi,n+p-1,0];
Array[ti,n+p-1,0];
Zi[n]=(y[n]/.t-> ti[n]);
For [i=n,i<p+n-1,i++,(Zi[i+1]=0;For [k=0,k<n,k++,(Zi[i+1]=Zi[i+1]+D[Zi[i],y[k]]*(y[k+1]/.t-> ti[i+1]));]);];
For[i=n,i<= p+n-1,i++,(Wi[i]=Zi[i];For[k=0,k<n,k++,(Wi[i]=(Wi[i]/.y[k]-> y0[k]))];)];
For[i=n,i<= p+n-1,i++,(Yi[i]=Wi[i];For[i1=i,i1>=1,i1--,(Yi[i1-1]=Integrate[Yi[i1],{ti[i1],0,ti[i1-1]}])];Y=Y+(Yi[0]/.ti[0]-> t);)];     This is integration loop which iteratively integrates starting from zero. replace zero with initial time when it is different from zero.
ZAns=Collect[Y,t,Simplify]//PolynomialForm[#,TraditionalOrder -> False]&
.
I think your ODE(2) can be very easily solved. ODE(1) seems singular at t=0 but can possibly be solved by changing starting value to some positive time and specifying initial conditions there.
Please feel free to use the mathematica ODE solver above and come back with your comments.
Sorry, I have been pretty busy with work on my trading algorithms and want to apply them on NASDAQ stocks so little time left for other things. Will come back soon to tackle the discussion.
In my paper, I have given a mathematica routine that can solve arbitrary systems of nth order ODEs (at least in theory). I have also given an example explaining the algorithm where I used the algorithm to solve a system of three 2nd order non-linear ODEs. 
I have non-hands on experience of Mathematica from a previous article with Alan and Paul.
AFAIR the routine "Integrate" is a magic black box ODE solver and has hidden know-how on "tricky" ODEs. It does noty bring you closer to an understanding of the essential complexity. It is also slow when used as a method of lines (MOL) , and relatively few quants use it (?). For speed and insights, I reckon hand-crafted C++ code is needed. Quants  want to know what is happening under the radar.
They sought it with thimbles, they sought it with care; They pursued it with forks and hope; They threatened its life with a railway-share; They charmed it with smiles and soap.
 
User avatar
Amin
Topic Author
Posts: 2012
Joined: July 14th, 2002, 3:00 am

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

December 4th, 2022, 3:37 pm

Cuch, Integrate command in my code is nothing more than symbolic integration. It gives an analytic answer which is very helpful. It is not at all what you suspect and nothing more than symbolic computing. The above code solves nth order ODEs. So in order to automate everything, there are loops that solve iterated integrals but only symbolic integration is performed in the final line with loops.
Another thing when we specify ODE
y[1]:= y[0]^2; in the code means dy/dt=y(t)^2;
you could have specified an ODE with t in it as

y[1]:=t y[0]^2 -t y[0]; which would mean the ODE is: dy/dt=t y(t)^2 - t y(t);

It is this specified ODE that is symbolically integrated in iterated integrals using loops. There are loops both to tackle change of variables and also iterated symbolic integrations. The final answer of mathematica code is a symbolic series expansion of the true solution of the ODE. It is not anything like numerical integration or MOL or automatic ODE solutions etc. as you fear.
Mathematic code is a wizard. Paste it in mathematica with your own ODE and initial conditions and you get the series expansion solution right away.


Friends, I went back to my work on solution of SDEs with addition of cumulants. I had earlier solved one dimensional SDEs with this method but Z-series I used was only up till 5th power of Z. The previous solution also had relatively poor convergence because root-search was done over cumulants. With only five powers in the Z-series solution done in Bessel coordinates meant when we converted the solution into Z-series coefficients in original coordinates, the approximation in tails was poor at best and it could be easily seen that Z-series coefficients in original coordinates diverged from true solution to quite large extent in both tails.

Now, I have re-done the above solution of SDEs with a seven power Z-series. I have also updated the root-search and changed that from cumulants to moments which is very helpful. Z-series root-search solution almost always converges(though there are relatively rare exceptions when cumulants/moments change too fast and root-search solution to moments does not converge but I know that can be tackled by adaptively decreasing the step size when moments/cumulants change too fast. I have not introduced this adaptive step size setting in present program but want to do it sometime in the future soon.)
But I am absolutely in love with the new program when I noticed that converting Z-series SDE coefficients in Bessel coordinates to Z-series coefficients in original coordinates very faithfully mimics the true solution. When solution is plotted from original coordinates Z-series, right tail is almost always indistinguishable from the true right tail till the end. Left tail close to zero sometimes sneaks into negative territory but that is mostly very little and very insignificant and very close to true solution.
This is very important since all the stochastic volatility analytics are done on Z-series representation of SDE in original coordinates. And if conversion of Z-series coefficients from Bessel coordinates to original coordinates is not faithful, all analytics will be very approximate which is luckily not the case when Z-series coefficients are taken to seventh power of Z in Bessel coordinates to solve the SDE.
I will post my new programs for one dimensional SDEs that solve the SDE in Bessel coordinates using 7th order Z-series and also find analytic Z-series in original coordinates in another day or two latest by Tuesday. Again I am very happy that we have been able to get a faithful analytical representation of SDEs in original coordinates. 
And we will need this analytical representation of stochastic volatility SDE in original coordinates again and again when we would solve for the SV dependent asset SDE using new methods.
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
 
User avatar
Cuchulainn
Posts: 18230
Joined: July 16th, 2004, 7:38 am
Location: Llanfairpwllgwyngyllgogerychwyrndrobwllllantysiliogogogoch

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

December 4th, 2022, 4:23 pm

 It is not anything like numerical integration or MOL or automatic ODE solutions etc. as you fear.

Me fear?  :roll:
Actually, I was confused with NSolve. 
It is possible ODE solvers

https://library.wolfram.com/infocenter/ ... aPart1.pdf

And probably better performance than Shlow Symbolix.
They sought it with thimbles, they sought it with care; They pursued it with forks and hope; They threatened its life with a railway-share; They charmed it with smiles and soap.
 
User avatar
Amin
Topic Author
Posts: 2012
Joined: July 14th, 2002, 3:00 am

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

December 6th, 2022, 12:14 pm

Friends, I am posting my code for solution for densities of one dimensional SDEs from addition of cumulants and calculation of their variable representation as a Z-series which is a series in powers of standard normal with variable coefficients (which we calculate for every specific SDE).
I want to submit this work for publication in Wilmott magazine. I am still unsure whether to submit it as one large article that covers all the work up to solution of system of stochastic volatility SDEs or as in various smaller articles. But here is the code.
.
.
function [] = FPESeriesCoeffsFromMomentsAndVarianceAdditionAdaptStep()

%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334 
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=mu1 x(t)^beta1 dt + mu2 x(t)^beta2 dt +sigma x(t)^gamma dz(t)

%I have simulated the transformed 
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coordinates.


dt=.125/2/2*1;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.
Tt=128/2/2*2;%*4*4*1;%*4;%128*2*2*2;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year; 
T=Tt*dt;

%Below I have done calculations for smaller step size at start.
ds(1:3)=dt/4;

%Below order for monte carlo
OrderA=4;  %
OrderM=4;  %

if(T>=1)
dtM=1/32;%.0625/1;
TtM=T/dtM;
else
dtM=T/16;%.0625/1;
TtM=T/dtM;
end



dNn=.1/2;   % Normal density subdivisions width. would change with number of subdivisions
Nn=45*4;  % No of normal density subdivisions

x0=1;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.95;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=3.2;%.950;   %mean reversion parameter.
theta=.25;%mean reversion target
sigma0=1.0;%Volatility value


%StepSize Management. decrease the step size if MomentsIncreaseRatio is larger
    %than  MaxMomentRatio and increase the step size if 
    %MomentsIncreaseRatio is smaller than  MinMomentRatio
    %Look at the code in the body and feel free to play with stepsize
    %variables
MaxMomentRatio=1.01;
MinMomentRatio=1.002;
MaxStepSize=1/128;
MinStepSize=1/512/16;



%you can specify any general mu1 and mu2 and beta1 and beta2.
mu1=1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
%mu1=0;
%mu2=0;

alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;%This is for expansion of integrals for calculation of drift 
%and volatility coefficients
%yy(1:Nn)=x0;                
w0=x0^(1-gamma)/(1-gamma);


NnMid=((1+Nn)/2)*dNn;
Z(1:Nn)=(((1:Nn)*dNn)-NnMid);
%Z(1:Nn)=(((1:Nn)-20.5)*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);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This block is analytics for higher order monte carlo. This is explained in
%my thread on wilmott.com

sigma11(1:OrderA+1)=0;
mu11(1:OrderA+1)=0;
mu22(1:OrderA+1)=0;
sigma22(1:OrderA+1)=0;
% index 1 correponds to zero level since matlab indexing starts at one. 
sigma11(1)=1;
mu11(1)=1;
mu22(1)=1;
sigma22(1)=1;

for k=1:(OrderA+1)
    if sigma0~=0
        sigma11(k)=sigma0^(k-1);
    end
    if mu1 ~= 0
        mu11(k)=mu1^(k-1);
    end
    if mu2 ~= 0
        mu22(k)=mu2^(k-1);
    end
    if sigma0~=0
        sigma22(k)=sigma0^(2*(k-1));
    end
end
%Ft(1:TtM+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0; %General time powers on hermite polynomials
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers on coefficients of hermite polynomials.
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers for bessel transformed coordinates.

%YCoeff0 and YCoeff are coefficents for original coordinates monte carlo.
%YqCoeff0 and YqCoeff are bessel/lamperti version monte carlo.

YCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
YqCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
%Pre-compute the time and power exponent values in small multi-dimensional arrays
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma); %expand y^alpha where alpha=1;
YqCoeff = ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);%expand y^alpha1 where alpha1=(1-gamma)
YqCoeff=YqCoeff/(1-gamma); %Transformed coordinates coefficients have to be 
%further divided by (1-gamma)

for k = 0 : (OrderA)
    for m = 0:k
        l4 = k - m + 1;
        for n = 0 : m
            l3 = m - n + 1;
            for j = 0:n
                l2 = n - j + 1;
                l1 = j + 1;
                %Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
                Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
                    - (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
                
                YCoeff0(l1,l2,l3,l4) =YCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
                YqCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
            end
        end
    end
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Below analytics to calculate the Z-series of SDE at different time steps. 
%We use the principle of addition of cumulants to find the resulting
%density of SDE at each step.
ExpnOrder=7;
SeriesOrder=7;
NMoments=8;
wnStart=1;
tic
%Tx is running time.
%T is terminal time.
%tt is time index.
Tx=0;
tt=0;
while(Tx<T)
    tt=tt+1;
   % t1=(tt-1)*dt;
   % t2=(tt)*dt;
    if(tt==1)
 
        [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),8);
 
        c0=wMu0dt+w0;
        c(1)=sigma0*sqrt(ds(tt))
        c(2:10)=0.0;
        w0=c0;
        
        Tx=Tx+ds(tt);
    
    elseif(tt==2) 

    %dt=ds(tt);    
        
    %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),8);
    [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),8);
     [b0,b] = CalculateDriftbCoeffs08(wMu0dt,dwMu0dtdw,c,SeriesOrder);
     
     
     c0=c0+b0;
     c(1)=b(1)+sqrt(c(1)^2+sigma0^2*ds(tt));
     c(2)=b(2);
     c(3)=b(3);
     c(4)=b(4);
     c(5)=b(5);
     c(6)=b(6);
     c(7)=b(7);
     
     Tx=Tx+ds(tt);       
            
    else
        
    
    %below wMu0dta is drift and dwMu0dtdwa are its derivatives.
    %wMu1dt is its first derivative and dwMu1dtdw are derivatives of first
    %derivative.


     %[wMu0dta,dwMu0dtdwa,wMu1dt,dwMu1dtdw] = BesselDriftAndDerivatives08(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,10);
     
     %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),8);
     [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),8);
     [b0,b] = CalculateDriftbCoeffs08(wMu0dt,dwMu0dtdw,c,SeriesOrder);

%Below Calculate original terms(associated with first hermite and second
%hermite polynomial in volatility)and its eight derivatives    
    [wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),8);
    [wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),8);

%Below Calculate Z-Series expansions of volatility terms associated with 
%first hermite polynomial and 2nd hermite polynomial.
    [g0,g] = CalculateDriftbCoeffs08(wVol0dt,dwVol0dtdw,c,SeriesOrder);
    [h0,h] = CalculateDriftbCoeffs08(wVol2dt,dwVol2dtdw,c,SeriesOrder);
%Principal term in first hermite was not included in expansion. It is 
%added here
    g0=g0+sigma0*sqrt(ds(tt));
    
    %Below function calculates Moments and Cumulants of diffusion term of
    %the SDE
    [Moments,kk] = CalculateVolTermMoments(g0,g,h0,h,SeriesOrder,NMoments);
    
    
 %Below drift series is linearly added to SDE series to get new intermediate 
 %SDE series.     
    c0=c0+b0;
    c(1:7)=c(1:7)+b(1:7);
    
 %Below calculate cumulants and moments of intermediate SDE series   
    [CC] = CalculateCumulants(c0,c,SeriesOrder,NMoments);
    [MM1] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
    
    
    %Below Calculate Cumulants of target density by arithmetic addition of
    %first eight cumulants of diffusion term and intermediate SDE term.
    k1=CC+kk;
    
    %Below Calculate Standardized Cumulants of the target density. 
    
    k(1)=0;
    k(2)=1;
    k(3)=k1(3)/k1(2).^1.5;
    k(4)=k1(4)/k1(2).^2.0;
    k(5)=k1(5)/k1(2).^2.5;
    k(6)=k1(6)/k1(2).^3.0;
    k(7)=k1(7)/k1(2).^3.5;
    k(8)=k1(8)/k1(2).^4.0;
    
    
    %Below Calculate Standardized Moments of target density.
    rmu(1)=k(1);
    rmu(2)=k(1)*rmu(1)+k(2);
    rmu(3)=k(1)*rmu(2)+2*k(2)*rmu(1)+k(3);
    rmu(4)=k(1)*rmu(3)+3*k(2)*rmu(2)+3*k(3)*rmu(1)+k(4);    
    rmu(5)=k(1)*rmu(4)+4*k(2)*rmu(3)+6*k(3)*rmu(2)+4*k(4)*rmu(1)+k(5);
    rmu(6)=k(1)*rmu(5)+5*k(2)*rmu(4)+10*k(3)*rmu(3)+10*k(4)*rmu(2)+5*k(5)*rmu(1)+k(6);
    rmu(7)=k(1)*rmu(6)+6*k(2)*rmu(5)+15*k(3)*rmu(4)+20*k(4)*rmu(3)+15*k(5)*rmu(2)+6*k(6)*rmu(1)+k(7);
    rmu(8)=k(1)*rmu(7)+7*k(2)*rmu(6)+21*k(3)*rmu(5)+35*k(4)*rmu(4)+35*k(5)*rmu(3)+21*k(6)*rmu(2)+7*k(7)*rmu(1)+k(8);

    %We provide initial guess for calibration procedure
    
    if(tt==3)
        a0Guess=(c0-k1(1))/sqrt(CC(2));
        aGuess(1)=c(1)/sqrt(CC(2));
        aGuess(2:7)=c(2:7)/sqrt(CC(2));
    else
        a0Guess=(c0-k1(1))/sqrt(CC(2));
        aGuess(1:7)=(c(1:7))/sqrt(CC(2));
    end
        
    %We calculate coefficients of Standardized target density
    [c0,c] = CalculateDensityFromStandardizedMoments(rmu,a0Guess,aGuess);
    
    %From standardized Z-seriesCoefficients of target density, we change to
    %actual non-standardized Z-series coefficients.
    c0=c0*sqrt(k1(2))+k1(1);
    c=c*sqrt(k1(2));
 
    %Actual Raw Moments after the Calculation of target density.
    [MM2] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
    %Calculate the Ratio of Moments between intermediate density and the
    %new target density
    MomentsIncreaseRatio=max(MM2(2:8)./MM1(2:8));
    
    %Below decrease the step size if MomentsIncreaseRatio is larger
    %than  MaxMomentRatio and increase the step size if 
    %MomentsIncreaseRatio is smaller than  MinMomentRatio

    IncreaseRatioFlag=0;
    DecreaseRatioFlag=0;
    if(MomentsIncreaseRatio>MaxMomentRatio)
        if (ds(tt)>MinStepSize)
            DecreaseRatioFlag=1;
        end
    elseif (MomentsIncreaseRatio<MinMomentRatio)   
        if (ds(tt)<MaxStepSize)
            IncreaseRatioFlag=1;
        end
    end
    
    if(DecreaseRatioFlag==1)
        ds(tt+1)=ds(tt)/2;
    
    elseif(IncreaseRatioFlag==1)
        ds(tt+1)=ds(tt)*2;
    else
        ds(tt+1)=ds(tt);
    end    
    
    %Tx is running time.
    Tx=Tx+ds(tt);
    
    %At last step make sure that step size is appropriate.
    if(Tx+ds(tt+1)>T)
        ds(tt+1)=T-Tx;
    end
    

    %Assign the median of density to w0
    w0=c0;
    
    end

end

ttFirst=tt;
wnStart=1;


%Below w which is the SDE in bessel coordinates and it is calulated from
%Z-series in Bessel coordinates.
w(1:Nn)=c0;
for nn=1:ExpnOrder
    w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^nn;
end
Flag=0;
for nn=Nn-1:-1:1
    if((w(nn)<0)&&(Flag==0))
        wnStart=nn+1;
        Flag=1;
    end
end


%Below we directly transform the density of SDE in bessel coordinates into
%density of SDE in original coordinates.

yy(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
%yy0=((1-gamma).*w0).^(1/(1-gamma));

%Below we make calculations of density in original coordinates by doing 
%change of probability derivative with respect to gaussian density.
Dfyy(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    
    Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    
    %Change of variable derivative for densities
end
Dfyy(Nn)=Dfyy(Nn-1);
Dfyy(1)=Dfyy(2);


pyy(1:Nn)=0;
for nn = wnStart:Nn
    
    pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
end


toc

%Below we do calculations for change of Z-series coefficients from
%coordinates in Bessel to Z-series coefficients in original coordinates of
%the SDE.
%The new series is expanded around c0 which is median of the density in
%Bessel coordinates.
   Bmu=c0;%+c(2)+3*c(4)+15*c(6)
   Bmu
   
%Calculate the first seven derivatives of transformation from Bessel
%coordinates to original coordinates of SDE at the median value Bmu=c0.
   
  yB0=((1-gamma).*Bmu).^(1/(1-gamma));
  dyB0(1)=((1-gamma).*Bmu).^(1/(1-gamma)-1);
  dyB0(2)=(1/(1-gamma)-1).*((1-gamma).*Bmu).^(1/(1-gamma)-2)*(1-gamma);
  dyB0(3)=(1/(1-gamma)-1).*(1/(1-gamma)-2).*((1-gamma).*Bmu).^(1/(1-gamma)-3)*(1-gamma)^2;
  dyB0(4)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*((1-gamma).*Bmu).^(1/(1-gamma)-4)*(1-gamma)^3;
  dyB0(5)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*((1-gamma).*Bmu).^(1/(1-gamma)-5)*(1-gamma)^4;
  dyB0(6)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*((1-gamma).*Bmu).^(1/(1-gamma)-6)*(1-gamma)^5;
  dyB0(7)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*(1/(1-gamma)-6)*((1-gamma).*Bmu).^(1/(1-gamma)-7)*(1-gamma)^6;
%  dyB0(8)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*(1/(1-gamma)-6)*(1/(1-gamma)-7)*((1-gamma).*Bmu).^(1/(1-gamma)-8)*(1-gamma)^7;
%  dyB0(9)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*(1/(1-gamma)-6)*(1/(1-gamma)-7)*(1/(1-gamma)-8)*((1-gamma).*Bmu).^(1/(1-gamma)-9)*(1-gamma)^8;
%  dyB0(10)=(1/(1-gamma)-1).*(1/(1-gamma)-2)*(1/(1-gamma)-3)*(1/(1-gamma)-4)*(1/(1-gamma)-5)*(1/(1-gamma)-6)*(1/(1-gamma)-7)*(1/(1-gamma)-8)*(1/(1-gamma)-9)*((1-gamma).*Bmu).^(1/(1-gamma)-10)*(1-gamma)^9;

c(8:10)=0;

%Below finally calculate the Z-series coefficients from Bessel coordinates 
%to original coordinates using derivatives calculated at previous step.
   [d0,d] = CalculateDriftbCoeffs08(yB0,dyB0,c,7);
 
%Below calculate SDE variable in original coordinates after expanding the 
%Z-series coefficients in original coordinates.   
Yw(1:Nn)=d0;
for nn=1:7
    Yw(1:Nn)=Yw(1:Nn)+d(nn)*Z(1:Nn).^nn;
end

%Below calculate the density from Z-series Coefficients in original
%coordinates. We want to see how exact they are since we would need to use
%them in analytics later when calculating integrals of variance for
%stochastic volatility models.

DfYw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    
    DfYw(nn) = (Yw(nn + 1) - Yw(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    
    %Change of variable derivative for densities
end
DfYw(Nn)=DfYw(Nn-1);
DfYw(1)=DfYw(2);


pYw(1:Nn)=0;
for nn = wnStart:Nn
    
    pYw(nn) = (normpdf(Z(nn),0, 1))/abs(DfYw(nn));
end

%AnalyticMean1 is calculated from direct transformation of density in Bessel coordinates
AnalyticMean1=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from transformed density
%AnalyticMean2 is calculated from Z-series coordinates of density in Original coordinates
AnalyticMean2=d0+d(2)+3*d(4)+15*d(6) %Original process average from coordinates 
disp('true Mean only applicable to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average
%yy0
%yB0

%theta1=1;

%Monte Carlo starts here. For details see my thread on wilmott.com 
%where I gave details of higher order monte carlo.

rng(29079137, 'twister')
paths=50000;
YY(1:paths)=x0;  %Original process monte carlo.
Random1(1:paths)=0;
for tt=1:TtM
    
    
    Random1=randn(size(Random1));
    HermiteP1(1,1:paths)=1;
    HermiteP1(2,1:paths)=Random1(1:paths);
    HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
    HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
    HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;

 
    YY(1:paths)=YY(1:paths) + ...
        (YCoeff0(1,1,2,1).*YY(1:paths).^Fp(1,1,2,1)+ ...
        YCoeff0(1,2,1,1).*YY(1:paths).^Fp(1,2,1,1)+ ...
        YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1))*dtM + ...
        (YCoeff0(1,1,3,1).*YY(1:paths).^Fp(1,1,3,1)+ ...
        YCoeff0(1,2,2,1).*YY(1:paths).^Fp(1,2,2,1)+ ...
        YCoeff0(2,1,2,1).*YY(1:paths).^Fp(2,1,2,1)+ ...
        YCoeff0(1,3,1,1).*YY(1:paths).^Fp(1,3,1,1)+ ...
        YCoeff0(2,2,1,1).*YY(1:paths).^Fp(2,2,1,1)+ ...
        YCoeff0(3,1,1,1).*YY(1:paths).^Fp(3,1,1,1))*dtM^2 + ...
        ((YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
        (YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+ ...
        YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
        YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5) .*HermiteP1(2,1:paths) + ...
        ((YCoeff0(1,1,1,3).*YY(1:paths).^Fp(1,1,1,3) *dtM) + ...
        (YCoeff0(1,1,2,3).*YY(1:paths).^Fp(1,1,2,3)+ ...
        YCoeff0(1,2,1,3).*YY(1:paths).^Fp(1,2,1,3)+ ...
        YCoeff0(2,1,1,3).*YY(1:paths).^Fp(2,1,1,3)).*dtM^2).*HermiteP1(3,1:paths) + ...
        ((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )).*HermiteP1(4,1:paths) + ...
        (YCoeff0(1,1,1,5).*YY(1:paths).^Fp(1,1,1,5)*dtM^2.0).*HermiteP1(5,1:paths);
    
    
end




YY(YY<0)=0;
disp('Original process average from monte carlo');
MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
disp('Original process average from our simulation');
ItoHermiteMean=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates 
 
disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average


MaxCutOff=100;
NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too 
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
plot(Yw(wnStart+1:Nn-1),pYw(wnStart+1:Nn-1),'b',yy(wnStart+1:Nn-1),pyy(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 %plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');

 
dtAvg=T/ttFirst; 
 
title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dtAvg=%.5f,M=%.4f,TM=%.4f', x0,theta,kappa,gamma,sigma0,T,dtAvg,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
legend({'SDE density From Z-series in Original Coordinates','SDE density Transformed Directly From Bessel Coordinates','Monte Carlo Density'},'Location','northeast')
 




str=input('red line is density of SDE from Addition of Cumulants, green is monte carlo.');


plot((wnStart+1:Nn-1)*dNn-NnMid,Yw(wnStart+1:Nn-1),'b',(wnStart+1:Nn-1)*dNn-NnMid,yy(wnStart+1:Nn-1),'r');

title('SDE Stochastic Variable as a Function of Standard Normal');%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
legend({'SDE Variable From Z-series in Original Coordinates','SDE Variable Transformed Directly From Bessel Coordinates'},'Location','northwest')

str=input('red line is density of SDE from addition of Cumulants method, green is monte carlo.');



end
.
.
.
function [Y] = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma)


%In the coefficient calculation program which calculates Y(l1,l2,l3,l4), 
%I have used four levels of looping each for relevant expansion order. 
%The first loop takes four values and second loop takes 16 values and 
%third loop takes 64 values and so on. And then each coefficient 
%term can be individually calculated while carefully accounting 
%for path dependence. 
%So for example in a nested loop structure 
%m1= 1:mDim
% m2=1:mDim
% m3=1:mDim
%    l(m1)=l(m1)+1;
%    l(m2)=l(m2)+1;
%    l(m3)=l(m3)+1;

%in the above looping loop takes values from one to four with one 
%indicating the first drift term, two indicating the second drift term 
%and three indicating quadratic variation term and 
%four indicating the volatility term. And with this looping structure 
%we can So in the above looping m1=1 would mean that all terms are 
%descendent of first drift term and m2=4 would mean that all terms are 
%descendent of first drift term on first expansion order and descendent 
%of volatility term on the second order and so we can keep track of path 
%dependence perfectly. 
%And on each level, we individually calculate the descendent terms. While 
%keeping track of path dependence and calculating the coefficients with 
%careful path dependence consideration, we update the appropriate element 
%in our polynomial like expansion coefficient array 

%explaining the part of code
%m1= 1:mDim
% m2=1:mDim
% m3=1:mDim
%    l(m1)=l(m1)+1;
%    l(m2)=l(m2)+1;
%    l(m3)=l(m3)+1;
%Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2);

%Here l(1) denotes l1 but written as l(1) so it can be conveniently 
%updated with the loop variable when the loop variable takes value one 
%indicating first drift term . And l(2) could be conveniently updated when 
%the loop variable takes value equal to two indicating second 
%drift term and so on.
%Here is the part of code snippet for that

%for m1=1:mDim
%    l(1)=1;
%    l(2)=1;
%    l(3)=1;
%    l(4)=1;
%    l(m1)=l(m1)+1;
%CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
%    CoeffDX2 = CoeffDX1 - 1;
%    ArrIndex0=m1;
%    ArrIndex=(m1-1)*mDim;
%    Coeff1st=Y1(ArrIndex0)*CoeffDX1;
%    Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
%    Y2(1+ArrIndex)=Coeff1st;
%    Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
%    Y2(2+ArrIndex)=Coeff1st;
%    Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
%    Y2(3+ArrIndex)=Coeff2nd;
%    Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,2,n1(m1));
%    Y2(4+ArrIndex)=Coeff1st;
%    Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,3,n1(m1));

 
%The first four lines update the array indices according to the parent term. 
%And then CoeffDX1 and CoeffDX2 are calculated according to algebraic exponents on parent terms.    
%ArrIndex0=m1; calculates the array index of the parent term 
%And ArrIndex=(m1-1)*mDim; calculates the array index of the descendent terms
%And coefficient of the drift and volatility descendent terms is 
%calculated by multiplying the coefficient of the parent term by 
%Coeff1st=Y1(ArrIndex0)*CoeffDX1;
%And coefficient of the quadratic variation descendent terms is 
%calculated by multiplying the coefficient of the parent term by 
%Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
%And then each of the four descendent terms are updated with Coeff1st 
%if they are drift or volatility descendent terms or Coeff2nd if 
%they are quadratic variation descendent terms.
%Here Y1 indicates the temporary coefficient array with parent terms on 
%first level. And Y2 denotes temporary coefficient array with parent terms 
%on second level and Y3 indicates temporary coefficient array with parent terms on third level.

[IntegralCoeff,IntegralCoeffdt,IntegralCoeffdz] = ComputeIntegralCoeffs();

n1(1)=2;
n1(2)=2;
n1(3)=2;
n1(4)=3;
%n1(1), n1(2), n1(3) are drift and quadratic variation term variables 
%and take a value equal to 2 indicating a dt integral.
%n1(4) is volatility term variable and indicates a dz-integral by taking a
%value of 3. 

mDim=4; % four descendent terms in each expansion
Y(1:5,1:5,1:5,1:5)=0;
Y1(1:mDim)=0;
Y2(1:mDim*mDim)=0;
Y3(1:mDim*mDim*mDim)=0;


%First Ito-hermite expansion level starts here. No loops but four
%descendent terms.
l(1)=1;
l(2)=1;
l(3)=1;
l(4)=1;
CoeffDX1 = alpha;
CoeffDX2 = CoeffDX1 - 1;
Coeff1st=CoeffDX1;
Coeff2nd=.5*CoeffDX1*CoeffDX2;

Y1(1)=Coeff1st;
Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2);
Y1(2)=Coeff1st;
Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,1,2);
Y1(3)=Coeff2nd;
Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,1,2);
Y1(4)=Coeff1st;
Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,1,3);

%Second Ito-hermite expansion level starts. It has a loop over four parent
%terms and there are four descendent terms for each parent term.
%The coefficient terms are then lumped in a polynomial-like expansion
%array of coefficients.
for m1=1:mDim
    l(1)=1;
    l(2)=1;
    l(3)=1;
    l(4)=1;
    l(m1)=l(m1)+1;
    CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
    CoeffDX2 = CoeffDX1 - 1;
    ArrIndex0=m1;
    ArrIndex=(m1-1)*mDim;
    Coeff1st=Y1(ArrIndex0)*CoeffDX1;
    Coeff2nd=Y1(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
    Y2(1+ArrIndex)=Coeff1st;
    Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
    Y2(2+ArrIndex)=Coeff1st;
    Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,1,2,n1(m1));
    Y2(3+ArrIndex)=Coeff2nd;
    Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,1,2,n1(m1));
    Y2(4+ArrIndex)=Coeff1st;
    Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,1,3,n1(m1));
    %Third Ito-hermite expansion level starts and it is a nested loop with
    %a total of sixteen parents and each parent takes four descendent
    %terms.
    for m2=1:mDim
        l(1)=1;
        l(2)=1;
        l(3)=1;
        l(4)=1;
        l(m1)=l(m1)+1;
        l(m2)=l(m2)+1;
        CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
        CoeffDX2=CoeffDX1-1;            
        ArrIndex0=(m1-1)*mDim+m2;
        ArrIndex=((m1-1)*mDim+(m2-1))*mDim;
        Coeff1st=Y2(ArrIndex0)*CoeffDX1;
        Coeff2nd=Y2(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
        Y3(1+ArrIndex)=Coeff1st;
        Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(1,2,n1(m2),n1(m1));
        Y3(2+ArrIndex)=Coeff1st;
        Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(1,2,n1(m2),n1(m1));
        Y3(3+ArrIndex)=Coeff2nd;
        Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(1,2,n1(m2),n1(m1));
        Y3(4+ArrIndex)=Coeff1st;
        Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(1,3,n1(m2),n1(m1));
        %fourht Ito-hermite expansion level starts and it is a triply-nested loop with
        %a total of sixteen parents and each parent takes four descendent
        %terms. We then lump the terms in a relatively sparse polynomial 
        %like expansion coefficient array that has smaller number of
        %non-zero terms.

        for m3=1:mDim
            l(1)=1;
            l(2)=1;
            l(3)=1;
            l(4)=1;
            l(m1)=l(m1)+1;
            l(m2)=l(m2)+1;
            l(m3)=l(m3)+1;
            CoeffDX1 = alpha + (l(1)-1) *beta1 + (l(2)-1) *beta2 + (l(3)-1) *2*gamma + (l(4)-1)*gamma ...
                        - (l(1)-1) - (l(2)-1) - 2*(l(3)-1) - (l(4)-1);
            CoeffDX2=CoeffDX1-1;
            ArrIndex0=((m1-1)*mDim+(m2-1))*mDim+m3;
            %ArrIndex=(((m1-1)*mDim+(m2-1))*mDim+(m3-1))*mDim;
            Coeff1st=Y3(ArrIndex0)*CoeffDX1;
            Coeff2nd=Y3(ArrIndex0)*.5*CoeffDX1*CoeffDX2;
            %Y4(1+ArrIndex)=Coeff1st;
            Y(l(1)+1,l(2),l(3),l(4))=Y(l(1)+1,l(2),l(3),l(4))+Coeff1st*IntegralCoeff(2,n1(m3),n1(m2),n1(m1));
            %Y4(2+ArrIndex)=Coeff1st;
            Y(l(1),l(2)+1,l(3),l(4))=Y(l(1),l(2)+1,l(3),l(4))+Coeff1st*IntegralCoeff(2,n1(m3),n1(m2),n1(m1));
            %Y4(3+ArrIndex)=Coeff2nd;
            Y(l(1),l(2),l(3)+1,l(4))=Y(l(1),l(2),l(3)+1,l(4))+Coeff2nd*IntegralCoeff(2,n1(m3),n1(m2),n1(m1));
            %Y4(4+ArrIndex)=Coeff1st;
            Y(l(1),l(2),l(3),l(4)+1)=Y(l(1),l(2),l(3),l(4)+1)+Coeff1st*IntegralCoeff(3,n1(m3),n1(m2),n1(m1));
         end
    end
end
        
        
        

end

.
.
.
function [IntegralCoeff0,IntegralCoeffdt,IntegralCoeffdz] = ComputeIntegralCoeffs()


IntegralCoeff(1:3,1:3,1:3,1:3)=0;
IntegralCoeffdt(1:3,1:3,1:3,1:3)=0;
IntegralCoeffdz(1:3,1:3,1:3,1:3)=0;

IntegralCoeff0(1:3,1:3,1:3,1:3)=0;

IntegralCoeff0(1,1,1,2)=1;
IntegralCoeff0(1,1,1,3)=1;

IntegralCoeff0(1,1,2,2)=1/2;
IntegralCoeff0(1,1,3,2)=1-1/sqrt(3);
IntegralCoeff0(1,1,2,3)=1/sqrt(3);
IntegralCoeff0(1,1,3,3)=1/2;

IntegralCoeff0(1,2,2,2)=1/6;
IntegralCoeff0(1,3,2,2)=(1-1/sqrt(3))*1/2*(1-1/sqrt(5));
IntegralCoeff0(1,2,3,2)=1/sqrt(3)/2*(1-1/sqrt(5));
IntegralCoeff0(1,3,3,2)=1/2*(1-sqrt(2)/2);
IntegralCoeff0(1,2,2,3)=1/2*1/sqrt(5);
IntegralCoeff0(1,3,2,3)=(1-1/sqrt(3))*1/sqrt(2)*1/2;
IntegralCoeff0(1,2,3,3)=1/sqrt(3)/sqrt(2)*1/(2);
IntegralCoeff0(1,3,3,3)=1/6;

IntegralCoeff0(2,2,2,2)=1/24;
IntegralCoeff0(2,2,3,2)=1/2*1/sqrt(5)*1/3*(1-1/sqrt(7));
IntegralCoeff0(2,3,2,2)=1/sqrt(3)*1/2*(1-1/sqrt(5))*1/3*(1-1/sqrt(7));
IntegralCoeff0(3,2,2,2)=(1-1/sqrt(3))*1/2*(1-1/sqrt(5))*1/3*(1-1/sqrt(7));
IntegralCoeff0(3,3,2,2)=1/2*(1-sqrt(2)/2)*1/2*(1-sqrt(2)/sqrt(6));
IntegralCoeff0(2,3,3,2)=1/sqrt(3)*1/sqrt(2)*1/2*1/2*(1-sqrt(2)/sqrt(6));
IntegralCoeff0(3,2,3,2)=(1-1/sqrt(3))*1/sqrt(2)*1/2*1/2*(1-sqrt(2)/sqrt(6));
IntegralCoeff0(3,3,3,2)=1/6*(1-sqrt(3)/sqrt(5));

IntegralCoeff0(2,2,2,3)=1/6*1/sqrt(7);
IntegralCoeff0(2,2,3,3)=1/2*1/sqrt(5)*1/sqrt(2)*1/sqrt(6);
IntegralCoeff0(2,3,2,3)=1/sqrt(3)*1/2*(1-1/sqrt(5))*1/sqrt(2)*1/sqrt(6);
IntegralCoeff0(3,2,2,3)=(1-1/sqrt(3))*1/2*(1-1/sqrt(5))*1/sqrt(2)*1/sqrt(6);
IntegralCoeff0(3,3,2,3)=1/2*(1-sqrt(2)/2)*1/sqrt(3)*1/sqrt(5);
IntegralCoeff0(2,3,3,3)=1/sqrt(3)*1/sqrt(2)*1/sqrt(4)*1/sqrt(3)*(1/sqrt(5));
IntegralCoeff0(3,2,3,3)=(1-1/sqrt(3))*1/sqrt(2)*1/2*1/sqrt(3)*1/sqrt(5);
IntegralCoeff0(3,3,3,3)=1/24;


%Can also be calculated with algorithm below.
%here IntegralCoeffdt indicates the coefficients of a dt-integral.
%This dt-integral means that you are calculating the Ito-hermite expansion
%of  Integral X(t)^alpha dt with X(t) dynamics given by the SDE 
%here IntegralCoeffdz indicates the coefficients of a dz-integral.
%This dz-integral means that you are calculating the Ito-hermite expansion
%of  Integral X(t)^alpha dz(t) with X(t) dynamics given by the SDE 

%IntegralCoeff is is associated with expansion of X(t)^alpha.
%IntegralCoeff below is not returned and IntegralCoeff0 manually calculated
%above is returned but both are the same.

l0(1:2)=1;

for m4=1:2
    l0(1)=1;
    l0(2)=1;
    
    %IntegralCoeff4(m4,1,1,1)=1;
    %IntegralCoeff4(m4,1,1,1)=1;
    %1 is added to m4 since index 1 stands for zero, 2 for one and three
    %for two.
    IntegralCoeff(1,1,1,m4+1)=1;
    l0(m4)=l0(m4)+1;
    IntegralCoeffdt(1,1,1,m4+1)=IntegralCoeff(1,1,1,m4+1)* ...
        1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
    IntegralCoeffdz(1,1,1,m4+1)= IntegralCoeff(1,1,1,m4+1)* ...
        1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
    for m3=1:2
        l0(1)=1;
        l0(2)=1;
        l0(m4)=l0(m4)+1;
        
        if(m3==1)
            IntegralCoeff(1,1,m4+1,m3+1)=1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
        end
        if(m3==2)
            IntegralCoeff(1,1,m4+1,m3+1)= 1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
        end
        l0(m3)=l0(m3)+1;
        %IntegralCoeff(1,1,m4+1,m3+1)=IntegralCoeff4(m4,m3,1,1);
        IntegralCoeffdt(1,1,m4+1,m3+1)=IntegralCoeff(1,1,m4+1,m3+1)* ...
            1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
        IntegralCoeffdz(1,1,m4+1,m3+1)= IntegralCoeff(1,1,m4+1,m3+1)* ...
            1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
        for m2=1:2
            l0(1)=1;
            l0(2)=1;
            l0(m4)=l0(m4)+1;
            l0(m3)=l0(m3)+1;
            
            if(m2==1)
                IntegralCoeff(1,m4+1,m3+1,m2+1)=IntegralCoeff(1,1,m4+1,m3+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
            end
            if(m2==2)
                IntegralCoeff(1,m4+1,m3+1,m2+1)= IntegralCoeff(1,1,m4+1,m3+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
            end
            l0(m2)=l0(m2)+1;
            %IntegralCoeff(1,m4+1,m3+1,m2+1)=IntegralCoeff4(m4,m3,m2,1);
            IntegralCoeffdt(1,m4+1,m3+1,m2+1)=IntegralCoeff(1,m4+1,m3+1,m2+1)* ...
                1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
            IntegralCoeffdz(1,m4+1,m3+1,m2+1)= IntegralCoeff(1,m4+1,m3+1,m2+1)* ...
                1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);

            for m1=1:2
                l0(1)=1;
                l0(2)=1;
                l0(m4)=l0(m4)+1;
                l0(m3)=l0(m3)+1;
                l0(m2)=l0(m2)+1;
                if(m1==1)
                    IntegralCoeff(m4+1,m3+1,m2+1,m1+1)=IntegralCoeff(1,m4+1,m3+1,m2+1)*1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                end
                if(m1==2)
                    IntegralCoeff(m4+1,m3+1,m2+1,m1+1)= IntegralCoeff(1,m4+1,m3+1,m2+1)*1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
                end
                l0(m1)=l0(m1)+1;
                %IntegralCoeff(m4+1,m3+1,m2+1,m1+1)=IntegralCoeff4(m4,m3,m2,m1);
                IntegralCoeffdt(m4+1,m3+1,m2+1,m1+1)=IntegralCoeff(m4+1,m3+1,m2+1,m1+1)* ...
                    1/(l0(1)-1+1)*(1-sqrt(l0(2)-1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+2));
                IntegralCoeffdz(m4+1,m3+1,m2+1,m1+1)= IntegralCoeff(m4+1,m3+1,m2+1,m1+1)* ...
                    1/sqrt(l0(2)-1+1)/sqrt(2*(l0(1)-1)+(l0(2)-1)+1);
            end
        end
    end
end

end
.
.
.
function [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder)


NoOfTerms=19;
%NoOfTerms=9;

YqCoeffa(1:NoOfTerms)=0.0;
%Fp1
Fp1=Fp1/(1-gamma);
%Fp1
%gamma
YqCoeffa(1)=YqCoeff0(1,1,2,1)*dt;%*(1-gamma)^Fp1(1,1,2,1)*dt;
YqCoeffa(2)=YqCoeff0(1,2,1,1)*dt;%*(1-gamma)^Fp1(1,2,1,1)*dt;
YqCoeffa(3)=YqCoeff0(2,1,1,1)*dt;%*(1-gamma)^Fp1(2,1,1,1)*dt;
YqCoeffa(4)=YqCoeff0(1,1,3,1)*dt^2;%*(1-gamma)^Fp1(1,1,3,1)*dt^2;
YqCoeffa(5)=YqCoeff0(1,2,2,1)*dt^2;%*(1-gamma)^Fp1(1,2,2,1)*dt^2;
YqCoeffa(6)=YqCoeff0(2,1,2,1)*dt^2;%*(1-gamma)^Fp1(2,1,2,1)*dt^2;
YqCoeffa(7)=YqCoeff0(1,3,1,1)*dt^2;%*(1-gamma)^Fp1(1,3,1,1)*dt^2;
YqCoeffa(8)=YqCoeff0(2,2,1,1)*dt^2;%*(1-gamma)^Fp1(2,2,1,1)*dt^2;
YqCoeffa(9)=YqCoeff0(3,1,1,1)*dt^2;%*(1-gamma)^Fp1(3,1,1,1)*dt^2;
YqCoeffa(10)=YqCoeff0(1,1,4,1)*dt^3;%;*(1-gamma)^Fp1(1,1,4,1)*dt^3;
YqCoeffa(11)=YqCoeff0(1,2,3,1)*dt^3;%;*(1-gamma)^Fp1(1,2,3,1)*dt^3;
YqCoeffa(12)=YqCoeff0(2,1,3,1)*dt^3;%*(1-gamma)^Fp1(2,1,3,1)*dt^3;
YqCoeffa(13)=YqCoeff0(1,3,2,1)*dt^3;%*(1-gamma)^Fp1(1,3,2,1)*dt^3;
YqCoeffa(14)=YqCoeff0(2,2,2,1)*dt^3;%*(1-gamma)^Fp1(2,2,2,1)*dt^3;
YqCoeffa(15)=YqCoeff0(3,1,2,1)*dt^3;%*(1-gamma)^Fp1(3,1,2,1)*dt^3;
YqCoeffa(16)=YqCoeff0(1,4,1,1)*dt^3;%*(1-gamma)^Fp1(1,4,1,1)*dt^3;
YqCoeffa(17)=YqCoeff0(2,3,1,1)*dt^3;%*(1-gamma)^Fp1(2,3,1,1)*dt^3;
YqCoeffa(18)=YqCoeff0(3,2,1,1)*dt^3;%*(1-gamma)^Fp1(3,2,1,1)*dt^3;
YqCoeffa(19)=YqCoeff0(4,1,1,1)*dt^3;%*(1-gamma)^Fp1(4,1,1,1)*dt^3;

Fp2(1)=Fp1(1,1,2,1);
Fp2(2)=Fp1(1,2,1,1);
Fp2(3)=Fp1(2,1,1,1);
Fp2(4)=Fp1(1,1,3,1);
Fp2(5)=Fp1(1,2,2,1);
Fp2(6)=Fp1(2,1,2,1);
Fp2(7)=Fp1(1,3,1,1);
Fp2(8)=Fp1(2,2,1,1);
Fp2(9)=Fp1(3,1,1,1);
Fp2(10)=Fp1(1,1,4,1);
Fp2(11)=Fp1(1,2,3,1);
Fp2(12)=Fp1(2,1,3,1);
Fp2(13)=Fp1(1,3,2,1);
Fp2(14)=Fp1(2,2,2,1);
Fp2(15)=Fp1(3,1,2,1);
Fp2(16)=Fp1(1,4,1,1);
Fp2(17)=Fp1(2,3,1,1);
Fp2(18)=Fp1(3,2,1,1);
Fp2(19)=Fp1(4,1,1,1);

%YqCoeffa
%Fp2
%str=input('Look at numbers');
wMu0dt0=0;
dwMu0dt(1:ExpnOrder)=0.0;

wMu0dt=0;
dwMu0dtdw(1:ExpnOrder)=0.0;

for mm=1:NoOfTerms

    wMu0dt0=YqCoeffa(mm).*((1-gamma)*w0).^Fp2(mm);
    for nn=1:ExpnOrder
        if(nn==1)
            dwMu0dt(nn)=wMu0dt0*Fp2(mm)/w0;
        else
        dwMu0dt(nn)=dwMu0dt(nn-1)*(Fp2(mm)-(nn-1))/w0;
        end
    end
    wMu0dt=wMu0dt+wMu0dt0;
    for nn=1:ExpnOrder
        dwMu0dtdw(nn)=dwMu0dtdw(nn)+dwMu0dt(nn);
    end
        
    
end



end

.
.
.
function [b0,b] = CalculateDriftbCoeffs08(wMu0dt,dwMu0dtdw,a,SeriesOrder)



b0=wMu0dt;

b(1)=a(1) *dwMu0dtdw(1);

b(2)=1/2*(2* a(2) *dwMu0dtdw(1)+a(1)^2 *dwMu0dtdw(2));

b(3)=1/6*(6 *a(3) *dwMu0dtdw(1)+6 *a(1)* a(2) *dwMu0dtdw(2)+a(1)^3 *dwMu0dtdw(3));

b(4)=1/24*(24* a(4) *dwMu0dtdw(1)+12* a(2)^2 *dwMu0dtdw(2)+24 *a(1)* a(3) *dwMu0dtdw(2)+ ...
    12* a(1)^2* a(2) *dwMu0dtdw(3)+a(1)^4 *dwMu0dtdw(4));

if(SeriesOrder>=5)
b(5)=1/120*(120* a(5) *dwMu0dtdw(1)+120* a(2) *a(3) *dwMu0dtdw(2)+120* a(1)* a(4) *dwMu0dtdw(2)+ ...
    60 *a(1)* a(2)^2 *dwMu0dtdw(3)+60* a(1)^2* a(3) *dwMu0dtdw(3)+ ...
    20 *a(1)^3* a(2) *dwMu0dtdw(4)+a(1)^5 *dwMu0dtdw(5));
end
if(SeriesOrder>=6)
b(6)=1/720*(720* a(6) *dwMu0dtdw(1)+360* a(3)^2 *dwMu0dtdw(2)+720* a(2)* a(4) *dwMu0dtdw(2)+720* a(1)* a(5) *dwMu0dtdw(2)+ ...
    120* a(2)^3 *dwMu0dtdw(3)+720* a(1)* a(2)* a(3) *dwMu0dtdw(3)+360* a(1)^2* a(4) *dwMu0dtdw(3)+ ...
    180* a(1)^2* a(2)^2 *dwMu0dtdw(4)+120* a(1)^3* a(3) *dwMu0dtdw(4)+ ...
    30* a(1)^4 *a(2) *dwMu0dtdw(5)+a(1)^6 *dwMu0dtdw(6));
end


if(SeriesOrder>=7)
  b(7)=1/(720*7)*(5040* a(7)*dwMu0dtdw(1)+5040* a(3)* a(4)*dwMu0dtdw(2)+5040* a(2)* a(5) *dwMu0dtdw(2)+5040 *a(1)* a(6) *dwMu0dtdw(2)+ ...
      2520* a(2)^2* a(3) *dwMu0dtdw(3)+2520* a(1)* a(3)^2 *dwMu0dtdw(3)+5040 *a(1)* a(2)* a(4) *dwMu0dtdw(3)+2520* a(1)^2 *a(5) *dwMu0dtdw(3)+ ...
      840* a(1)* a(2)^3 *dwMu0dtdw(4)+2520* a(1)^2* a(2)* a(3) *dwMu0dtdw(4)+840* a(1)^3* a(4) *dwMu0dtdw(4)+ ...
      420* a(1)^3* a(2)^2 *dwMu0dtdw(5)+210* a(1)^4* a(3) *dwMu0dtdw(5)+ ...
      42* a(1)^5* a(2) *dwMu0dtdw(6)+a(1)^7 *dwMu0dtdw(7));
  
end
  
if(SeriesOrder>=8)
  b(8)=1/(720*7*8)*(40320* a(8) *dwMu0dtdw(1)+20160* a(4)^2 *dwMu0dtdw(2)+40320* a(3)* a(5) *dwMu0dtdw(2)+40320* a(2) *a(6) *dwMu0dtdw(2)+ ...
      40320* a(1) *a(7) *dwMu0dtdw(2)+ ...
      20160* a(2)* a(3)^2 *dwMu0dtdw(3)+20160* a(2)^2* a(4) *dwMu0dtdw(3)+40320* a(1)* a(3)* a(4) *dwMu0dtdw(3)+40320* a(1)* a(2)* a(5) *dwMu0dtdw(3)+ ...
      20160* a(1)^2* a(6) *dwMu0dtdw(3)+ ...
      1680* a(2)^4 *dwMu0dtdw(4)+20160* a(1)* a(2)^2* a(3) *dwMu0dtdw(4)+10080* a(1)^2* a(3)^2 *dwMu0dtdw(4)+20160* a(1)^2* a(2) *a(4) *dwMu0dtdw(4)+ ...
      6720* a(1)^3* a(5) *dwMu0dtdw(4)+ ...
      3360* a(1)^2* a(2)^3 *dwMu0dtdw(5)+6720* a(1)^3 *a(2)* a(3) *dwMu0dtdw(5)+1680* a(1)^4* a(4) *dwMu0dtdw(5)+ ...
      840* a(1)^4* a(2)^2 *dwMu0dtdw(6)+336* a(1)^5* a(3) *dwMu0dtdw(6)+ ...
      56* a(1)^6* a(2) *dwMu0dtdw(7)+a(1)^8 *dwMu0dtdw(8));
end

if(SeriesOrder>=9)
b(9)=1/(720*7*8*9)*(362880* a(9) *dwMu0dtdw(1)+ ...
    362880* a(4)* a(5) *dwMu0dtdw(2)+ 362880* a(3)* a(6) *dwMu0dtdw(2) +362880* a(2)* a(7) *dwMu0dtdw(2) +362880* a(1) *a(8) *dwMu0dtdw(2) + ...
    60480* a(3)^3 *dwMu0dtdw(3)+362880* a(2)* a(3)* a(4) *dwMu0dtdw(3)+181440* a(1)* a(4)^2 *dwMu0dtdw(3)+181440* a(2)^2* a(5) *dwMu0dtdw(3) + ...
    362880* a(1)* a(3)* a(5) *dwMu0dtdw(3)+362880 *a(1)* a(2)* a(6) *dwMu0dtdw(3) +181440 *a(1)^2* a(7) *dwMu0dtdw(3)+ ...
    60480* a(2)^3* a(3) *dwMu0dtdw(4)+181440* a(1)* a(2)* a(3)^2 *dwMu0dtdw(4)+181440* a(1)* a(2)^2 *a(4) *dwMu0dtdw(4)+ ...
    181440* a(1)^2* a(3)* a(4) *dwMu0dtdw(4)+181440* a(1)^2* a(2)* a(5) *dwMu0dtdw(4)+60480* a(1)^3* a(6) *dwMu0dtdw(4)+ ...
    15120* a(1)* a(2)^4 *dwMu0dtdw(5)+90720* a(1)^2* a(2)^2* a(3) *dwMu0dtdw(5)+30240* a(1)^3* a(3)^2 *dwMu0dtdw(5)+ ...
    60480 *a(1)^3* a(2)* a(4) *dwMu0dtdw(5)+15120* a(1)^4* a(5) *dwMu0dtdw(5)+ ...
    10080* a(1)^3* a(2)^3 *dwMu0dtdw(6)+15120* a(1)^4 *a(2)* a(3) *dwMu0dtdw(6)+3024* a(1)^5* a(4) *dwMu0dtdw(6)+ ...
    1512* a(1)^5* a(2)^2 *dwMu0dtdw(7)+504* a(1)^6* a(3) *dwMu0dtdw(7)+ ...
    72* a(1)^7 *a(2) *dwMu0dtdw(8)+a(1)^9 *dwMu0dtdw(9));
end


if(SeriesOrder>=10)
b(10)=1/(720*7*8*9*10)*(3628800* a(10)* dwMu0dtdw(1)+ ...
    1814400 * a(5)^2 *dwMu0dtdw(2)+3628800 *a(4)* a(6)* dwMu0dtdw(2)+3628800* a(3)* a(7)* dwMu0dtdw(2)+3628800 *a(2)* a(8)* dwMu0dtdw(2)+ ...
    3628800* a(1)* a(9)* dwMu0dtdw(2)+ ...
    1814400* a(3)^2* a(4)* dwMu0dtdw(3)+1814400* a(2)* a(4)^2* dwMu0dtdw(3)+3628800* a(2)* a(3)* a(5)* dwMu0dtdw(3)+3628800* a(1)* a(4)* a(5) *dwMu0dtdw(3)+ ...
    1814400* a(2)^2* a(6)* dwMu0dtdw(3)+3628800* a(1)* a(3)* a(6)* dwMu0dtdw(3)+3628800 *a(1)* a(2)* a(7)* dwMu0dtdw(3)+1814400* a(1)^2* a(8)* dwMu0dtdw(3)+ ...
    907200* a(2)^2 *a(3)^2* dwMu0dtdw(4)+604800* a(1)* a(3)^3* dwMu0dtdw(4)+604800* a(2)^3* a(4)* dwMu0dtdw(4)+3628800* a(1)* a(2)* a(3)* a(4)* dwMu0dtdw(4)+ ...
    907200* a(1)^2* a(4)^2* dwMu0dtdw(4)+1814400* a(1)* a(2)^2* a(5)* dwMu0dtdw(4)+1814400* a(1)^2* a(3)* a(5)* dwMu0dtdw(4)+ ...
    1814400* a(1)^2* a(2)* a(6)* dwMu0dtdw(4)+ 604800* a(1)^3* a(7) *dwMu0dtdw(4)+ ...
    30240 *a(2)^5* dwMu0dtdw(5)+604800* a(1)* a(2)^3* a(3)* dwMu0dtdw(5)+907200 *a(1)^2* a(2)* a(3)^2* dwMu0dtdw(5)+907200* a(1)^2* a(2)^2* a(4) *dwMu0dtdw(5)+ ...
    604800* a(1)^3* a(3)* a(4)* dwMu0dtdw(5)+604800* a(1)^3* a(2)* a(5)* dwMu0dtdw(5)+151200* a(1)^4* a(6)* dwMu0dtdw(5)+ ...
    75600* a(1)^2* a(2)^4* dwMu0dtdw(6)+302400 *a(1)^3* a(2)^2* a(3)* dwMu0dtdw(6)+75600* a(1)^4 *a(3)^2* dwMu0dtdw(6)+151200* a(1)^4* a(2)* a(4) *dwMu0dtdw(6)+ ...
    30240* a(1)^5* a(5)* dwMu0dtdw(6)+ ...
    25200* a(1)^4* a(2)^3* dwMu0dtdw(7)+30240* a(1)^5* a(2)* a(3)* dwMu0dtdw(7)+5040* a(1)^6* a(4)* dwMu0dtdw(7)+ ...
    2520* a(1)^6* a(2)^2* dwMu0dtdw(8)+720* a(1)^7* a(3)* dwMu0dtdw(8)+ ...
    90 *a(1)^8* a(2)* dwMu0dtdw(9)+ ...
    a(1)^10* dwMu0dtdw(10));

end

end

.
.
.
function [wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder)




% c1(wnStart:Nn)=YqCoeff0(1,1,1,2).*yy(wnStart:Nn).^Fp1(1,1,1,2).*sqrt(dt)+ ...
%     (YqCoeff0(1,1,2,2).*yy(wnStart:Nn).^Fp1(1,1,2,2)+YqCoeff0(1,2,1,2).*yy(wnStart:Nn).^Fp1(1,2,1,2)+ ...
%     YqCoeff0(2,1,1,2).*yy(wnStart:Nn).^Fp1(2,1,1,2)).*dt^1.5+ ...
%     (YqCoeff0(1,1,3,2).*yy(wnStart:Nn).^Fp1(1,1,3,2)+YqCoeff0(1,2,2,2).*yy(wnStart:Nn).^Fp1(1,2,2,2)+ ...
%     YqCoeff0(2,1,2,2).*yy(wnStart:Nn).^Fp1(2,1,2,2)+YqCoeff0(1,3,1,2).*yy(wnStart:Nn).^Fp1(1,3,1,2)+ ...
%     YqCoeff0(2,2,1,2).*yy(wnStart:Nn).^Fp1(2,2,1,2)+YqCoeff0(3,1,1,2).*yy(wnStart:Nn).^Fp1(3,1,1,2)).*dt^2.5;

NoOfTerms=9;%excluding dt^3 terms
YqCoeffa(1:NoOfTerms)=0.0;

Fp1=Fp1/(1-gamma);
YqCoeffa(1)=YqCoeff0(1,1,2,2)*(1-gamma)^Fp1(1,1,2,2)*dt^1.5;
YqCoeffa(2)=YqCoeff0(1,2,1,2)*(1-gamma)^Fp1(1,2,1,2)*dt^1.5;
YqCoeffa(3)=YqCoeff0(2,1,1,2)*(1-gamma)^Fp1(2,1,1,2)*dt^1.5;
YqCoeffa(4)=YqCoeff0(1,1,3,2)*(1-gamma)^Fp1(1,1,3,2)*dt^2.5;
YqCoeffa(5)=YqCoeff0(1,2,2,2)*(1-gamma)^Fp1(1,2,2,2)*dt^2.5;
YqCoeffa(6)=YqCoeff0(2,1,2,2)*(1-gamma)^Fp1(2,1,2,2)*dt^2.5;
YqCoeffa(7)=YqCoeff0(1,3,1,2)*(1-gamma)^Fp1(1,3,1,2)*dt^2.5;
YqCoeffa(8)=YqCoeff0(2,2,1,2)*(1-gamma)^Fp1(2,2,1,2)*dt^2.5;
YqCoeffa(9)=YqCoeff0(3,1,1,2)*(1-gamma)^Fp1(3,1,1,2)*dt^2.5;

Fp2(1)=Fp1(1,1,2,2);
Fp2(2)=Fp1(1,2,1,2);
Fp2(3)=Fp1(2,1,1,2);
Fp2(4)=Fp1(1,1,3,2);
Fp2(5)=Fp1(1,2,2,2);
Fp2(6)=Fp1(2,1,2,2);
Fp2(7)=Fp1(1,3,1,2);
Fp2(8)=Fp1(2,2,1,2);
Fp2(9)=Fp1(3,1,1,2);


wVol0dt0=0;
dwVol0dt(1:ExpnOrder)=0.0;

wVol0dt=0;
dwVol0dtdw(1:ExpnOrder)=0.0;

for mm=1:NoOfTerms

    wVol0dt0=YqCoeffa(mm).*w0.^Fp2(mm);
    for nn=1:ExpnOrder
        if(nn==1)
            dwVol0dt(nn)=wVol0dt0*Fp2(mm)*1/w0;
        else
        dwVol0dt(nn)=dwVol0dt(nn-1)*(Fp2(mm)-(nn-1))/w0;
        end
    end
    wVol0dt=wVol0dt+wVol0dt0;
    for nn=1:ExpnOrder
        dwVol0dtdw(nn)=dwVol0dtdw(nn)+dwVol0dt(nn);
    end
        
    
end



end
.
.
.
function [wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder)


% c2(wnStart:Nn)=YqCoeff0(1,1,1,3).*yy(wnStart:Nn).^Fp1(1,1,1,3) *dt + ...
%     (YqCoeff0(1,1,2,3).*yy(wnStart:Nn).^Fp1(1,1,2,3)+YqCoeff0(1,2,1,3).*yy(wnStart:Nn).^Fp1(1,2,1,3)+ ...
%     YqCoeff0(2,1,1,3).*yy(wnStart:Nn).^Fp1(2,1,1,3)).*dt^2+ ...
%     (YqCoeff0(1,1,3,3).*yy(wnStart:Nn).^Fp1(1,1,3,3)+YqCoeff0(1,2,2,3).*yy(wnStart:Nn).^Fp1(1,2,2,3)+ ...
%     YqCoeff0(2,1,2,3).*yy(wnStart:Nn).^Fp1(2,1,2,3) + YqCoeff0(1,3,1,3).*yy(wnStart:Nn).^Fp1(1,3,1,3)+ ...
%     YqCoeff0(2,2,1,3).*yy(wnStart:Nn).^Fp1(2,2,1,3)+YqCoeff0(3,1,1,3).*yy(wnStart:Nn).^Fp1(3,1,1,3)).*dt^3;

Fp1=Fp1/(1-gamma);
 
NoOfTerms=9;%excluding dt^3 terms
YqCoeffa(1:NoOfTerms)=0.0;
YqCoeffa(1)=YqCoeff0(1,1,2,3)*(1-gamma)^Fp1(1,1,2,3)*dt^2;
YqCoeffa(2)=YqCoeff0(1,2,1,3)*(1-gamma)^Fp1(1,2,1,3)*dt^2;
YqCoeffa(3)=YqCoeff0(2,1,1,3)*(1-gamma)^Fp1(2,1,1,3)*dt^2;
YqCoeffa(4)=YqCoeff0(1,1,3,3)*(1-gamma)^Fp1(1,1,3,3)*dt^3;
YqCoeffa(5)=YqCoeff0(1,2,2,3)*(1-gamma)^Fp1(1,2,2,3)*dt^3;
YqCoeffa(6)=YqCoeff0(2,1,2,3)*(1-gamma)^Fp1(2,1,2,3)*dt^3;
YqCoeffa(7)=YqCoeff0(1,3,1,3)*(1-gamma)^Fp1(1,3,1,3)*dt^3;
YqCoeffa(8)=YqCoeff0(2,2,1,3)*(1-gamma)^Fp1(2,2,1,3)*dt^3;
YqCoeffa(9)=YqCoeff0(3,1,1,3)*(1-gamma)^Fp1(3,1,1,3)*dt^3;


Fp2(1)=Fp1(1,1,2,3);
Fp2(2)=Fp1(1,2,1,3);
Fp2(3)=Fp1(2,1,1,3);
Fp2(4)=Fp1(1,1,3,3);
Fp2(5)=Fp1(1,2,2,3);
Fp2(6)=Fp1(2,1,2,3);
Fp2(7)=Fp1(1,3,1,3);
Fp2(8)=Fp1(2,2,1,3);
Fp2(9)=Fp1(3,1,1,3);


wVol2dt0=0;
dwVol2dt(1:ExpnOrder)=0.0;

wVol2dt=0;
dwVol2dtdw(1:ExpnOrder)=0.0;

for mm=1:NoOfTerms

    wVol2dt0=YqCoeffa(mm).*w0.^Fp2(mm);
    for nn=1:ExpnOrder
        if(nn==1)
            dwVol2dt(nn)=wVol2dt0*Fp2(mm)*1/w0;
        else
        dwVol2dt(nn)=dwVol2dt(nn-1)*(Fp2(mm)-(nn-1))/w0;
        end
    end
    wVol2dt=wVol2dt+wVol2dt0;
    for nn=1:ExpnOrder
        dwVol2dtdw(nn)=dwVol2dtdw(nn)+dwVol2dt(nn);
    end
        
    
end



end
.
.
.
function [Moments,k] = CalculateVolTermMoments(a0,a,b0,b,SeriesOrder,NMoments)


[EA,EB,EAB] = CalculateSeriesProducts(a0,a,b0,b,SeriesOrder,NMoments);


c0=0;
c(1)=1;
c(2)=0;
d0=-1;
d(1)=0;
d(2)=1;

[EZ1,EZ2,EZ1Z2] = CalculateSeriesProducts(c0,c,d0,d,2,8);


Moments(1)=EA(1)*EZ1(1)+EB(1)*EZ2(1);
Moments(2)=EA(2)*EZ1(2)+2*EAB(1,1)*EZ1Z2(1,1)+EB(2)*EZ2(2);
Moments(3)=EA(3)*EZ1(3)+3*EAB(2,1)*EZ1Z2(2,1)+3*EAB(1,2)*EZ1Z2(1,2)+EB(3)*EZ2(3);
Moments(4)=EA(4)*EZ1(4)+4*EAB(3,1)*EZ1Z2(3,1)+6*EAB(2,2)*EZ1Z2(2,2)+4*EAB(1,3)*EZ1Z2(1,3)+EB(4)*EZ2(4);
Moments(5)=EA(5)*EZ1(5)+5*EAB(4,1)*EZ1Z2(4,1)+10*EAB(3,2)*EZ1Z2(3,2)+10*EAB(2,3)*EZ1Z2(2,3)+5*EAB(1,4)*EZ1Z2(1,4)+EB(5)*EZ2(5);
Moments(6)=EA(6)*EZ1(6)+6*EAB(5,1)*EZ1Z2(5,1)+15*EAB(4,2)*EZ1Z2(4,2)+20*EAB(3,3)*EZ1Z2(3,3)+15*EAB(2,4)*EZ1Z2(2,4)+6*EAB(1,5)*EZ1Z2(1,5)+EB(6)*EZ2(6);
Moments(7)=EA(7)*EZ1(7)+7*EAB(6,1)*EZ1Z2(6,1)+21*EAB(5,2)*EZ1Z2(5,2)+35*EAB(4,3)*EZ1Z2(4,3)+35*EAB(3,4)*EZ1Z2(3,4)+21*EAB(2,5)*EZ1Z2(2,5)+7*EAB(1,6)*EZ1Z2(1,6)+EB(7)*EZ2(7);
Moments(8)=EA(8)*EZ1(8)+8*EAB(7,1)*EZ1Z2(7,1)+28*EAB(6,2)*EZ1Z2(6,2)+56*EAB(5,3)*EZ1Z2(5,3)+70*EAB(4,4)*EZ1Z2(4,4)+56*EAB(3,5)*EZ1Z2(3,5)+28*EAB(2,6)*EZ1Z2(2,6)+8*EAB(1,7)*EZ1Z2(1,7)+EB(8)*EZ2(8);



k(1)=0;
k(2)=Moments(2);
k(3)=Moments(3);
k(4)=Moments(4)-3*Moments(2).^2;
k(5)=Moments(5)-10*Moments(3)*Moments(2);
k(6)=Moments(6)-15*Moments(4)*Moments(2)-10*Moments(3).^2+30*Moments(2).^3;
k(7)=Moments(7)- 21* Moments(2) *Moments(5)- 35* Moments(3) *Moments(4) + 210* Moments(2).^2* Moments(3);
k(8)=Moments(8)- 7 *k(2)* Moments(6) - 21* k(3)* Moments(5) - 35* k(4)* Moments(4) - ... 
  35* k(5)* Moments(3) - 21* k(6)* Moments(2);


end

.
.
.
%function [u1,u2,u3,u4,u5,u6] = CalculateCumulants(a0,a,SeriesOrder,NMoments)
function [CC] = CalculateCumulants(a0,a,SeriesOrder,NMoments)
%[CC] = CalculateCumulants(a0,a,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;

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

%u1=EXZ(2,1);

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


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

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;

CC(1)=k1;
CC(2)=k2;
CC(3)=k3;
CC(4)=k4;

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;
CC(5)=k5;
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;
CC(6)=k6;
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;
CC(7)=k7;
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;
CC(8)=k8;
end


end


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


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

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


for mm=1:NMoments
    if(mm>1)
        %[b0,b] =SeriesProductLogarithmic(a0,a,b0,b,SeriesOrder*mm);
        
        [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*mm);
        
        %b0-bb0
        
        %b-bb
        
        %str=input('Look at two products')
        
        
        b(SeriesOrder*mm+1:SeriesOrder*NMoments+1)=0;
    end
   % Logb=log(abs(b));
   % Signb=sign(b);
    EXZ(mm,1)=b0;
    for pp2=1:SeriesOrder*mm
        if rem(pp2,2)==0
            %EXZ(mm,1)=EXZ(mm,1)+exp(Logb(pp2)+LogEZ(pp2)).*Signb(pp2);
            EXZ(mm,1)=EXZ(mm,1)+b(pp2).*EZ(pp2);
 %           b(pp2).*EZ(pp2)
 %           exp(Logb(pp2)+LogEZ(pp2)).*Signb(pp2)
 %           mm
 %           pp2
            %str=input('Look at moment values')
        end
    end
end

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



end

.
.
.
function [b0,b] = CalculateDensityFromStandardizedMoments(Moments0,b0Guess,bGuess)

% cmu(1)=0;
% cmu(2)=1;
% cmu(3)=Moments0(3)/Moments0(2).^1.5;
% cmu(4)=Moments0(4)/Moments0(2).^2.0;
% cmu(5)=Moments0(5)/Moments0(2).^2.50;
% cmu(6)=Moments0(6)/Moments0(2).^3.0;
% cmu(7)=Moments0(7)/Moments0(2).^3.50;
% cmu(8)=Moments0(8)/Moments0(2).^4.0;

cmu=Moments0;

%[c0,c] =PreSmoothingGuess02(cmu);
%cc=c;
%cc0=-(cc(2)+cc(4)*3+cc(6)*15);

%cc=bGuess;
%cc0=b0Guess;

%[cc0,cc] =PreSmoothingGuess02(cmu);
[cc0,cc] =PreSmoothingGuess02NewFromGuess(cmu,bGuess);

SeriesOrder=7;
NMoments=8;
NZterms=8;




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

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



%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)
[F,dF] = CalculateMomentsAndDerivatives_0(cmu,cc0,cc,SeriesOrder,NZterms,NMoments);
da(1,1)=cc0;
da(2:SeriesOrder+1,1)=cc(1:SeriesOrder);



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


    
   ObjBest=100000*(abs(cmu(1)-Moments(1)))+abs(cmu(2)-Moments(2))+abs((cmu(3)-Moments(3))^(1/1.5))+abs((cmu(4)-Moments(4))^(1/2.0)) + ...
    abs((cmu(5)-Moments(5))^(1/2.0))+abs((cmu(6)-Moments(6))^(1/2.0))+abs((cmu(7)-Moments(7))^(1/2.0))+ ...
    abs((cmu(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)|| (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(cmu,b0,b,SeriesOrder,NZterms,NMoments);
[IsValidFlag] = CheckIsValidDensity(b0,b);

[Moments] = CalculateMomentsOfZSeries(b0,b,SeriesOrder,NMoments);
ObjNew=0;

    ObjNew=100000*(abs(cmu(1)-Moments(1)))+abs(cmu(2)-Moments(2))+abs((cmu(3)-Moments(3))^(1/1.5))+abs((cmu(4)-Moments(4))^(1/2.0)) + ...
    abs((cmu(5)-Moments(5))^(1/2.0))+abs((cmu(6)-Moments(6))^(1/2.0))+abs((cmu(7)-Moments(7))^(1/2.0))+ ...
    abs((cmu(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.
cmu
%str=input('Look at comparison of moments calibration--00');
%b0
%b
nn
%str=input('Look at comparison of moments calibration--0');


b0=b0*sqrt(Moments0(2));
b=b*sqrt(Moments0(2));
b0;
b;
%str=input('Look at comparison of moments calibration--1');
end

.
.
.
function [XDensity,IndexOut,IndexMax] = MakeDensityFromSimulation_Infiniti_NEW(X,Paths,NoOfBins,MaxCutOff )
%Processes monte carlo paths to return a series Xdensity as a function of IndexOut. IndexMax is the maximum value of index.
% 


Xmin=0;
Xmax=0;
for p=1:Paths
if(X(p)>MaxCutOff)
X(p)=MaxCutOff;
end
%if(X(p)<0)
%X(p)=0;
%end
if(Xmin>real(X(p)))
Xmin=real(X(p));
end
if(Xmax<real(X(p)))
Xmax=real(X(p));
end
end

%IndexMax=NoOfBins+1;
BinSize=(Xmax-Xmin)/NoOfBins;
%IndexMax=floor((Xmax-Xmin)/BinSize+.5)+1
IndexMax=floor((Xmax-Xmin)/BinSize+.5)+1

XDensity(1:IndexMax)=0.0;

for p=1:Paths
index=real(floor(real(X(p)-Xmin)/BinSize+.5)+1);
if(real(index)<1)
index=1;
end
if(real(index)>IndexMax)
index=IndexMax;
end

XDensity(index)=XDensity(index)+1.0/Paths/BinSize;

end

IndexOut(1:IndexMax)=Xmin+(0:(IndexMax-1))*BinSize;

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
 
User avatar
Amin
Topic Author
Posts: 2012
Joined: July 14th, 2002, 3:00 am

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

December 6th, 2022, 12:20 pm

Here is the graph output of the program.
The first graph shows density of the SDE from analytic methods VS from monte carlo.
Second graph shows how good reconstruction of the SDE variable is as a function of Standard normal when calculated from Z-series in original coordinates and it is compared with SDE variable from direct transformation of density in Bessel coordinates.
.
.
Image





Image
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
 
User avatar
Amin
Topic Author
Posts: 2012
Joined: July 14th, 2002, 3:00 am

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

December 6th, 2022, 12:46 pm

Sorry Missed a few matlab functions needed to run the program. Here they are:
If there is still any function missing, please mention and I will post it. Most of the functions (except a few new ones) have already been posted in the thread earlier and you can try finding them in the thread.
.
.
.
function [c0,c] = PreSmoothingGuess02NewFromGuess(cmu,cin)


Mul=1.0;
SeriesOrder=7;
NMoments=8;

c0=0;
c(1)=1.0;

%At this point, we have set first two coefficients of Z-series. Rest of the coefficients are zero.    
    
if(cmu(3)==0)
    c(2)=cin(2);
else
    if(cmu(3)>0)
        c(2)=.001;
    elseif(cmu(3)<0)
        c(2)=-.001;
    end
end


for nn=1:5*Mul
    [Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
    
    c(2)=c(2)*(cmu(3))/Moments(3);
    c0=-c(2);
end
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
c=c/sqrt(Moments(2));
c0=c0/sqrt(Moments(2));
%Above we have standardized the coefficients so that 2nd central moment is
%one.
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
Moments2=Moments;
%At this point first three coefficients have been set and the rest are
%zero.
%Moments2 above has been calculated with only c0, c1, c2 non-zero. All
%higher coefficients upto c7 are zero. Particularly in its calculation c3
%is zero.


%if(cmu(4)-Moments2(4)>0)
%    c(3)=.0001;
%elseif(cmu(4)-Moments2(4)<0)
%    c(3)=-.0001;
%end

c(3)=cin(3);

for nn=1:5*Mul
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%please recall that in calculation of Moments2, c3 is zero. We are taking
%slope for fourth moment from where the value of c3 is zero.
Ratio4=(cmu(4)-Moments2(4))/(Moments(4)-Moments2(4));
c(3)=c(3)*Ratio4;


if(cmu(4)<Moments2(4))
    c(3)=abs(c(3))*-1;
end
 %cmu(4)
 %Moments2(4)
 %Moments(4)
 %Ratio4
 %c(3)
% str=input('Look at fourth moment');
% 


end
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%Moments


%str=input('Look at moments--inside my function---M4');
c=c/sqrt(Moments(2));
c0=c0/sqrt(Moments(2));
%Above we have standardized the coefficients so that 2nd central moment is
%one.
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
Moments2=Moments;
%Moments2 above has been calculated with only c0, c1, c2,c3 non-zero. All
%higher coefficients upto c7 are zero. Particularly in its calculation c4
%is zero from where we will calculate slope in next loop..
%if(cmu(5)-Moments2(5)>0)
%    c(4)=.0001;
%elseif(cmu(5)-Moments2(5)<0)
%    c(4)=-.0001;
%end
    
c(4)=cin(4);

for nn=1:5*Mul
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%please recall that in calculation of Moments2, c4 is zero. We are taking
%slope for fifth moment from where the value of c4 is zero.
Ratio5=(cmu(5)-Moments2(5))/(Moments(5)-Moments2(5));
c(4)=c(4)*Ratio5;
c0=-c(2)-3*c(4);





end
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
c=c/sqrt(Moments(2));
c0=c0/sqrt(Moments(2));
%Above we have standardized the coefficients so that 2nd central moment is
%one.    
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);    
Moments2=Moments;    
%Moments2 above has been calculated with only c0, c1, c2,c3,c4 non-zero. All
%higher coefficients upto c7 are zero. Particularly in its calculation c5
%is zero from where we will calculate slope in next loop..    

%if(cmu(6)-Moments2(6)>0)
%    c(5)=.0001;
%elseif(cmu(6)-Moments2(6)<0)
%    c(5)=-.0001;
%end


c(5)=cin(5);

for nn=1:5*Mul
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%please recall that in calculation of Moments2, c5 is zero. We are taking
%slope for sixth moment from where the value of c5 is zero.
Ratio6=(cmu(6)-Moments2(6))/(Moments(6)-Moments2(6));
c(5)=c(5)*Ratio6;

end
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
c=c/sqrt(Moments(2));
c0=c0/sqrt(Moments(2));
%Above we have standardized the coefficients so that 2nd central moment is
%one.
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
Moments2=Moments;
%Moments2 above has been calculated with only c0, c1, c2,c3,c4,c5 non-zero. All
%higher coefficients upto c7 are zero. Particularly in its calculation c6
%is zero from where we will calculate slope in next loop..   

% c(6)=0;
% if(cmu(7)-Moments2(7)>0)
%     c(6)=.0001;
% 
% elseif(cmu(7)-Moments2(7)<0)
%     c(6)=-.0001;
% end


c(6)=cin(6);
    
for nn=1:5*Mul
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%please recall that in calculation of Moments2, c6 is zero. We are taking
%slope for seventh moment from where the value of c6 is zero.
Ratio7=(cmu(7)-Moments2(7))/(Moments(7)-Moments2(7));
c(6)=c(6)*Ratio7;
c0=-c(2)-3*c(4)-15*c(6);

if(cmu(6)<Moments2(6))
    c(5)=abs(c(5))*-1;
end
% cmu(6)
% Moments2(6)
% Moments(6)
% Ratio6
% c(5)
% str=input('Look at sixth moment');


end
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
c=c/sqrt(Moments(2));
c0=c0/sqrt(Moments(2));
%Above we have standardized the coefficients so that 2nd central moment is
%one.
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
Moments2=Moments;
%Moments2 above has been calculated with only c0, c1, c2,c3,c4,c5,c6 non-zero. All
%higher coefficients upto c7 are zero. Particularly in its calculation c7
%is zero from where we will calculate slope in next loop..
% c(7)=0;
% if(cmu(8)-Moments2(8)>0)
%     c(7)=.00005;
% 
% elseif(cmu(8)-Moments2(8)<0)
%     c(7)=-.00005;
% end


c(7)=cin(7);



for nn=1:5*Mul
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%please recall that in calculation of Moments2, c7 is zero. We are taking
%slope for eighth moment from where the value of c7 is zero.

Ratio8=(cmu(8)-Moments2(8))/(Moments(8)-Moments2(8));
c(7)=c(7)*(abs(Ratio8))*sign(Ratio8);
if(cmu(8)<Moments2(8))
    c(7)=abs(c(7))*-1;

% cmu(8)
% Moments2(8)
% Moments(8)
% Ratio8
% c(7)
% str=input('Look at eighth moment');
end
[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
c=c/sqrt(Moments(2));
c0=c0/sqrt(Moments(2));
%Above we have standardized the coefficients so that 2nd central moment is
%one.
end

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

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



%aa0=a0;
if(NMoments>8)
a0=-(a(2)+3*a(4)+15*a(6)+105*a(8));% ---1
end
if(NMoments<=8)
a0=-(a(2)+3*a(4)+15*a(6));% ---1
end

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

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


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

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

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

u1=a0+a(2);
if(SeriesOrder>=4)
u1=a0+a(2)+3*a(4);
end
if(SeriesOrder>=6)
u1=u1+15*a(6);
end
if(SeriesOrder>=8)
u1=u1+105*a(8);
end


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


du1(1)=1;%----2
% du1(2)=0;%----2
% du1(3)=1;%----2
% du1(4)=0;%----2
% du1(5)=1;%----2
% du1(6)=0;%----2
% du1(7)=15;%----2
% du1(8)=0;%----2
% du1(9)=105;%----2
% du1(10)=0;%----2
 du2(1)=2*EXZ(2,1);
 du3(1)=3*EXZ(3,1);
 du4(1)=4*EXZ(4,1);


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

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

end


if(NMoments>=5)
u5=EXZ(6,1);


%du5(1)=5*EXZ(5,1);

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

end

if(NMoments>=6)
u6=EXZ(7,1);
%du6(1)=0;

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

end


if(NMoments>=7)
u7=EXZ(8,1);
%str=input('Look at k7 and k71')

%du7(1)=0;

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

end
if(NMoments>=8)

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

end

if(NMoments>=9)

u9=EXZ(10,1);
for mm=1:SeriesOrder+1
du9(mm)=9*EXZ(9,mm);
end

end
if(NMoments>=10)

u10=EXZ(11,1);
for mm=1:SeriesOrder+1
du10(mm)=10*EXZ(10,mm);
end

end




 F(1,1)=u1-Ms(1);
 F(2,1)=u2-Ms(2);
 F(3,1)=u3-Ms(3);
 F(4,1)=u4-Ms(4);
 if(NMoments>=5)
 F(5,1)=u5-Ms(5);
 end
 if(NMoments>=6)
 F(6,1)=u6-Ms(6);
 end
 if(NMoments>=7)
 F(7,1)=u7-Ms(7);
 end
 if(NMoments>=8)
 F(8,1)=u8-Ms(8);
 end
 if(NMoments>=9)
 F(9,1)=u9-Ms(9);
 end
 if(NMoments>=10)
 F(10,1)=u10-Ms(10);
 end
 
     for mm=1:SeriesOrder+1
        dF(1,mm)=du1(mm);
        dF(2,mm)=du2(mm);
        dF(3,mm)=du3(mm);
        dF(4,mm)=du4(mm);
        if(NMoments>=5)
        dF(5,mm)=du5(mm);
        end
        if(NMoments>=6)
        dF(6,mm)=du6(mm);
        end
        if(NMoments>=7)
        dF(7,mm)=du7(mm);
        end
        if(NMoments>=8)
        dF(8,mm)=du8(mm);
        end
        if(NMoments>=9)
        dF(9,mm)=du9(mm);
        end
        if(NMoments>=10)
        dF(10,mm)=du10(mm);
        end
     end

end

.
.
.
function [EA,EB,EAB] = CalculateSeriesProducts(a0,a,b0,b,SeriesOrder,NMoments)


EA(1:NMoments)=0;
EB(1:NMoments)=0;
EAB(1:NMoments,1:NMoments)=0;
EA(1)=a0+a(2);
EB(1)=b0+b(2);
if(SeriesOrder>=4)
EA(1)=EA(1)+3*a(4);
EB(1)=EB(1)+3*b(4);

end
if(SeriesOrder>=6)
EA(1)=EA(1)+15*a(6);
EB(1)=EB(1)+15*b(6);
end


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



aa(1:NMoments,1:SeriesOrder*NMoments)=0;

aa0(1)=a0;
aa(1,1:SeriesOrder)=a(:);

bb(1:NMoments,1:SeriesOrder*NMoments)=0;
bb0(1)=b0;
bb(1,1:SeriesOrder)=b(:);




for mm=1:NMoments-1
        [ab0,ab] =SeriesProductHigher(a0,a,aa0(mm),aa(mm,:),SeriesOrder*(mm+1));
        aa0(mm+1)=ab0;
        aa(mm+1,1:(mm+1)*SeriesOrder)=ab(1:(mm+1)*SeriesOrder);
        
        [ac0,ac] =SeriesProductHigher(b0,b,bb0(mm),bb(mm,:),SeriesOrder*(mm+1));
        bb0(mm+1)=ac0;
        bb(mm+1,1:(mm+1)*SeriesOrder)=ac(1:(mm+1)*SeriesOrder);

        EA(mm+1)=aa0(mm+1);
        EB(mm+1)=bb0(mm+1);
        for pp=1:SeriesOrder*(mm+1)
        
            EA(mm+1)=EA(mm+1)+aa(mm+1,pp).*EZ(pp);
            EB(mm+1)=EB(mm+1)+bb(mm+1,pp).*EZ(pp);
            
        end
end
        %b(SeriesOrder*mm+1:NMoments*SeriesOrder+1)=0;
aa0;
aa;

for mm=1:NMoments
    for nn=1:NMoments-mm
        dd0=aa0(mm);
        dd=aa(mm,1:SeriesOrder*NMoments);
        ee0=bb0(nn);
        ee=bb(nn,1:SeriesOrder*NMoments);
        
        [ab0,ab] =SeriesProductHigher(dd0,dd,ee0,ee,SeriesOrder*(mm+nn));
        AB0(mm,nn)=ab0;
        AB(mm,nn,1:SeriesOrder*(mm+nn))=ab(1:SeriesOrder*(mm+nn));
        
        EAB(mm,nn)=AB0(mm,nn);
        for pp=1:SeriesOrder*(mm+nn)
        
            EAB(mm,nn)=EAB(mm,nn)+AB(mm,nn,pp).*EZ(pp);
        end
        
        
    end
        
end

end
.
.
.
function [c0,c] =SeriesProductHigher(a0,a,b0,b,SeriesOrder)

a(size(a,2)+1:SeriesOrder)=0;
b(size(b,2)+1:SeriesOrder)=0;
c0=a0*b0;
c(1:SeriesOrder)=0.0;
%a
%b
%str=input('Look at numbers');
for nn=1:SeriesOrder
    c(nn)=c(nn)+a0*b(nn);
    c(nn)=c(nn)+b0*a(nn);
    for kk=1:nn-1
        c(nn)=c(nn)+a(kk)*b(nn-kk);
    end
end

end

.
.
.
function [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 [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

.
.
.
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
 
User avatar
Amin
Topic Author
Posts: 2012
Joined: July 14th, 2002, 3:00 am

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

December 6th, 2022, 3:59 pm

Friends, almost everyday when I try to go to sleep, mind control agents leave a pungent like gas in my room. I always keep window of my room open to get fresh air, and as soon as the atmosphere in my room goes from fresh to pungent that feels like slightly acidic as I breathe it, I can right away tell they have released some gas in my room. The gas is so acidic that taste in my mouth changes as I breathe it and I can feel acidity in my mouth. I have been busy for past 4-5 days working on my programs and could not find time to mention it on the forum. But now when mind control agents released gas in my room, I decided to write about it immediately. Please force mind control crooks to stop using ionized, pungent and acidic gases on innocent civilians.
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
 
User avatar
Amin
Topic Author
Posts: 2012
Joined: July 14th, 2002, 3:00 am

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

December 7th, 2022, 9:49 am

Friends, I am travelling tomorrow. I am going from Lahore to Kot Addu where I might stay for a month or slightly more. I want to complete the work for my paper in the meantime. I want to put serious effort into trying to make it a good paper. Kot Addu is a very small city and they can easily drug food at just a few places where I can eat outside. I want to request friends to please ask mind control agencies to not drug any food in the small city. Also an antipsychotic injection is due next week and there is probably just one pharmacy in Kot Addu where I can get the particular injection and it would be very easy for them to give me an injection that is filled with mind control chemicals and my life after that would become extremely difficult. I want to request friends to please keep a vigil and force mind control agencies to not drug the food or injections. I hope that most friends would like if I can present some good research while working without mind control chemicals. 

Another interesting thing I want to do is to try to discover if we can find out Z-series representation of a stochastic variable once its closed form density is given. As opposed to finding Z-series from moments of the random variable. I have some ideas and I am very sure it is possible with a bit of hack. This should be very easy to find out if formulas for CDF of the random variable are given. I will come back in a two to three days with examples of lognormal diffusions, CEV diffusions and other cases where closed form formulas exist and will work out how to calculate their Z-series expansions. 
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
 
User avatar
Amin
Topic Author
Posts: 2012
Joined: July 14th, 2002, 3:00 am

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

December 9th, 2022, 9:40 am

Friends, I travelled from Lahore to Kot Addu yesterday. It is around six hours of travel. It is not very long but long enough to be stimulating. Another possible reason for stimulation for a mind control victim might be that it is difficult to control the victim and project waves on his mind with precision during long travel despite that most major highways and roads have mind control devices installed in Pakistan and Pak army has ensured that mind control waves could be projected when CIA mind control victims are travelling on those roads and highways. Especially when some victim connects his brain when mind control is weak and brain connects to parts unknown to mind control science, it can continue to create problems for mind control agents for several weeks in properly controlling the victim after the victim has travelled on a long journey. So when I have to travel, I like to make it stimulating while mind control agencies like to make it otherwise and sometimes drug the particular beverages like energy and caffeine drinks along the entire highway. 
So something very interesting happened during night before travel. I bought three red bulls (red bull has recently become good in Lahore) and a chicken grilled sandwich the night before travel and kept them in my room so that I would be able to eat and drink good food while travelling next day. I know they can open door bolts that are used to lock the doors from outside and before I sleep, I usually put compressed tissue in the bolt's path so that it could not be opened from outside using magnets or anything. Mostly I put enough tissue that it would be extremely difficult to open the door from outside. But no external intrusion/unlocking had happened for more than five months so I was not being extra careful. I had put compressed tissue beneath the bolt but it was not so super compressed and if enough force could be applied from outside, it would open. I was just not being careful enough since nothing of the sort had happened for months. When I woke up in the morning, compressed tissue was lying on the dressing table on one side(as nobody could place it back from outside using magnets or any devices after locking the door from outside).  I knew they had entered my room while I was sleeping to drug the red bulls and the sandwich. 
We (my parents and I) were supposed to start the travel well into the day since due to high fog motorways are mostly closed in the morning at this time of the year. I had enough time so I went out and bought a few red bulls and bottled water again quickly and came back to the house and I had good food while travelling.
The travel was very stimulating and I thought of several new ideas and gained clarity about many new things related to my research which I will continue to share with wilmott friends in coming weeks.
Today, after coming to Kot addy yesterday evening, I cleaned my room and removed huge amount of dust and then changed sheets and now I really like my clean and big airy room that has two large windows and lots of sunlight.  I will start working on my research tonight and hope that it would be a lot of fun working here in our Kot Addu home. I really enjoy coming to Kot Addu in winters. And drive around on the motorcycle in the morning around the countryside with lush green fields on both sides of the roads or streams. I was able to do very good research when I came here last winters.
I bought bottled water and some beverages from Kot Addu yesterday night and this morning and bottled water and many other drinks in the market are perfectly good. I am very afraid that mind control agencies would start drugging the food, drinks and bottled water in coming days to make life difficult for me. Please ask them to not drug and food or beverages in Kot Addu.
I will start working on my research in a few hours. Working at night is fun in mild winters. And first thing I Want to do is to see if we can easily find Z-series of a random variable when its closed form analytic density is given. And then present some simple examples from lognormal density and other well known densities. I hope some results will be ready in another day or two. So I hope to come back quickly with some new work in another 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
 
User avatar
Amin
Topic Author
Posts: 2012
Joined: July 14th, 2002, 3:00 am

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

December 9th, 2022, 2:31 pm

Friends, mind control agents are getting desperate to drug my food. And they are coming to my room even when I am present inside the house somewhere else. We have a large old house that was built in late eighties. Upper floor also has three bedrooms and a large living room and was constructed a decade later in late nineties. Stairs to upper floor are right next to main entrance. Main entrance (with stairs) leads to a gallery and other rooms and kitchen have approach from the gallery. There is no living room on ground floor. My room is on first floor and only I am living on first floor. It is a large 3.5 kanals house but since nobody lives here (except for my father) everything including lawns remains dark unless there is somebody living in a room for a few days or weeks.
Somewhat after writing previous post, I went out and brought some food. I quickly had some bit of food in my room on first floor and then went down to kitchen to make green tea in microwave before I start working on my research. OK, I very clearly remember this when I was leaving my room, I turned off lights of my room and then also turned off lights outside my room. And there was dark on first floor and the staircase and no light was coming out of my room even though its door was open. I clearly remember that I thought I can easily go down the stairs in the dark but it would be difficult to come up when I would have hot green tea in my hands and then I turned on the light outside my room. I recall this very clearly that there was dark all over the first floor and then I turned on outside light so I would not face difficulty going up the stairs when I come back. I remained in the kitchen for more than fifteen minutes. Anybody could come inside the house from the dark lawns and go to upper floor without my knowing about it at all since kitchen has approach from gallery and you could not view the entrance or the stairs from there. When I came upstairs after making green tea, not only outside light was on, the light of my room was also on which I knew very well that I had turned off. I was very surprised and came down and explicitly asked my brother, father and mother but they all refused that they went upstairs in my absence. My only guess is that somebody sneaked quickly from the entrance and took stairs to my room possibly to drug the food when I was in the kitchen making green tea. I am going to be careful with the food that was left in my room but I want friends to complain to mind control agencies to end their nasty tactics. I really do not understand what urgency they have that they have to resort to such desperate tactics.
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