Serving the Quantitative Finance Community

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

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

March 8th, 2022, 7:20 pm

Sorry friends, I prematurely hit the submit button and so could not complete the previous post. Here is the actual relevant post. Please disregard the previous incomplete post.
Friends, I wanted to share a rough sketch of the proof of the method I presented in above programs. I hope friends would like it. 
I am not writing Fokker-planck equation in Bessel coordinates and I am sure most friends are familiar with it. 
In our setting we have an SDE variable, [$]B[$], in Bessel coordinates on nth grid point [$]Z_n[$] of underlying Z-grid. This is represented as [$]B(Z_n)[$]. As a function of time, [$]B(Z_n)[$] is moving so as that CDF of the density at point [$]B(Z_n)[$] remains constant. This is written in equation form as below

[$]\frac{\partial}{\partial t} \big[\int_{-\infty}^{B(Z_n)} P(B) dB \big]=0[$]

applying the time derivative to terms in brackets, we get

[$]\int_{-\infty}^{B(Z_n)} \frac{\partial P(B)}{\partial t} dB +  \frac{\partial B}{\partial t} P(B(Z_n)) =0[$]

I will make it more formal but here, since we add drift separately and because we are in Bessel coordinates, I am just replacing [$]\frac{\partial P(B)}{\partial t}[$] by [$].5 {\sigma}^2 \frac{\partial^2 P(B)}{\partial B^2}[$] in the above equation.
Making the substitution in above equation, we get

[$]\int_{-\infty}^{B(Z_n)} \big[.5 {\sigma}^2 \frac{\partial^2 P(B)}{\partial B^2} \big] dB +  \frac{\partial B}{\partial t} P(B(Z_n)) =0[$]

Applying the integral on first term, we get

[$]\big[.5 {\sigma}^2 \frac{\partial P(B(Z_n))}{\partial B} \big] +  \frac{\partial B}{\partial t} P(B(Z_n)) =0[$]

Now we want to make a change from [$]P(B)[$] to [$]P(Z)[$]. In order to do that we have two equations
[$]P(B)=P(Z)\frac{\partial Z}{\partial B}[$] and
[$]\frac{\partial P(B)}{\partial B}=\frac{\partial P(Z)}{\partial Z} {(\frac{\partial Z}{\partial B})}^2 + P(Z) \frac{\partial^2 Z}{\partial B^2}[$]
Now substituting above two equations in previous equation, we get
[$].5 {\sigma}^2 \big[\frac{\partial P(Z_n)}{\partial Z} {(\frac{\partial Z}{\partial B})}^2 + P(Z_n) \frac{\partial^2 Z}{\partial B^2} \big] + \frac{\partial B}{\partial t} P(Z_n) \frac{\partial Z}{\partial B}=0 [$]
substituting [$]\frac{\partial^2 Z}{\partial B^2} =-{(\frac{\partial Z}{\partial B})}^3 \frac{\partial^2 B}{\partial Z^2}[$]
we get
[$].5 {\sigma}^2 \big[\frac{\partial P(Z_n)}{\partial Z} {(\frac{\partial Z}{\partial B})}^2 - P(Z_n) {(\frac{\partial Z}{\partial B})}^3 \frac{\partial^2 B}{\partial Z^2} \big] + \frac{\partial B}{\partial t} P(Z_n) \frac{\partial Z}{\partial B}=0 [$]
cancelling [$]\frac{\partial Z}{\partial B}[$] on both sides, we get
[$].5 {\sigma}^2 \big[\frac{\partial P(Z_n)}{\partial Z} {(\frac{\partial Z}{\partial B})} - P(Z_n) {(\frac{\partial Z}{\partial B})}^2 \frac{\partial^2 B}{\partial Z^2} \big] +\frac{\partial B}{\partial t} P(Z_n)=0 [$]
The first term in above is expansion that we have treated in the previous matlab program in a different analytic way by adding squared volatilities . Matching the coefficients of second and third term, we get

[$]\frac{\partial B}{\partial t} = .5 {\sigma}^2 \big[{(\frac{\partial Z}{\partial B})}^2 \frac{\partial^2 B}{\partial Z^2} \big][$]

the last equation is the analytic solution that I have used in my previous matlab programs to get the density of SDEs in Bessel coordinates right.
.
.
.
Here I want to explain how the above analytics could be applied to Fokker-planck equation in original coordinates (and it works seamlessly). 
First I write the fokker-planck equation in original coordinates as
[$]\frac{\partial P(x,t)}{\partial t} = -\frac{\partial \big[\mu(x)P(x,t) \big]}{\partial x}\, +.5 \, \frac{\partial^2 \big[{\sigma(x)}^2P(x,t) \big]}{\partial x^2}[$]

We have our partial differential equation on a grid and we want to determine the movement of an arbitrary grid point (boundary) along time so that mass within(total mass up till the grid point is conserved or associated CDF remains constant) that grid point (boundary) remains constant as a function of time. we denote this arbitrary grid point as [$]x_b[$].
We write the conservation of probability mass up till (or associated constant CDF) the grid point (boundary) as a function of time in equation form as

[$]\frac{\partial}{\partial t} \big[\int_{-\infty}^{x_b} P(x,t) dx \big]=0[$]

applying the time derivative to terms in brackets, we get

[$]\int_{-\infty}^{x_b} \frac{\partial P(x,t)}{\partial t} dx +  \frac{\partial x_b}{\partial t} P(x_b,t) =0[$]

But we know that term inside the integral is given by FP equation and is LHS of FP equation. We replace the time derivative inside the integral  by RHS of FP equation as

[$]\int_{-\infty}^{x_b} \big[   -\frac{\partial [\mu(x)P(x,t)]}{\partial x}\, + \, .5\frac{\partial^2 \big[{\sigma(x)}^2P(x,t) \big]}{\partial x^2}\big] dx +  \frac{\partial x_b}{\partial t} P(x_b,t) =0[$]

Applying the integral to complete differentials inside square brackets, we get

[$]-\mu(x_b)P(x_b,t)\, + \, .5 \, \frac{\partial \big[{\sigma(x_b)}^2P(x_b,t) \big]}{\partial x} +  \frac{\partial x_b}{\partial t} P(x_b,t) =0[$]
[$]=-\mu(x_b)P(x_b,t)\, + .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}P(x_b,t) + .5\,{\sigma(x_b)}^2 \, \frac{\partial P(x_b,t)}{\partial x} +\,  \frac{\partial x_b}{\partial t} P(x_b,t) =0[$]

We have solved the equation in original coordinates and we could do an ODE in [$]\frac{\partial x_b}{\partial t}[$] to find the evolution of every grid point so that probability mass till that grid point is conserved. But it will be a lot easier if we converted the above arbitrary density to underlying gaussian density as we have been doing and for that we write the following substitution equations from change of probability derivatives as([$]Z_b[$] is associated with each point [$]x_b[$] using common CDF as probability is being conserved.)

[$]P(x)=P(Z)\frac{\partial Z}{\partial x}[$] and

[$]\frac{\partial P(x)}{\partial x}=\frac{\partial P(Z)}{\partial Z} {(\frac{\partial Z}{\partial x})}^2 + P(Z) \frac{\partial^2 Z}{\partial x^2}[$]

Substituting the above two equations in previous equation, we get

[$]-\mu(x_b) \, P(Z_b)\, \frac{\partial Z}{\partial x} \,+ .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}P(Z_b) \frac{\partial Z}{\partial x} + .5\,{\sigma(x_b)}^2 \, \big[ \frac{\partial P(Z_b)}{\partial Z} {(\frac{\partial Z}{\partial x})}^2 + P(Z_b) \frac{\partial^2 Z}{\partial x^2}    \big] +\,  \frac{\partial x_b}{\partial t} P(Z_b) \frac{\partial Z}{\partial x}=0[$]

We make the following substitutions in above equation

[$]\frac{\partial^2 Z}{\partial x^2} =-{(\frac{\partial Z}{\partial x})}^3 \frac{\partial^2 x}{\partial Z^2}[$]

and

[$]\frac{\partial P(Z_b)}{\partial Z}=-Z_b \, P(Z_b)[$]


[$]-\mu(x_b) \, P(Z_b)\, \frac{\partial Z}{\partial x} \,+ .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}P(Z_b) \frac{\partial Z}{\partial x} + .5\,{\sigma(x_b)}^2 \, (-Z_b) \, P(Z_b) {(\frac{\partial Z}{\partial x})}^2 - .5\,{\sigma(x_b)}^2 P(Z_b) {(\frac{\partial Z}{\partial x})}^3 \frac{\partial^2 x}{\partial Z^2}  +\,  \frac{\partial x_b}{\partial t} P(Z_b) \frac{\partial Z}{\partial x}=0[$]

Cancelling  [$]P(Z_b) \frac{\partial Z}{\partial x}[$] throughout the equation and rearranging, we get the first order very simple ODE for our grid point that conserves the probability mass up till that grid point.

[$]  \frac{\partial x_b}{\partial t} =\mu(x_b) - .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}+ .5\,{\sigma(x_b)}^2 \, Z_b \, \frac{\partial Z}{\partial x} + \, .5\,{\sigma(x_b)}^2 \, {(\frac{\partial Z}{\partial x})}^2 \frac{\partial^2 x}{\partial Z^2}  [$]

The above is the first order ODE (that has to be solved in time) for the grid point with associated CDF or conserved mass and this point seamlessly moves in time as a function of above ODE when we solve it. I was able to solve the above ODE just as such (without any squared addition of volatilities separately) and got seamless results for the evolution of CDF(or probability mass up till  that point) conservation point in original coordinates.
When you evolve the SDE in original coordinates, you have to take a few steps with an alternative method (until second derivative is non-negligible) and then evolve with above method. I will be posting a program of simulation of densities of SDEs with above method in a few hours.

The above method is sister equation to continuity equation of mathematical physics. You can apply it to other first order equations or higher order equations. I am sure it can also be applied to Navier Stokes equation with some hack when some boundary is in motion and momentum is conserved. Continuity equation works when boundaries are given while this method works when boundary has to be solved.
.
.
Friends, below is the equation for evolution of SDE variable with constant CDF in original coordinates 
[$]  \frac{\partial x_b}{\partial t} =\mu(x_b) - .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}+ .5\,{\sigma(x_b)}^2 \, Z_b \, \frac{\partial Z}{\partial x} + \, .5\,{\sigma(x_b)}^2 \, {(\frac{\partial Z}{\partial x})}^2 \frac{\partial^2 x}{\partial Z^2}  [$]
Its Bessel coordinates version is given as
[$]  \frac{\partial x_b}{\partial t} =\mu(x_b) + .5\,{\sigma}^2 \, Z_b \, \frac{\partial Z}{\partial x} + \, .5\,{\sigma}^2 \, {(\frac{\partial Z}{\partial x})}^2 \frac{\partial^2 x}{\partial Z^2}  [$]
which is equivalent to 
[$]  \frac{\partial x_b}{\partial t} =\mu(x_b) + .5\,{\sigma}^2 \, Z_b \, \frac{\partial Z}{\partial x} - \, .5\,{\sigma}^2 \, \frac{\partial x}{\partial Z} \frac{\partial^2 Z}{\partial x^2}  [$]

I was playing with the above equation and realized that the above equation can be written in a surprisingly simple and very familiar form in a simple (partial)differential equation as

[$]  \frac{\partial x}{\partial t} =\mu(x) + .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \big[.5 Z^2-\,  \log \frac{\partial Z}{\partial x} \big]  [$]

or equivalently
[$]  \frac{\partial x}{\partial t}=\mu(x) + .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \big[.5 Z^2+\,  \log \frac{\partial x}{\partial Z} \big]   [$]

In the above equations I have removed subscript b that was meant to denote that we are dealing with a certain value of x as a function of Z.
I am writing these simpler equations in the hope that friends who are experts at differential equations might be able to find a simple analytical solution to above differential equation. It is also nice to see how these simplified equations intuitively relate to standard normal. I will also try to see on my end if this simplification could result in simpler and more robust methods to solve the SDEs and how these could be generalized to higher dimensions.

Interesting to note that the solution to evolution equation in original coordinates can also be given easily as
[$]  \frac{\partial x_b}{\partial t}=\mu(x_b) - .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}+ .5\,{\sigma(x_b)}^2 \, \frac{\partial }{\partial x} \big[.5 Z^2+\,  \log \frac{\partial x}{\partial Z} \big]   [$]
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: 1839
Joined: July 14th, 2002, 3:00 am

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

March 8th, 2022, 8:07 pm

.
.
Friends, below is the equation for evolution of SDE variable with constant CDF in original coordinates 
[$]  \frac{\partial x_b}{\partial t} =\mu(x_b) - .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}+ .5\,{\sigma(x_b)}^2 \, Z_b \, \frac{\partial Z}{\partial x} + \, .5\,{\sigma(x_b)}^2 \, {(\frac{\partial Z}{\partial x})}^2 \frac{\partial^2 x}{\partial Z^2}  [$]
Its Bessel coordinates version is given as
[$]  \frac{\partial x_b}{\partial t} =\mu(x_b) + .5\,{\sigma}^2 \, Z_b \, \frac{\partial Z}{\partial x} + \, .5\,{\sigma}^2 \, {(\frac{\partial Z}{\partial x})}^2 \frac{\partial^2 x}{\partial Z^2}  [$]
which is equivalent to 
[$]  \frac{\partial x_b}{\partial t} =\mu(x_b) + .5\,{\sigma}^2 \, Z_b \, \frac{\partial Z}{\partial x} - \, .5\,{\sigma}^2 \, \frac{\partial x}{\partial Z} \frac{\partial^2 Z}{\partial x^2}  [$]

I was playing with the above equation and realized that the above equation can be written in a surprisingly simple and very familiar form in a simple (partial)differential equation as

[$]  \frac{\partial x}{\partial t} =\mu(x) + .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \big[.5 Z^2-\,  \log \frac{\partial Z}{\partial x} \big]  [$]

or equivalently
[$]  \frac{\partial x}{\partial t}=\mu(x) + .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \big[.5 Z^2+\,  \log \frac{\partial x}{\partial Z} \big]   [$]

In the above equations I have removed subscript b that was meant to denote that we are dealing with a certain value of x as a function of Z.
I am writing these simpler equations in the hope that friends who are experts at differential equations might be able to find a simple analytical solution to above differential equation. It is also nice to see how these simplified equations intuitively relate to standard normal. I will also try to see on my end if this simplification could result in simpler and more robust methods to solve the SDEs and how these could be generalized to higher dimensions.

Interesting to note that the solution to evolution equation in original coordinates can also be given easily as
[$]  \frac{\partial x_b}{\partial t}=\mu(x_b) - .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}+ .5\,{\sigma(x_b)}^2 \, \frac{\partial }{\partial x} \big[.5 Z^2+\,  \log \frac{\partial x}{\partial Z} \big]   [$]
.
.
Taking this a bit further we have from above equation
[$]  \frac{\partial x}{\partial t} =\mu(x) + .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \big[.5 Z^2-\,  \log \frac{\partial Z}{\partial x} \big]  [$]
[$]  =\mu(x) - .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \big[-.5 Z^2+\,  \log \frac{\partial Z}{\partial x} \big]  [$]
[$]  =\mu(x) - .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \Big[\,  \log \big[ \exp(-.5 Z^2) \frac{\partial Z}{\partial x} \big]\Big]  [$]

derivative of log-likelihood multiplied by variance??
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: 1839
Joined: July 14th, 2002, 3:00 am

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

March 8th, 2022, 9:00 pm

.
.
Friends, below is the equation for evolution of SDE variable with constant CDF in original coordinates 
[$]  \frac{\partial x_b}{\partial t} =\mu(x_b) - .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}+ .5\,{\sigma(x_b)}^2 \, Z_b \, \frac{\partial Z}{\partial x} + \, .5\,{\sigma(x_b)}^2 \, {(\frac{\partial Z}{\partial x})}^2 \frac{\partial^2 x}{\partial Z^2}  [$]
Its Bessel coordinates version is given as
[$]  \frac{\partial x_b}{\partial t} =\mu(x_b) + .5\,{\sigma}^2 \, Z_b \, \frac{\partial Z}{\partial x} + \, .5\,{\sigma}^2 \, {(\frac{\partial Z}{\partial x})}^2 \frac{\partial^2 x}{\partial Z^2}  [$]
which is equivalent to 
[$]  \frac{\partial x_b}{\partial t} =\mu(x_b) + .5\,{\sigma}^2 \, Z_b \, \frac{\partial Z}{\partial x} - \, .5\,{\sigma}^2 \, \frac{\partial x}{\partial Z} \frac{\partial^2 Z}{\partial x^2}  [$]

I was playing with the above equation and realized that the above equation can be written in a surprisingly simple and very familiar form in a simple (partial)differential equation as

[$]  \frac{\partial x}{\partial t} =\mu(x) + .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \big[.5 Z^2-\,  \log \frac{\partial Z}{\partial x} \big]  [$]

or equivalently
[$]  \frac{\partial x}{\partial t}=\mu(x) + .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \big[.5 Z^2+\,  \log \frac{\partial x}{\partial Z} \big]   [$]

In the above equations I have removed subscript b that was meant to denote that we are dealing with a certain value of x as a function of Z.
I am writing these simpler equations in the hope that friends who are experts at differential equations might be able to find a simple analytical solution to above differential equation. It is also nice to see how these simplified equations intuitively relate to standard normal. I will also try to see on my end if this simplification could result in simpler and more robust methods to solve the SDEs and how these could be generalized to higher dimensions.

Interesting to note that the solution to evolution equation in original coordinates can also be given easily as
[$]  \frac{\partial x_b}{\partial t}=\mu(x_b) - .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}+ .5\,{\sigma(x_b)}^2 \, \frac{\partial }{\partial x} \big[.5 Z^2+\,  \log \frac{\partial x}{\partial Z} \big]   [$]
.
.
Taking this a bit further we have from above equation
[$]  \frac{\partial x}{\partial t} =\mu(x) + .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \big[.5 Z^2-\,  \log \frac{\partial Z}{\partial x} \big]  [$]
[$]  =\mu(x) - .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \big[-.5 Z^2+\,  \log \frac{\partial Z}{\partial x} \big]  [$]
[$]  =\mu(x) - .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \Big[\,  \log \big[ \exp(-.5 Z^2) \frac{\partial Z}{\partial x} \big]\Big]  [$]

derivative of log-likelihood multiplied by variance??
.
Taking the above equation
[$]  \mu(x) - .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \Big[\,  \log \big[ \exp(-.5 Z^2) \frac{\partial Z}{\partial x} \big]\Big]  [$]
but using change of densities formula (standard normal to SDE variable x)
[$]   \exp(-.5 Z^2) \frac{\partial Z}{\partial x} =\sqrt{2 \pi} p(x)  [$]

In the above equation p(x) means probability density function of SDE in bessel (or original) coordinates
which means

[$]  =\mu(x) - .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \Big[\,  \log \big[ \exp(-.5 Z^2) \frac{\partial Z}{\partial x} \big]\Big]  [$]
[$]  =\mu(x) - .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \Big[\,  \log \big[\sqrt{2 \pi} p(x)]\Big]  [$]
[$]  =\mu(x) - .5\,{\sigma}^2 \, \frac{ p'(x)}{p(x)}  [$]

so it seems that constant CDF lines of solution of SDEs density  move in time according to formula
[$]\frac{\partial x_b}{\partial t}  =\mu(x_b) - .5\,{\sigma}^2 \, \frac{ p'(x_b)}{p(x_b)}  [$]

where subscript b indicates that we are that we are solving for a particular point with associated constant CDF.


I am really not sure if the above equation is right and I might possibly have made some mistake but I hope this is interesting for friends. Please pardon if there is any error.
Interesting that the above formula can be derived directly from fokker-planck equation if we do not convert to Z-coordinates. 
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: 1839
Joined: July 14th, 2002, 3:00 am

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

March 8th, 2022, 10:25 pm

.
.
Friends, below is the equation for evolution of SDE variable with constant CDF in original coordinates 
[$]  \frac{\partial x_b}{\partial t} =\mu(x_b) - .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}+ .5\,{\sigma(x_b)}^2 \, Z_b \, \frac{\partial Z}{\partial x} + \, .5\,{\sigma(x_b)}^2 \, {(\frac{\partial Z}{\partial x})}^2 \frac{\partial^2 x}{\partial Z^2}  [$]
Its Bessel coordinates version is given as
[$]  \frac{\partial x_b}{\partial t} =\mu(x_b) + .5\,{\sigma}^2 \, Z_b \, \frac{\partial Z}{\partial x} + \, .5\,{\sigma}^2 \, {(\frac{\partial Z}{\partial x})}^2 \frac{\partial^2 x}{\partial Z^2}  [$]
which is equivalent to 
[$]  \frac{\partial x_b}{\partial t} =\mu(x_b) + .5\,{\sigma}^2 \, Z_b \, \frac{\partial Z}{\partial x} - \, .5\,{\sigma}^2 \, \frac{\partial x}{\partial Z} \frac{\partial^2 Z}{\partial x^2}  [$]

I was playing with the above equation and realized that the above equation can be written in a surprisingly simple and very familiar form in a simple (partial)differential equation as

[$]  \frac{\partial x}{\partial t} =\mu(x) + .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \big[.5 Z^2-\,  \log \frac{\partial Z}{\partial x} \big]  [$]

or equivalently
[$]  \frac{\partial x}{\partial t}=\mu(x) + .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \big[.5 Z^2+\,  \log \frac{\partial x}{\partial Z} \big]   [$]

In the above equations I have removed subscript b that was meant to denote that we are dealing with a certain value of x as a function of Z.
I am writing these simpler equations in the hope that friends who are experts at differential equations might be able to find a simple analytical solution to above differential equation. It is also nice to see how these simplified equations intuitively relate to standard normal. I will also try to see on my end if this simplification could result in simpler and more robust methods to solve the SDEs and how these could be generalized to higher dimensions.

Interesting to note that the solution to evolution equation in original coordinates can also be given easily as
[$]  \frac{\partial x_b}{\partial t}=\mu(x_b) - .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}+ .5\,{\sigma(x_b)}^2 \, \frac{\partial }{\partial x} \big[.5 Z^2+\,  \log \frac{\partial x}{\partial Z} \big]   [$]
.
.
Taking this a bit further we have from above equation
[$]  \frac{\partial x}{\partial t} =\mu(x) + .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \big[.5 Z^2-\,  \log \frac{\partial Z}{\partial x} \big]  [$]
[$]  =\mu(x) - .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \big[-.5 Z^2+\,  \log \frac{\partial Z}{\partial x} \big]  [$]
[$]  =\mu(x) - .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \Big[\,  \log \big[ \exp(-.5 Z^2) \frac{\partial Z}{\partial x} \big]\Big]  [$]

derivative of log-likelihood multiplied by variance??
.
Taking the above equation
[$]  \mu(x) - .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \Big[\,  \log \big[ \exp(-.5 Z^2) \frac{\partial Z}{\partial x} \big]\Big]  [$]
but using change of densities formula (standard normal to SDE variable x)
[$]   \exp(-.5 Z^2) \frac{\partial Z}{\partial x} =\sqrt{2 \pi} p(x)  [$]

In the above equation p(x) means probability density function of SDE in bessel (or original) coordinates
which means

[$]  =\mu(x) - .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \Big[\,  \log \big[ \exp(-.5 Z^2) \frac{\partial Z}{\partial x} \big]\Big]  [$]
[$]  =\mu(x) - .5\,{\sigma}^2 \, \frac{\partial }{\partial x} \Big[\,  \log \big[\sqrt{2 \pi} p(x)]\Big]  [$]
[$]  =\mu(x) - .5\,{\sigma}^2 \, \frac{ p'(x)}{p(x)}  [$]

so it seems that constant CDF lines of solution of SDEs density  move in time according to formula
[$]\frac{\partial x_b}{\partial t}  =\mu(x_b) - .5\,{\sigma}^2 \, \frac{ p'(x_b)}{p(x_b)}  [$]

where subscript b indicates that we are that we are solving for a particular point with associated constant CDF.


I am really not sure if the above equation is right and I might possibly have made some mistake but I hope this is interesting for friends. Please pardon if there is any error.
Interesting that the above formula can be derived directly from fokker-planck equation if we do not convert to Z-coordinates. 
.
.
In order to use this equation.
[$]\frac{\partial x_b}{\partial t}  =\mu(x_b) - .5\,{\sigma}^2 \, \frac{ p'(x_b)}{p(x_b)}  [$]
Since [$]Z_b[$] associated with every constant cdf point on grid is known
use
[$] p(x_b) =1.0/\sqrt{2 \pi}  \exp(-.5 {Z_b}^2) \frac{\partial Z_b}{\partial x_b}    [$]
for above we have to calculate change of derivative with respect to [$]Z_b[$] on our grid.
Once [$] p(x_b) [$] is known, directly differentiate it on non-uniform grid using high order discretization for first derivative to find [$]p'(x_b) [$]
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: 1839
Joined: July 14th, 2002, 3:00 am

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

March 12th, 2022, 4:58 pm

Friends, I drove from Islamabad to Kot Addu today. It took me around eight hours to reach kot Addu. Long drives are very stimulating and I thought of several new ideas both about SDEs and stocks trading. But when I reached Kot Addu, I had extreme uneasiness in my penis since mind control agents continued to charge my penis on the way as I had to be in the same accurately predictable position for many hours sitting in the driving seat. After reaching home, I urinated again and again to feel better but it did not work. I also realized that my urine was reddish as there was blood mixed with urine as I urinated. 
When they force electromagnetic waves inside the penis, the body tissues inside the penis tube get damaged and blood starts to seep through the tissues and gets mixed with urine. Forcing waves into the penis is a very important part of mind control since it is impossible to control adrenaline and dopamine without controlling the sex organ. Adrenaline and dopamine give us drive to concentrate and force ourselves to work hard so mind control agents really have to control these and several other similar neurotransmitters and they have to force waves into our penis to be able to do that. 
I normally do not like to urinate in the wash room since when I urinate there, it is very easy for them to force high frequency waves inside the penis and it feels like a needle is going into the penis when urine is coming out and high frequency waves are going inside. I would rather take a mineral water bottle and urine into it in my room and drop the urine in the wash room and wash the bottle to use again. But if I would remain in a single place in my room to urinate, they install gadgets and urinating in my room is no different from urinating in the wash room as they can again force waves into the penis and the target victim would have needle like feeling going into the penis so I have to continue to change places randomly in my room to urinate. 
Anyway they were doing the same thing again and blood was coming out of my penis mixed with urine as I supposed the inside of tube that carries urine out of our body was damaged due to high frequency waves. This is very serious thing so I decided to write about it so as to deter the mind control crooks who are bringing a great name to American people by such feats. But these hardliner conservative crooks in army will continue their ugly antics since their desire to enter the paradise of Jewish billionaire godfather is stronger than regard for any country or any human ethics.
Previously I had told friends several times (may be eight or nine times in total) that they force chemicals into my mouth when I am sleeping. This is sheer science though I do not know underlying exact details but I know they can create potential at a precise place for special types of ions(neurotransmitters) emitted from elsewhere with lasers and EM waves. Well I continued to mention that they charge my mouth every night with special mind control chemicals that settle on my tongue and inside of upper jaw and my mouth is very bitter and I have to wash it or brush inside of my mouth to feel better. I continued to mention it again and again on this forum but they never stopped so I stopped writing about it.
I want to tell Americans that Americans are mostly not bad people at all but these hardliner conservative crooks in army and mind control who are dying to enter the paradise of the jewish billionaire godfather are not bringing any good name to Americans with their inhuman machinations.
Another thing I have mentioned several times before is that I find likes on several linkedin posts (in my activity) that I never liked. They know my passwords and have enough control over my computer and they repeatedly like something I never liked (I never like anything unless I strongly feel about it which is usually relatively rare). This is something that gets known to me since there is a log of every like in my Linkedin activity. Imagine what else they are doing from my internet traffic that I never came to know.
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: 1839
Joined: July 14th, 2002, 3:00 am

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

March 16th, 2022, 9:24 am

Friends I recently posted a cumulant matching program that advances the density of SDE diffusions (represented in the form of a Z-series) by adding the cumulants of existing density with cumulants of innovations(diffusion term) and then finding another SDE density in Z-series so that its cumulants match the added cumulants of innovations and existing density. This technique and its derivatives can be applied in interesting ways at a lot of places in applied probability theory and I am sure many friends would like to use it extensively. Apart from a lot of simple applied uses, it can be used to solve for example the densities where the driver noise in not normal but for example is a gamma density variable.   
However backing out coefficients of Z-series by Newton method so that its cumulants match with desired cumulants might not always work easily for many diffusions. I was thinking of an alternative easier method that can find the Z-series representation with given target cumulants when a "nearby density" whose Z-series representation is known is input to the numerical method. By "nearby density" I mean a density that is not very different in cumulants/moments from target density and shares the support of the target density. Sharing the support is important since we want to find a new density by multiplying the original density with a polynomial and if original density has smaller tails/support, there is no way to find a good density beyond those tails by any possible polynomial multiplication where the support of original input density is zero. But we describe how to increase the support of the original input density by linearly scaling it when target density has far larger variance so the method could still be used. 
We want to multiply the original "nearby density" with a polynomial whose coefficients are as yet unknown and then numerically calculate first few moments by integrating over the density. Variables get integrated and by equating each moment, we are left with a linear equation for each moment in terms of as yet unknown coefficients. Zeroth moment is the probability distribution itself whose integral is equated to unity to find an extra linear equation in terms of unknown coefficients. We can solve for example six linear equations by gaussian elimination to retrieve the as yet unknown coefficients of the polynomial that multiplies the original "nearby density" and results in a density with deaired cumulants. 
Though I explain things again, please read this post 1023 on page 69 of this thread where I have used this method earlier: Breakthrough in the theory of stochastic differential equations and their simulation - Page 69 - wilmott.com.

Suppose we have an original nearby density f(x) on a grid. We can easily calculate all moments/cumulants by numerically integrating over the density. We want to find a target density g(x) that has certain cumulants usually different from the cumulants of nearby density. We multiply the original density with a polynomial with as yet unknown coefficients to find the target density.

[$]g(x) \,=  \, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4] [$]

We do not know the coefficient c's and want to find them so that cumulants or equivalently the moments we have backed out from cumulants are matched with cumulants of g(x), the target density. 
We write as

Since resulting density g(x) has to be a valid probability distribution, our first condition is that integral of density should  to one.

[$]  \int  \, g(x)  \, dx =\int  \, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]   \, dx = \, 1  \, [$]   Equation(1)

if  [$]\mu_1[$] is the value of first raw moment backed out from target cumulants that we want to match, we must have 

[$]  \int  \,x \, g(x)  \, dx =\int  \, x\, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]   \, dx = \,\mu_1  \, [$]   Equation(2)

similarly if if  [$]\mu_2[$] is the value of second raw moment backed out from target cumulants, we must have 

[$]  \int  \,x^2 \, g(x)  \, dx =\int  \, x^2\, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]   \, dx = \,\mu_2  \, [$]   Equation(3)

if  [$]\mu_3[$] is the value of second raw moment backed out from target cumulants, we must have 

[$]  \int  \,x^3 \, g(x)  \, dx =\int  \, x^3\, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]   \, dx = \,\mu_3  \, [$]   Equation(3)

similarly we can write equations of desired raw moments backed out from central moments to fourth or sixth or eighth order.
From equation (1), we have 

[$]  \int  \, g(x)  \, dx =\int  \, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]   \, dx = \, 1  \, [$]   Equation(1)
[$] = c_0  \, \int  \, f(x)  \,  \, dx+ c_1  \int  \, f(x) \, x  \, dx + c_2  \int  \, f(x) \, x^2 \, dx + c_3 \int  \, f(x)  \, x^3 \, dx+ c_4  \int  \, f(x) \, x^4 \, dx  = \, 1  \, [$]

if we use the notation [$]E[x^n]=\int  \, f(x) \, x^n \, dx =  \zeta_n [$], we can write the moment equations as

[$]  c_0  \, + c_1 \, \zeta_1 + c_2 \, \zeta_2 + c_3 \, \zeta_3+ c_4 \, \zeta_4  = \, 1  \, [$]                            Equation(1)
[$]  c_0  \, \zeta_1 + c_1 \, \zeta_2 + c_2 \, \zeta_3 + c_3 \, \zeta_4+ c_4 \, \zeta_5  = \,\mu_1  \, \, [$]      Equation(2)
[$]  c_0  \, \zeta_2 + c_1 \, \zeta_3 + c_2 \, \zeta_4 + c_3 \, \zeta_5+ c_4 \, \zeta_6  = \,\mu_2  \, \, [$]      Equation(3)
[$]  c_0  \, \zeta_3 + c_1 \, \zeta_4 + c_2 \, \zeta_5 + c_3 \, \zeta_6+ c_4 \, \zeta_7  = \,\mu_3  \, \, [$]      Equation(4)

Since [$]E[x^n]=  \zeta_n [$]  are all known from numerical integration on original nearby density, we are left with a number of linear equations that we need to solve by Gaussian elimination to find out the coefficients c's of multiplying polynomial such that resulting density has cumulants equated with target cumulants. Once these equations are solved we find the resulting density as

[$]g(x) \,=  \, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4] [$]

There is one check that you will need. The density has to be positive everywhere. For well-stated moments with proper density support, this should not be a problem at all. But when bad cumulants are specified or support of the density does not match the support of target density, there could be a possibility that you would get a target density that is negative at some places as a result of solution of  above linear equations.

I will follow with one or two posts to complete description of the method. In new posts, I want to further explain how to use the above method when we have our original density f(x) as a Z_series and we want to find the Z-Series of target density g(x) such that cumulants of target density match the desired cumulants. I also want to mention simple scaling that can be used to increase the support of original density when its variance is smaller so its support matches the support of target density with larger variance. 
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: 1839
Joined: July 14th, 2002, 3:00 am

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

March 16th, 2022, 3:27 pm

As I stated in the previous post that our original intent is to take an initial "nearby density" that has a Z-series representation and convert it into a target density Z-series representation in a way that cumulants of target density match the desired cumulants. We named initial nearby density function as f(x) and target density function as g(x). Both pdfs of initial nearby density, f(x) and target density, g(x) are related according to formula

[$]g(x) \,=  \, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4] [$]

but two probability density functions are related by change of density derivative as
[$]g(x) \,=  \, f(x)  \, \frac{dx_f}{dx_g} [$]
which means that change of density derivative between two densities is given as
[$]\frac{dx_f}{dx_g} \,=   [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4] [$]
The above series is in x-coordinates and we want to convert it into Z-coordinates
Initially we know the Z-series representation of input nearby density as
[$] x_f(Z) \,= \, a_0 \,+ \,a_1 \,Z \,+\,a_2\, Z^2\,+\, a_3\, Z^3 \,+\,a_4 \,Z^4 [$]

We want to represent our series [$]\frac{dx_f}{dx_g} \,[$] which is a power series in x-coordinates to a power series in Z-coordinates using composition of series since we know Z-series representation of original density [$]x_f[$] which is given in equation above. Friends like me who are not experts in series can find about composition of series in mathematica or on the web (math.stackexchange must have examples).

Now the original nearby density with pdf f(x) is represented in normal coordinates as

[$]f(x(Z)) \,= \, p(Z) \, \frac{dZ}{dx_f} \,[$]

Similarly the target density in normal coordinates is represented as

[$]g(x(Z)) \,= \, p(Z) \, \frac{dZ}{dx_g} \,=\, p(Z) \, \frac{dZ}{dx_f} \,\frac{dx_f}{dx_g} \,[$]

which further means that 

[$]\frac{dx_g}{dZ} \,= \, \frac{dx_f}{dZ} \,\frac{dx_g}{dx_f} \,[$]

The above equation involves reciprocal of our series [$]\frac{dx_f}{dx_g} \,[$].

suppose after all the above calculations, the above term [$]\frac{dx_g}{dZ} \,[$] is a Z-series given as

[$]\frac{dx_g}{dZ} \,= b_1 \,+ \,b_2 \,Z \,+\,b_3\, Z^2\,+\, b_4\, Z^3 \,+\,b_5 \,Z^4 [$]
We can integrate the above series to find an expression of target density variable as a Z-series.

[$] x_g(Z) \,= IntegrationConstant + \,b_1 \,Z \,+\,b_2\, Z^2/2+\, b_3\, Z^3/3+\,b_4 \,Z^4/4 [$]

Integration constant in above series is median and can easily be found from equation of 1st cumulant/mean since target mean is known. If we absorb the numbers in the denominator of Z-series in its coefficients and represent  [$]x_g(Z)[$] as 

[$] x_g(Z) \,= \, b_0 \,+ \,b_1 \,Z \,+\,b_2\, Z^2\,+\, b_3\, Z^3 \,+\,b_4 \,Z^4 [$]
which means
[$]E[x_g]=\, b_0 \,+\,b_2\,+3\,b_4 = \mu_1  [$]
or
[$] b_0 = \mu_1 \,-\,b_2\,-3\,b_4  [$]

 all the terms in rhs of above equation are known which gives us median [$]b_0[$] and we know  Z-series of [$]x_g(Z)[$] since other coefficients are obtained from integration of [$]\frac{dx_g}{dZ} \,[$]

I will post a working program on this in a few days.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
 
User avatar
Amin
Topic Author
Posts: 1839
Joined: July 14th, 2002, 3:00 am

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

March 16th, 2022, 5:37 pm

Friends, I have some more ideas. Sorry for this abrupt interruption but please let me compose myself. I will come back in another day or two with better ideas. I think problem is even simpler than what I had previously considered but I want to check my ideas on computer for fear of making errors. I want to request friends to give me some time and I will soon come back with an easier solution.
My first post today is perfectly valid and useful but I am thinking of using the same ideas slightly differently.
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: 1839
Joined: July 14th, 2002, 3:00 am

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

March 18th, 2022, 12:12 pm

Friends I recently posted a cumulant matching program that advances the density of SDE diffusions (represented in the form of a Z-series) by adding the cumulants of existing density with cumulants of innovations(diffusion term) and then finding another SDE density in Z-series so that its cumulants match the added cumulants of innovations and existing density. This technique and its derivatives can be applied in interesting ways at a lot of places in applied probability theory and I am sure many friends would like to use it extensively. Apart from a lot of simple applied uses, it can be used to solve for example the densities where the driver noise in not normal but for example is a gamma density variable.   
However backing out coefficients of Z-series by Newton method so that its cumulants match with desired cumulants might not always work easily for many diffusions. I was thinking of an alternative easier method that can find the Z-series representation with given target cumulants when a "nearby density" whose Z-series representation is known is input to the numerical method. By "nearby density" I mean a density that is not very different in cumulants/moments from target density and shares the support of the target density. Sharing the support is important since we want to find a new density by multiplying the original density with a polynomial and if original density has smaller tails/support, there is no way to find a good density beyond those tails by any possible polynomial multiplication where the support of original input density is zero. But we describe how to increase the support of the original input density by linearly scaling it when target density has far larger variance so the method could still be used. 
We want to multiply the original "nearby density" with a polynomial whose coefficients are as yet unknown and then numerically calculate first few moments by integrating over the density. Variables get integrated and by equating each moment, we are left with a linear equation for each moment in terms of as yet unknown coefficients. Zeroth moment is the probability distribution itself whose integral is equated to unity to find an extra linear equation in terms of unknown coefficients. We can solve for example six linear equations by gaussian elimination to retrieve the as yet unknown coefficients of the polynomial that multiplies the original "nearby density" and results in a density with deaired cumulants. 
Though I explain things again, please read this post 1023 on page 69 of this thread where I have used this method earlier: Breakthrough in the theory of stochastic differential equations and their simulation - Page 69 - wilmott.com.

Suppose we have an original nearby density f(x) on a grid. We can easily calculate all moments/cumulants by numerically integrating over the density. We want to find a target density g(x) that has certain cumulants usually different from the cumulants of nearby density. We multiply the original density with a polynomial with as yet unknown coefficients to find the target density.

[$]g(x) \,=  \, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4] [$]

We do not know the coefficient c's and want to find them so that cumulants or equivalently the moments we have backed out from cumulants are matched with cumulants of g(x), the target density. 
We write as

Since resulting density g(x) has to be a valid probability distribution, our first condition is that integral of density should  to one.

[$]  \int  \, g(x)  \, dx =\int  \, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]   \, dx = \, 1  \, [$]   Equation(1)

if  [$]\mu_1[$] is the value of first raw moment backed out from target cumulants that we want to match, we must have 

[$]  \int  \,x \, g(x)  \, dx =\int  \, x\, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]   \, dx = \,\mu_1  \, [$]   Equation(2)

similarly if if  [$]\mu_2[$] is the value of second raw moment backed out from target cumulants, we must have 

[$]  \int  \,x^2 \, g(x)  \, dx =\int  \, x^2\, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]   \, dx = \,\mu_2  \, [$]   Equation(3)

if  [$]\mu_3[$] is the value of second raw moment backed out from target cumulants, we must have 

[$]  \int  \,x^3 \, g(x)  \, dx =\int  \, x^3\, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]   \, dx = \,\mu_3  \, [$]   Equation(3)

similarly we can write equations of desired raw moments backed out from central moments to fourth or sixth or eighth order.
From equation (1), we have 

[$]  \int  \, g(x)  \, dx =\int  \, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]   \, dx = \, 1  \, [$]   Equation(1)
[$] = c_0  \, \int  \, f(x)  \,  \, dx+ c_1  \int  \, f(x) \, x  \, dx + c_2  \int  \, f(x) \, x^2 \, dx + c_3 \int  \, f(x)  \, x^3 \, dx+ c_4  \int  \, f(x) \, x^4 \, dx  = \, 1  \, [$]

if we use the notation [$]E[x^n]=\int  \, f(x) \, x^n \, dx =  \zeta_n [$], we can write the moment equations as

[$]  c_0  \, + c_1 \, \zeta_1 + c_2 \, \zeta_2 + c_3 \, \zeta_3+ c_4 \, \zeta_4  = \, 1  \, [$]                            Equation(1)
[$]  c_0  \, \zeta_1 + c_1 \, \zeta_2 + c_2 \, \zeta_3 + c_3 \, \zeta_4+ c_4 \, \zeta_5  = \,\mu_1  \, \, [$]      Equation(2)
[$]  c_0  \, \zeta_2 + c_1 \, \zeta_3 + c_2 \, \zeta_4 + c_3 \, \zeta_5+ c_4 \, \zeta_6  = \,\mu_2  \, \, [$]      Equation(3)
[$]  c_0  \, \zeta_3 + c_1 \, \zeta_4 + c_2 \, \zeta_5 + c_3 \, \zeta_6+ c_4 \, \zeta_7  = \,\mu_3  \, \, [$]      Equation(4)

Since [$]E[x^n]=  \zeta_n [$]  are all known from numerical integration on original nearby density, we are left with a number of linear equations that we need to solve by Gaussian elimination to find out the coefficients c's of multiplying polynomial such that resulting density has cumulants equated with target cumulants. Once these equations are solved we find the resulting density as

[$]g(x) \,=  \, f(x)  \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4] [$]

There is one check that you will need. The density has to be positive everywhere. For well-stated moments with proper density support, this should not be a problem at all. But when bad cumulants are specified or support of the density does not match the support of target density, there could be a possibility that you would get a target density that is negative at some places as a result of solution of  above linear equations.

I will follow with one or two posts to complete description of the method. In new posts, I want to further explain how to use the above method when we have our original density f(x) as a Z_series and we want to find the Z-Series of target density g(x) such that cumulants of target density match the desired cumulants. I also want to mention simple scaling that can be used to increase the support of original density when its variance is smaller so its support matches the support of target density with larger variance. 
.
.
This was a good post. I would request friends to please disregard the post 1372 I made after this. But I thought more about the ideas and have something to share with friends.
In the copied post, we start with an initial density and we solve a set of linear equations to find a target density that matches target moments. I want to try to represent the target density variable as a transformation (given by power series) of the initial input density variable. Once we know this transformation, hopefully, we would even be able to evaluate expectations of functions of target density variable on initial density and  vice versa.

Suppose the initial density of variable X is given 
[$]p_X(x) =  f(x)  \, [$]

and we have the new target density 

[$]p_Y(x) =  f(x)   \, [c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4]=  f(x)   \, Q(x) [$] 
where
[$]Q(x)=[c_0  \, + c_1  \, x  + c_2  \, x^2 + c_3  \, x^3 + c_4  \, x^4][$]

We suppose that we can expand both f(x) and Q(x) as power series.
We want to find a transformation from X to Y,  [$]Y(X)[$] so that given any value(realization) of random variable X, the CDF of [$]Y(X)[$] on target density exactly matches the CDF of X on initial density. 
For that we write the CDF equations as

[$]  \int_{-\infty}^X  \, f(x) \, dx =  \, \int_{-\infty}^{Y(X)}  \, f(x)  \,Q(x) \, dx [$]  Equation(1)

X is natural variable for initial density so [$]\frac{dX}{dx}=1[$]

Differentiating our integral equations by x on both sides, we get

[$]    \, f(X) \, =  \,   \, f(Y(X))  \,Q(Y(X)) \, \frac{dY}{dX} [$]     Equation (A)

We had supposed that we could find power series of  f(X) and Q(X). Now we suppose that transformation relationship between X and Y is also given by a power series as
[$]Y(X)\,=\,a_0 \,+\, a_1  \,X \,+\, a_2\, X^2 \,+\, a_3\, X^3\, +\, a_4 \,X^4[$] Equation (B)

We notice that all terms in Equation(A) can be expanded as a power series. Equality between both sides of the equation would hold when power series expansions on both sides have equal coefficients on powers of X in which we are expanding the terms. We want to turn this equality of coefficients into equations(for each power of X) in terms of unknown coefficients [$]a_0\,, \,a_1 , \,a_2 , \,a_3 , \,a_4[$] . These equations would be algebraic expressions that can be solved possibly by root search. We can write rhs of equation (A) as an expansion as

[$]  \,   \, f(Y(X))  \,Q(Y(X)) \, \frac{dY}{dX} = \alpha_0(a_0,a_1,a_2,a_3,a_4) + \alpha_1(a_0,a_1,a_2,a_3,a_4) X + \alpha_2(a_0,a_1,a_2,a_3,a_4) X^2[$]
[$]+ \alpha_3(a_0,a_1,a_2,a_3,a_4) X^3+ \alpha_4(a_0,a_1,a_2,a_3,a_4) X^4 [$]  

while lhs of equation(A) can be expanded as

[$]  \,   \, f(X) = b_0 + \,b_1 X + \,b_2 X^2+ \,b_3 X^3+ \,b_4 X^4 [$]
which gives us five(or more) equations by matching coefficients of each power of X as  

[$]\alpha_0(a_0,a_1,a_2,a_3,a_4)=b_0 [$]
[$]\alpha_1(a_0,a_1,a_2,a_3,a_4)=b_1 [$]
[$]\alpha_2(a_0,a_1,a_2,a_3,a_4)=b_2 [$]
[$]\alpha_3(a_0,a_1,a_2,a_3,a_4)=b_3 [$]
[$]\alpha_4(a_0,a_1,a_2,a_3,a_4)=b_4 [$]
Each alpha above is an algebraic equation.
Once we solve the above equations by root search, we can find the transformation in the form of a power series that gives relationship between variable X of initial distribution with the variable Y of the target distribution as 

[$]Y(X)\,=\,a_0 \,+\, a_1  \,X \,+\, a_2\, X^2 \,+\, a_3\, X^3\, +\, a_4 \,X^4[$] Equation (B)

Once this relationship is known it could be possible to evaluate expectation (and its moments)  of variable Y(X) of the target density by integrating its series transformation relationship on the initial density as

[$]E[Y]=  \int_{-\infty}^{\infty} \,Y(X) \, f(X) \, dX =  \, \int_{-\infty}^{\infty} \, X\, f(X)  \,Q(X) \, dX[$]
or for its moments
[$]E[Y^n]=  \int_{-\infty}^{\infty} \,[Y(X)]^n \, f(X) \, dX =  \, \int_{-\infty}^{\infty} \, X^n\, f(X)  \,Q(X) \, dX[$]

So in the copied post we understood how to find a new target density whose cumulants match the target cumulants. And in this post, we have seen how to find a transformation between the two variables related to initial density and a different target density respectively.
I hope this post will be useful but please pardon any errors. 
In a few days I want to come back with a program that uses this method to simulate the densities of SDEs in time.
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: 1839
Joined: July 14th, 2002, 3:00 am

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

March 21st, 2022, 12:56 pm

Friends, I have remained busy with work on solution of densities of SDEs within Cumulant matching framework. I wrote a new program for cumulant matching. This program does not borrow equations from mathematica. To calculate cumulants and their first derivative with respect to series coefficients, I simply calculate the moments and their first derivative and convert them into cumulants and derivatives with respect to cumulants using moment to cumulant relationships. For moments, I have not expanded the series explicitly but used loops to multiply series repeatedly retaining all coefficients and then substituting the expected value of powers of normal random variable associated with each series term. I am going to post a full working but slightly preliminary program in half an hour. It is a working program but I want to add comments and streamline some more things so that you can specify arbitrary series order and related moment order that calculations have to take into account. All moment calculations are done algorithmically with loops (that are slow in matlab but faster in C++ etc.). You can even use the program for expansions with some other probability distribution different from normal and all you have to specify is to change the part of program that once calculates the expectations of powers of normal random random variable. Once you replace this small part of the program with your own modification to calculate powers of different random variable, you can easily calculate cumulants and their derivatives with a different base probability distribution other than normal.  In the older program, I had used two hermite polynomials for diffusion terms and each hermite polynomial was taking coefficients that was a series in a different independent normal random variable. Now you can take up to four hermite polynomials in cumulant calculations and you have to specify a series expansion in independent normal for each of the hermite polynomials. This was a bit involved but I algorithmically calculated the diffusion cumulants and you can now calculate cumulants of the diffusion term to any order after specifying coefficients(which are series themselves) of four hermite polynomials. 
The current program takes six cumulants and returns a series in a constant coefficient  and coefficients of five powers of the standard normal variable. Since everything is algorithmic, you can change the number of terms and cumulants used (but that will be assured in next version in a few days. You can try but there might be programming issues like different series order in different parts of the program and such small things that I will streamline in the next version.) You do not have to specify any median since constant coefficient of the series is calculated directly from cumulants (six cumulants and six coefficients upto fifth order after including a constant.)  
I will come back with very detailed explanation of the programs since otherwise it will be difficult for friends to understand the program. This program is very robust and better than all previous programs. Densities that go to zero are captured well  once they are initially going into zero but once a substantial portion of density is at zero, the program does not run properly (since we need to modify our calculations of cumulants with a delta mass at zero instead of totally continuous distributions something we have not done.)

A few graphs (For some reason my matlab was not changing the axes of graphs so I could not focus on the main body.). Five graphs are for densities of SDEs at .5 years and last five graphs are out to  five years.


Image

Image

Image

Image

Image

Image

Image

Image

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: 1839
Joined: July 14th, 2002, 3:00 am

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

March 21st, 2022, 1:30 pm

This is a full working but still preliminary program and I will update it with a better, neater and commented version streamlining some of the functionality that already exists. I will explain the program in total detail. Once friends understand the logic and nomenclature of of program, it is very easy and simple.
I am posting this program rather quickly since I am travelling tomorrow and then I may not have time to work properly for a few days. It is completed and working. All I need to do is to streamline a few things like parameters going into different functions and then you can possibly use it for different number of series terms and cumulants. I will try that friends could use it for up to eight cumulants. Some things have to be added for example turning moments to cumulants have equations to sixth order only and they are hardcopied as opposed to algorithmically generated.

[code]
function [] = FPESeriesCoeffsFromGaussiansWilmottNew02A()

%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 not directly simulated the SDE but simulated the transformed 
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.

dt=.125/2/2;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.
Tt=32*5;%16*2*4;%*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:64)=dt/16;
ds(65:128)=dt/4;
if(Tt<=4)
Ts=(Tt*16);%+(64-16);
ds(1:Ts)=dt/16;
end
if(Tt>4)&&(Tt<=20)
Ts=(64)+((Tt-4)*4);
ds(1:64)=dt/16;
ds(65:Ts)=dt/4;
end
if(Tt>20)
Ts=(64)+((20-4)*4)+(Tt-20);    
ds(1:64)=dt/16;
ds(65:128)=dt/4;
ds(129:Ts)=dt;
end
T
sum(ds(1:Ts))
Ts
%str=input('Check if time is right');

OrderA=4;  %
OrderM=4;  %

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



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

NnMid=4.0;

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


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

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ExpnOrder=5;
%SeriesOrder=4;

wnStart=1;
tic

for tt=1:Ts
    tt

    if(tt==1)
        %dt=ds(1);
        %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder);
 %       [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),ExpnOrder);
        
        [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
    %    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
        
      %  dwMu0dtdw
      %  dwMu0dtdwb
      %  wMu0dt
      %  wMu0dtb
      %  str=input('Look at drift');
        c0=wMu0dt+w0;
        c(1)=sigma0*sqrt(ds(tt));
        c(2:10)=0.0;
        w0=c0;
    
    elseif(tt<=6) 

    %dt=ds(tt);    
        
    %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
     %[b0,b] = CalculateDriftbCoeffs08O4(wMu0dt,dwMu0dtdw,c);
     [b0,b] = CalculateDriftbCoeffs08(wMu0dt,dwMu0dtdw,c,gamma);
     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);
     
     w0=c0;
            
    else
        
      
        
        
        
        
    %below wMu0dt is drift and dwMu0dtdw are its derivatives.
     [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
     [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
     [b0,b] = CalculateDriftbCoeffs08(wMu0dt,dwMu0dtdw,c,gamma);
     
     size(b)
     ba(1:6)=b(1:6);
     
     
    [wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),8);
    
    [wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),8);
    
    [wVol3dt,dwVol3dtdw] = BesselVolAndDerivativesH3(w0,YqCoeff0,Fp1,gamma,ds(tt),8);

    [g10,g1] = CalculateDriftbCoeffs08O4(wVol0dt,dwVol0dtdw,c);
    
    g1(5:6)=0;
    
    
    [g20,g2] = CalculateDriftbCoeffs08O4(wVol2dt,dwVol2dtdw,c);
 g2(5:6)=0;
 
    [g30,g3] = CalculateDriftbCoeffs08O4(wVol3dt,dwVol3dtdw,c);
 
 g3(5:6)=0;
 
 g40=0;
 g4(1:6)=0;

 a0=c0+b0;
 
 a(1:5)=c(1:5)+b(1:5);
    
     [a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A(a0,a,g10,g1,g20,g2,g30,g3,g40,g4,sigma0,ds(tt));
     %[a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A(a0,a,b0,ba,g0,g,h0,h,sigma0,ds(tt));
     
     
 c0=a0;
 c(1)=a1;
 c(2)=a2;
 c(3)=a3;
 c(4)=a4;
 c(5)=a5;
    
    w0=c0;
    
    end

end
wnStart=1;


yy0=((1-gamma).*w0).^(1/(1-gamma));

 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=ceil(Nn/2)-1:-1:1
     if((w(nn)<0)&&(Flag==0))
         wnStart=nn+1;
         Flag=1;
     end
 end
 
 
 %c0
 %c
 
 
 %str=input('Look at numbers')
 yy0(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));

 yy=yy0;
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

ItoHermiteMean=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %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


theta1=1;
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=30;
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(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');
 
title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
 
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
 

end[/code]

.
.
2nd function.

[code]%function [a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A(a0,a,c0,c,g0,g,h0,h,sigma0,dt)
function [a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A(a0,a,g10,g1,g20,g2,g30,g3,g40,g4,sigma0,dt)


a1=a(1);
a2=a(2);
a3=a(3);
a4=a(4);
a5=a(5);
[F] = CalculateSixCumulantsParamsFive(a0,a1,a2,a3,a4,a5);

C1=F(1,1);
C2=F(2,1);
C3=F(3,1);
C4=F(4,1);
C5=F(5,1);
C6=F(6,1);

g10=g10+sigma0*sqrt(dt);
%[dC2,dC3,dC4] = CalculateCumulantsHighOrder(g0,g,h0,h);

 SeriesOrder=6;
      ab00=0;
      ab0(1:SeriesOrder)=0;
%      b30=0;
%      b3(1:SeriesOrder)=0;
%      b40=0;
%      b4(1:SeriesOrder)=0;
 [dC1,dC2,dC3,dC4,dC5,dC6] = CalculateDiffusionCumulants(ab00,ab0,g10,g1,g20,g2,g30,g3,g40,g4,4,6);
 %[dC1,dC2,dC3,dC4,dC5,dC6] = CalculateDiffusionCumulants(ab00,ab0,g0,g,h0,h,b30,b3,b40,b4,4,6);


%dC1=0;
C1=C1+dC1;
C2=C2+dC2;
C3=C3+dC3;
C4=C4+dC4;
C5=C5+dC5;
C6=C6+dC6;



a1=a(1);
a2=a(2);
a3=a(3);
a4=a(4);
a5=a(5);

DZa0=a(1);
DZa(1)=2*a(2);
DZa(2)=3*a(3);
DZa(3)=4*a(4);
DZa(4)=5*a(5);

[DZaI0,DZaI] = SeriesReciprocal(DZa0,DZa,4);


a1=a1+.5*g10^2*DZaI0;
a2=a2+.5*g10^2*DZaI(1)/1;
a3=a3+.5*g10^2*DZaI(2)/1;
a4=a4+.5*g10^2*DZaI(3)/1;
a5=a5+.5*g10^2*DZaI(4)/1;

% a0
% a1
% a2
% a3
% a4
% a5
%str=input('Look at coefficients');
da(1,1)=a0;
da(2,1)=a1;
da(3,1)=a2;
da(4,1)=a3;
da(5,1)=a4;
da(6,1)=a5;

a(1)=a1;
a(2)=a2;
a(3)=a3;
a(4)=a4;
a(5)=a5;

%a0=0;
%[Fa,dFa] = CalculateSixCumulantsParamsFiveAndDerivs(C1,C2,C3,C4,C5,C6,a0,a1,a2,a3,a4,a5);
%Fa
%[F,dF] = CalculateCumulantsAndDerivativesFromMoments(C1,C2,C3,C4,C5,C6,a0,a,5,6,6);
[F,dF] = CalculateCumulantsAndDerivativesFromMoments(C1,C2,C3,C4,C5,C6,a0,a1,a2,a3,a4,a5,5,6,6);
%F-Fa
%dF-dFa
%str=input('Look at coefficients');
%[F,dF] = CalculateCumulantsAndDerivativesFromMoments(C1,C2,C3,C4,C5,C6,a0,a,SeriesOrder,NZterms,NMoments)


ObjBest=F(1,1)^2+F(2,1)^2+F(3,1)^2+F(4,1)^2+F(5,1)^2+F(6,1)^2;

a0Best=a0;
a1Best=a1;
a2Best=a2;
a3Best=a3;
a4Best=a4;
a5Best=a5;

nn=0;
while((nn<5)&&((abs(F(1,1))>.000000000001) || (abs(F(2,1))>.000000000001) || (abs(F(3,1))>.000000000001) || (abs(F(4,1))>.00000000001)|| (abs(F(5,1))>.0000000001)|| (abs(F(6,1))>.000000001)))
nn=nn+1;
da=da-dF\F;

%inv(dF)

a0=da(1,1);
a1=da(2,1);
a2=da(3,1);
a3=da(4,1);
a4=da(5,1);
a5=da(6,1);

%[F,dF] = CalculateCumulantsAndDerivatives(dC1,dC2,dC3,dC4,a1,a2,a3,a4,da0,da1,da2,da3,da4);
%[F,dF] = CalculateSixCumulantsParamsFiveAndDerivs(C1,C2,C3,C4,C5,C6,a0,a1,a2,a3,a4,a5);
%function [F,dF] = CalculateCumulantsAndDerivativesFromMoments(C1,C2,C3,C4,C5,C6,a0,a1,a2,a3,a4,a5,SeriesOrder,NZterms,NMoments)
[F,dF] = CalculateCumulantsAndDerivativesFromMoments(C1,C2,C3,C4,C5,C6,a0,a1,a2,a3,a4,a5,5,6,6);

 ObjNew=F(1,1)^2+F(2,1)^2+F(3,1)^2+F(4,1)^2+F(5,1)^2+F(6,1)^2;
 
  if(ObjBest<ObjNew)
      a0=a0Best;
      a1=a1Best;
      a2=a2Best;
      a3=a3Best;
      a4=a4Best;
      a5=a5Best;
  else
      ObjBest=ObjNew;
      a0Best=a0;
      a1Best=a1;
      a2Best=a2;
      a3Best=a3;
      a4Best=a4;
      a5Best=a5;
  end

da(1,1)=a0;
da(2,1)=a1;
da(3,1)=a2;
da(4,1)=a3;
da(5,1)=a4;
da(6,1)=a5;



end
  a0=a0Best;   
  a1=a1Best;
  a2=a2Best;
  a3=a3Best;
  a4=a4Best;
  a5=a5Best;

  
% a0
% a1
% a2
% a3
% a4
nn
%str=input('Look at results');

end

[/code]
.
.
.
[code]function [F,dF] = CalculateCumulantsAndDerivativesFromMoments(C1,C2,C3,C4,C5,C6,a0,a1,a2,a3,a4,a5,SeriesOrder,NZterms,NMoments)


a(1)=a1;
a(2)=a2;
a(3)=a3;
a(4)=a4;
a(5)=a5;


aa0=a0;
a0=0;

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:50)=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:50)=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);
u5=EXZ(6,1);
u6=EXZ(7,1);


% u1
% u2
% u3
% u4
% u5
% u6

%str=input('Look at numbers')

du1(1)=0;
du2(1)=0;
du3(1)=0;
du4(1)=0;
du5(1)=0;
du6(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);
du5(mm)=5*EXZ(5,mm);
du6(mm)=6*EXZ(6,mm);
end

%du1(1)=1;



k1=aa0+a(2)+3*a(4);
%k1=u1;
k2=u2-u1^2;
k3=u3-3*u2*u1+2*u1^3;
k4=u4-4*u3*u1-3*u2^2+12*u2*u1^2-6*u1^4;
k5=u5-5*u4*u1-10*u3*u2+20*u3*u1^2+30*u2^2*u1-60*u2*u1^3+24*u1^5;
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;


% k1
% k2
% k3
% k4
% k5
% k6

% str=input('Look at numbers--2')
dk1(1)=1.0;
dk1(2)=0.0;
dk1(3)=1.0;
dk1(4)=0.0;
dk1(5)=3.0;
dk1(6)=0.0;

dk2(1)=0;
dk3(1)=0;
dk4(1)=0;
dk5(1)=0;
dk6(1)=0;

for mm=2:SeriesOrder+1
%dk1(mm)=du1(mm);
dk2(mm)=-2.* u1.*du1(mm)+du2(mm);
dk3(mm)=du3(mm)-3*du1(mm)*u2-3*u1*du2(mm)+2*3*u1^2*du1(mm);
dk4(mm)=du4(mm)-4*du3(mm)*u1-4*u3*du1(mm)-3*2*u2*du2(mm)+12*du2(mm)*u1^2+12*2*u2*u1*du1(mm)-6*4*u1^3*du1(mm);
dk5(mm)=du5(mm)-5*du4(mm)*u1-5*u4*du1(mm)-10*du3(mm)*u2-10*u3*du2(mm)+20*du3(mm)*u1^2+20*2*u3*u1*du1(mm)+ ...
    30*2*u2*du2(mm)*u1+30*u2^2*du1(mm)-60*du2(mm)*u1^3-60*3*u2*u1^2*du1(mm)+24*5*u1^4*du1(mm);
dk6(mm)=du6(mm)-6*du5(mm)*u1-6*u5*du1(mm)-15*du4(mm)*u2-15*u4*du2(mm)+30*du4(mm)*u1^2+30*2*u4*u1*du1(mm)-10*2*u3*du3(mm)+ ...
    120*du3(mm)*u2*u1-120*u3*du2(mm)*u1-120*u3*u2*du1(mm)-120*du3(mm)*u1^3-120*3*u3*u1^2*du1(mm)+30*3*u2^2*du2(mm) - ...
    270*2*u2*du2(mm)*u1^2+270*2*u2^2*u1*du1(mm)+360*du2(mm)*u1^4+360*4*u2*u1^3*du1(mm)-120*6*u1^5*du1(mm);

end





F(1,1)=k1-C1;
 F(2,1)=k2-C2;
 F(3,1)=k3-C3;
 F(4,1)=k4-C4;
 F(5,1)=k5-C5;
 F(6,1)=k6-C6;

 
     for mm=1:SeriesOrder+1
        dF(1,mm)=dk1(mm);
        dF(2,mm)=dk2(mm);
        dF(3,mm)=dk3(mm);
        dF(4,mm)=dk4(mm);
        dF(5,mm)=dk5(mm);
        dF(6,mm)=dk6(mm);
     end

end

[/code]
.
.
.
[code]function [dC1,dC2,dC3,dC4,dC5,dC6] = CalculateDiffusionCumulants(d00,d0,d10,d1,d20,d2,d30,d3,d40,d4,SeriesOrder,NMoments)



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;

%[dC2,dC3,dC4] = CalculateCumulantsHighOrder(g0,g,h0,h)

a0=d0-d2+3*d4;
a00=d00-d20+3*d40;
a0(SeriesOrder+1:SeriesOrder*NMoments)=0;

%size(d1)
%size(d3)

a1=d1-3*d3;
a10=d10-3*d30;
a1(SeriesOrder+1:SeriesOrder*NMoments)=0;

a2=d2-6.*d4;
a20=d20-6.*d40;
a2(SeriesOrder+1:SeriesOrder*NMoments)=0;

a3=d3;
a30=d30;
a3(SeriesOrder+1:SeriesOrder*NMoments)=0;

a4=d4;
a40=d40;
a4(SeriesOrder+1:SeriesOrder*NMoments)=0;







%                     str=input('Look at indices');
T0(1:NMoments,1:NMoments*SeriesOrder)=0;                    
T0(1,1:SeriesOrder)=a0(1:SeriesOrder);
T00(1)=a00;

T1(1:NMoments,1:NMoments*SeriesOrder)=0;                    
T1(1,1:SeriesOrder)=a1(1:SeriesOrder);
T10(1)=a10;

T2(1:NMoments,1:NMoments*SeriesOrder)=0; 
T2(1,1:SeriesOrder)=a2(1:SeriesOrder);
T20(1)=a20;

T3(1:NMoments,1:NMoments*SeriesOrder)=0; 
T3(1,1:SeriesOrder)=a3(1:SeriesOrder);
T30(1)=a30;

T4(1:NMoments,1:NMoments*SeriesOrder)=0; 
T4(1,1:SeriesOrder)=a4(1:SeriesOrder);
T40(1)=a40;
for k = 2 : (NMoments)                    
b=T0(k-1,1:SeriesOrder*k);
b0=T00(k-1);
[b0,b] =SeriesProduct(a00,a0,b0,b,SeriesOrder*k);                    
T0(k,1:SeriesOrder*k)=b(1:SeriesOrder*k);
T00(k)=b0;  
    
    
b=T1(k-1,1:SeriesOrder*k);
b0=T10(k-1);
[b0,b] =SeriesProduct(a10,a1,b0,b,SeriesOrder*k);                    
T1(k,1:SeriesOrder*k)=b(1:SeriesOrder*k);
T10(k)=b0;  

b=T2(k-1,1:SeriesOrder*k);
b0=T20(k-1);
[b0,b] =SeriesProduct(a20,a2,b0,b,SeriesOrder*k);                    
T2(k,1:SeriesOrder*k)=b(1:SeriesOrder*k);  
T20(k)=b0;

b=T3(k-1,1:SeriesOrder*k);
b0=T30(k-1);
[b0,b] =SeriesProduct(a30,a3,b0,b,SeriesOrder*k);                    
T3(k,1:SeriesOrder*k)=b(1:SeriesOrder*k);  
T30(k)=b0;

b=T4(k-1,1:SeriesOrder*k);
b0=T40(k-1);
[b0,b] =SeriesProduct(a40,a4,b0,b,SeriesOrder*k);                    
T4(k,1:SeriesOrder*k)=b(1:SeriesOrder*k);  
T40(k)=b0;

end

CC(1,1,1,1,1)=1;


Moment1(1:NMoments)=0;
for k = 0 : (NMoments)
    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;
                for p = 0:j
                    l1=j-p+1;
                    l_0=p+1;
                    if(k>0)
                        b(1:SeriesOrder*k)=0;
                        b0=1;
                        %CC(l_0,l1,l2,l3,l4)=1.0;
                        if(l_0>1)
                            a=T0(l_0-1,1:SeriesOrder*k);
                            a0=T00(l_0-1);
                            [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*k);
                            %CC(l_0,l1,l2,l3,l4)=CC(l_0,l1,l2,l3,l4).*4^(l_0-1);
                        end  
                        if(l1>1)
                            a=T1(l1-1,1:SeriesOrder*k);
                            a0=T10(l1-1);
                            [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*k);
                            %CC(l_0,l1,l2,l3,l4)=CC(l_0,l1,l2,l3,l4).*2^(l1-1);
                        end    
                        if(l2>1)
                            a=T2(l2-1,1:SeriesOrder*k);
                            a0=T20(l2-1);
                            [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*k);
                            %CC(l_0,l1,l2,l3,l4)=CC(l_0,l1,l2,l3,l4).*2^(l2-1);
                        end
                        if(l3>1)
                            a=T3(l3-1,1:SeriesOrder*k);
                            a0=T30(l3-1);
                            [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*k);
                            %CC(l_0,l1,l2,l3,l4)=CC(l_0,l1,l2,l3,l4).*0^(l3-1);
                        end
                        if(l4>1)
                            a=T4(l4-1,1:SeriesOrder*k);
                            a0=T40(l4-1);
                            [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*k);
                            %CC(l_0,l1,l2,l3,l4)=CC(l_0,l1,l2,l3,l4).*0^(l4-1);
                        end

                      %  CC(l_0,l1,l2,l3,l4)=CC(l_0,l1,l2,l3,l4).*4^(l_0-1).*2^(l1-1)*2^(l2-1)*0^(l3-1)*0^(l4-1);
                        CC(l_0,l1,l2,l3,l4)=b0;
                        for qq=1:SeriesOrder*k
                            CC(l_0,l1,l2,l3,l4)=CC(l_0,l1,l2,l3,l4)+b(qq)*EZ(qq);
                        end
                        
                        
%                     Moment1(k)=Moment1(k)+factorial(k)/(factorial(l_0-1).*factorial(l1-1).*factorial(l2-1).*factorial(l3-1).*factorial(l4-1)).*c0(l_0).*CC(l1,l2,l3,l4);%.*EZ((l1-1)+2*(l2-1)+3*(l3-1)+4*(l4-1));
%                     k
%                     l_0
%                     l1
%                     l2
%                     l3
%                     l4
%                     
%                     str=input('Look at indices');
                    end
                end
            end
        end
    end
end



Moment1(1:NMoments)=0;
for k = 0 : (NMoments)
    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;
                for p = 0:j
                    l1=j-p+1;
                    l_0=p+1;
                    if((k>0)&&((l_0+l1+l2+l3+l4)>5))
                        if((l1+l2+l3+l4==4))
                             Moment1(k)=Moment1(k)+factorial(k)/(factorial(l_0-1).*factorial(l1-1).*factorial(l2-1).*factorial(l3-1).*factorial(l4-1)).*CC(l_0,l1,l2,l3,l4);%.*EZ((l1-1)+2*(l2-1)+0*3*(l3-1)+0*4*(l4-1));
                        else
                             Moment1(k)=Moment1(k)+factorial(k)/(factorial(l_0-1).*factorial(l1-1).*factorial(l2-1).*factorial(l3-1).*factorial(l4-1)).*CC(l_0,l1,l2,l3,l4).*EZ((l1-1)+2*(l2-1)+3*(l3-1)+4*(l4-1));
                        end
                        %                     k
                    % l_0
                    % l1
                    % l2
                    % l3
                    % l4
                    %Moment1(k)=Moment1(k)+factorial(k)/(factorial(l_0-1).*factorial(l1-1).*factorial(l2-1).*factorial(l3-1).*factorial(l4-1)).*CC(l_0,l1,l2,l3,l4);%.*EZ((l1-1)+2*(l2-1)+0*3*(l3-1)+0*4*(l4-1));
%                     k
%                     l_0
%                     l1
%                     l2
%                     l3
%                     l4
%                     
%                     str=input('Look at indices');
                    end
                end
            end
        end
    end
end

%Moment1
%str=input('Look at Moments');

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



dC1=u1;
dC2=u2-u1^2;
dC3=u3-3*u1*u2+2*u1^3;
dC4=u4-4*u3*u1-3*u2^2+12*u2*u1^2-6*u1^4;
dC5=u5-5*u4*u1-10*u3*u2+20*u3*u1^2+30*u2^2*u1-60*u2*u1^3+24*u1^5;
dC6=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;
%dC5=0;
%dC6=0;




end

[/code]
.
.
.
[code]function [F] = CalculateSixCumulantsParamsFive(a0,a1,a2,a3,a4,a5)
%[F,dF] = CalculateCumulantsAndDerivatives(dC1,dC2,dC3,dC4,a1,a2,a3,a4,da0,da1,da2,da3,da4);

k1=a0+a2+3*a4;

k2=a1^2+2*a2^2+6*a1*a3+15*a3^2+24*a2*a4+96*a4^2+ ... 
30*(a1+7*a3)*a5+945*a5^2;
 
k3=8*a2^3+216*a2^2*a4+6*a1^2*(a2+6*a4)+72*a1*a3*(a2+8*a4)+ ...
540*a1*(a2+10*a4)*a5+ ...
18*a2*(15*a3^2+128*a4^2+280*a3*a5+1575*a5^2)+ ...
108*a4*(25*a3^2+88*a4^2+560*a3*a5+3675*a5^2);
 
k4=24*(2*a2^4+96*a2^3*a4+a1^3*(a3+10*a5)+ ...
30*a2^2*(6*a3^2+64*a4^2+140*a3*a5+945*a5^2)+ ...
a1^2*(2*(a2^2+9*a3^2+18*a2*a4+96*a4^2)+375*a3*a5+ ...
2250*a5^2)+ ...
36*a2*a4*(125*a3^2+528*a4^2+3360*a3*a5+25725*a5^2)+ ...
9*(45*a3^4+3560*a3^2*a4^2+8704*a4^4+ ...
35*a3*(69*a3^2+3104*a4^2)*a5+ ...
700*(79*a3^2+1332*a4^2)*a5^2+632625*a3*a5^3+ ...
3018750*a5^4)+ ...
3*a1*(12*a2^2*(a3+10*a5)+8*a2*a4*(32*a3+375*a5)+ ...
3*(15*a3^3+528*a3*a4^2+530*a3^2*a5+7120*a4^2*a5+ ...
7175*a3*a5^2+36750*a5^3)));

k5=24*(16*a2^5+5*a1^4*a4+1200*a2^4*a4+ ...
10*a1^3*(3*a2*a3+42*a3*a4+40*a2*a5+570*a4*a5)+ ...
300*a2^3*(10*a3^2+128*a4^2+280*a3*a5+2205*a5^2)+ ...
1080*a2^2*a4*(125*a3^2+616*a4^2+3920*a3*a5+34300*a5^2)+ ...
30*a1*(a3*(16*a2^3+225*a2*a3^2+640*a2^2*a4+3330*a3^2*a4+ ...
9504*a2*a4^2+52224*a4^3)+ ...
10*(20*a2^3+954*a2*a3^2+900*a2^2*a4+15735*a3^2*a4+ ...
14952*a2*a4^2+91296*a4^3)*a5+ ...
3675*a3*(41*a2+750*a4)*a5^2+15750*(56*a2+1129*a4)*a5^3)+ ...
27*a4*(14775*a3^4+507200*a3^2*a4^2+950784*a4^4+ ...
700*a3*(1469*a3^2+27904*a4^2)*a5+ ...
1750*(16973*a3^2+118128*a4^2)*a5^2+419002500*a3*a5^3+ ...
2419134375*a5^4)+ ...
10*a1^2*(2*a2^3+72*a2^2*a4+ ...
3*a2*(24*a3^2+320*a4^2+625*a3*a5+4500*a5^2)+ ...
9*a4*(109*a3^2+528*a4^2+3130*a3*a5+24975*a5^2))+ ...
90*a2*(270*a3^4+69632*a4^4+16905*a3^3*a5+8391600*a4^2*a5^2+ ...
30187500*a5^4+280*a3^2*(89*a4^2+1580*a5^2)+ ...
35*a3*(24832*a4^2*a5+162675*a5^3)));

k6=120*(32*a2^6+3456*a2^5*a4+6*a1^5*a5+ ...
3*a1^4*(9*a3^2+8*a4*(2*a2+21*a4)+330*a3*a5+2850*a5^2)+ ...
720*a2^4*(15*a3^2+224*a4^2+490*a3*a5+4410*a5^2)+ ...
6048*a2^3*a4*(125*a3^2+704*a4^2+4480*a3*a5+44100*a5^2)+ ...
12*a1^3*(12*a3*(a2^2+6*a3^2+35*a2*a4+298*a4^2)+ ...
10*(20*a2^2+351*a3^2+684*a2*a4+6132*a4^2)*a5+ ...
60450*a3*a5^2+374625*a5^3)+ ...
216*a2^2*(945*a3^4+99680*a3^2*a4^2+313344*a4^4+ ...
420*a3*(161*a3^2+9312*a4^2)*a5+ ...
25200*(79*a3^2+1665*a4^2)*a5^2+28468125*a3*a5^3+ ...
166031250*a5^4)+ ...
1296*a2*a4*(6*(985*a3^4+38040*a3^2*a4^2+79232*a4^4)+ ...
35*a3*(13221*a3^2+279040*a4^2)*a5+ ...
1925*(7715*a3^2+59064*a4^2)*a5^2+230451375*a3*a5^3+ ...
1451480625*a5^4)+ ...
18*a1*(4*a3*(20*a2^4+1280*a2^3*a4+ ...
27*a2^2*(25*a3^2+1232*a4^2)+ ...
6*a2*(3885*a3^2*a4+69632*a4^3)+ ...
9*(135*a3^4+24280*a3^2*a4^2+237696*a4^4))+ ...
15*(80*a2^4+5600*a2^3*a4+56*a2^2*(159*a3^2+2848*a4^2)+ ...
64*a2*(5245*a3^2*a4+34236*a4^3)+ ...
3*(9717*a3^4+1145552*a3^2*a4^2+4076032*a4^4))*a5+ ...
2100*a3*(1148*a2^2+8199*a3^2+47250*a2*a4+ ...
526740*a4^2)*a5^2+ ...
126000*(126*a2^2+2935*a3^2+5645*a2*a4+68241*a4^2)*a5^3+ ...
4314397500*a3*a5^4+21772209375*a5^5)+ ...
27*(9945*a3^6+2894760*a3^4*a4^2+61140480*a3^2*a4^4+ ...
92930048*a4^6+ ...
210*a3*(5553*a3^4+1177072*a3^2*a4^2+13526272*a4^4)*a5+ ...
2100*(29649*a3^4+4101044*a3^2*a4^2+17011776*a4^4)*a5^2+ ...
3500*a3*(551015*a3^2+41249952*a4^2)*a5^3+ ...
39375*(923683*a3^2+24833432*a4^2)*a5^4+ ...
394058306250*a3*a5^5+1907724656250*a5^6)+ ...
6*a1^2*(8*a2^4+480*a2^3*a4+ ...
180*a2^2*(4*a3^2+64*a4^2+125*a3*a5+1050*a5^2)+ ...
72*a2*a4*(327*a3^2+1848*a4^2+10955*a3*a5+99900*a5^2)+ ...
9*(225*a3^4+22824*a3^2*a4^2+69632*a4^4+15090*a3^3*a5+ ...
830240*a3*a4^2*a5+412925*a3^2*a5^2+8239800*a4^2*a5^2+ ...
5475750*a3*a5^3+29636250*a5^4)));




 F(1,1)=k1;
 F(2,1)=k2;
 F(3,1)=k3;
 F(4,1)=k4;
 F(5,1)=k5;
 F(6,1)=k6;
       
end

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


NoOfTerms=19;

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

[/code]
.
.
.
[code]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[/code]
.
.
.
[code]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[/code]
.
.
.
[code]function [wVol3dt,dwVol3dtdw] = BesselVolAndDerivativesH3(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder)


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

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



Fp2(1)=Fp1(1,1,2,4);
Fp2(2)=Fp1(1,2,1,4);
Fp2(3)=Fp1(2,1,1,4);


wVol3dt0=0;
dwVol3dt(1:ExpnOrder)=0.0;

wVol3dt=0;
dwVol3dtdw(1:ExpnOrder)=0.0;

for mm=1:NoOfTerms

    wVol3dt0=YqCoeffa(mm).*w0.^Fp2(mm);
    for nn=1:ExpnOrder
        if(nn==1)
            dwVol3dt(nn)=wVol3dt0*Fp2(mm)*1/w0;
        else
        dwVol3dt(nn)=dwVol3dt(nn-1)*(Fp2(mm)-(nn-1))/w0;
        end
    end
    wVol3dt=wVol3dt+wVol3dt0;
    for nn=1:ExpnOrder
        dwVol3dtdw(nn)=dwVol3dtdw(nn)+dwVol3dt(nn);
    end
        
    
end



end[/code]
.
.
.
[code]function [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,ExpnOrder)



NoOfTerms=9;%excluding dt^3 terms

YqCoeff(1:NoOfTerms)=0.0;
Fp(1:NoOfTerms)=0.0;
a0=((beta1-gamma)/(1-gamma));
b0=((beta2-gamma)/(1-gamma));
A=mu1;
B=mu2;
C=-.5*sigma0^2*gamma;
dt2=dt^2/2;

YqCoeff(1:9)=0;

YqCoeff(1)=A*dt;
YqCoeff(2)=B*dt;
YqCoeff(3)=C*dt;

 YqCoeff(4)=A^2*a0*dt2*(1-gamma);
 YqCoeff(5)=B^2*b0*dt2*(1-gamma);
 YqCoeff(6)=1*1.5*.66*C^2*(-1)*dt2*(1-gamma)+1*1.5*.66* C*(-1)*(-2)*.5*sigma0^2 *dt2*(1-gamma)^2;
 YqCoeff(7)=A*B*(a0+b0)*dt2*(1-gamma);
 YqCoeff(8)=1*1.5*.66*A*C*(a0-1)*dt2*(1-gamma)^1+1*1.5*.66*.5*sigma0^2*A*a0*(a0-1)*dt2*(1-gamma)^2;
 YqCoeff(9)=1*1.5*.66*B*C*(b0-1)*dt2*(1-gamma)^1+1*1.5*.66*.5*sigma0^2*B*b0*(b0-1)*dt2*(1-gamma)^2;
 



Fp(1)=(beta1-gamma)/(1-gamma);
Fp(2)=(beta2-gamma)/(1-gamma);
Fp(3)=-1;

Fp(4)=2*(beta1-gamma)/(1-gamma)-1;
Fp(5)=2*(beta2-gamma)/(1-gamma)-1;
Fp(6)=-3;
Fp(7)=(beta1-gamma)/(1-gamma)+(beta2-gamma)/(1-gamma)-1;
Fp(8)=(beta1-gamma)/(1-gamma)-2;
Fp(9)=(beta2-gamma)/(1-gamma)-2;


 Ratio=2/dt*(YqCoeff(4)*((1-gamma)*w0).^Fp(4)+YqCoeff(5)*((1-gamma)*w0).^Fp(5)+YqCoeff(6)*((1-gamma)*w0).^Fp(6)+ ...
     YqCoeff(7)*((1-gamma)*w0).^Fp(7)+YqCoeff(8)*((1-gamma)*w0).^Fp(8)+YqCoeff(9)*((1-gamma)*w0).^Fp(9))/ ...
     (YqCoeff(1)*((1-gamma)*w0).^Fp(1)+YqCoeff(2)*((1-gamma)*w0).^Fp(2)+YqCoeff(3)*((1-gamma)*w0).^Fp(3));

Fp2(1:9)=Fp(1:9);%/(1-gamma);

wMu0dt0=0;
dwMu0dt(1:ExpnOrder)=0.0;

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

for mm=1:NoOfTerms

    wMu0dt0=YqCoeff(mm).*((1-gamma)*w0).^Fp2(mm);
    for nn=1:ExpnOrder
        if(nn==1)
            dwMu0dt(nn)=wMu0dt0*Fp2(mm)*1/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[/code]
.
.
.
[code]function [b0,b] = CalculateDriftbCoeffs08(wMu0dt,dwMu0dtdw,a,gamma)



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


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


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




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

b(7:10)=0.0; 

Bound=abs(b(1)/b0)*gamma/(1-gamma);
%Bound=1/abs(b0/b(1))*gamma;%/(1-gamma);

for nn=1:9
   
    if(abs(b(nn+1))>abs(b(nn)*Bound))
  %      b(nn+1)=sign(b(nn+1))*abs(b(nn))*Bound;    
    end
end

%[b0,b] = ApplyBoundsOnSeriesCoeffs(b0,b);


 
 % 


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




% 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

[/code]
.
.
.
[code]function [b0,b] = CalculateDriftbCoeffs08O4(wMu0dt,dwMu0dtdw,a)



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



end

[/code]
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: 1839
Joined: July 14th, 2002, 3:00 am

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

March 21st, 2022, 1:45 pm

Sorry friends that I could not post the code properly in code boxes. I am using a new laptop with Edge browser and that might be the problem. I will try re-posting the matlab program later. 
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: 1839
Joined: July 14th, 2002, 3:00 am

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

March 21st, 2022, 2:01 pm

function [] = FPESeriesCoeffsFromGaussiansWilmottNew02A()

%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 not directly simulated the SDE but simulated the transformed 
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.

dt=.125/2/2;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.
Tt=32*5;%16*2*4;%*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:64)=dt/16;
ds(65:128)=dt/4;
if(Tt<=4)
Ts=(Tt*16);%+(64-16);
ds(1:Ts)=dt/16;
end
if(Tt>4)&&(Tt<=20)
Ts=(64)+((Tt-4)*4);
ds(1:64)=dt/16;
ds(65:Ts)=dt/4;
end
if(Tt>20)
Ts=(64)+((20-4)*4)+(Tt-20);    
ds(1:64)=dt/16;
ds(65:128)=dt/4;
ds(129:Ts)=dt;
end
T
sum(ds(1:Ts))
Ts
%str=input('Check if time is right');

OrderA=4;  %
OrderM=4;  %

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



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

NnMid=4.0;

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


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

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ExpnOrder=5;
%SeriesOrder=4;

wnStart=1;
tic

for tt=1:Ts
    tt

    if(tt==1)
        %dt=ds(1);
        %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder);
 %       [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),ExpnOrder);
        
        [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
    %    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
        
      %  dwMu0dtdw
      %  dwMu0dtdwb
      %  wMu0dt
      %  wMu0dtb
      %  str=input('Look at drift');
        c0=wMu0dt+w0;
        c(1)=sigma0*sqrt(ds(tt));
        c(2:10)=0.0;
        w0=c0;
    
    elseif(tt<=6) 

    %dt=ds(tt);    
        
    %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
     %[b0,b] = CalculateDriftbCoeffs08O4(wMu0dt,dwMu0dtdw,c);
     [b0,b] = CalculateDriftbCoeffs08(wMu0dt,dwMu0dtdw,c,gamma);
     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);
     
     w0=c0;
            
    else
        
      
        
        
        
        
    %below wMu0dt is drift and dwMu0dtdw are its derivatives.
     [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),6);
     [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
     [b0,b] = CalculateDriftbCoeffs08(wMu0dt,dwMu0dtdw,c,gamma);
     
     size(b)
     ba(1:6)=b(1:6);
     
     
    [wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),8);
    
    [wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),8);
    
    [wVol3dt,dwVol3dtdw] = BesselVolAndDerivativesH3(w0,YqCoeff0,Fp1,gamma,ds(tt),8);

    [g10,g1] = CalculateDriftbCoeffs08O4(wVol0dt,dwVol0dtdw,c);
    
    g1(5:6)=0;
    
    
    [g20,g2] = CalculateDriftbCoeffs08O4(wVol2dt,dwVol2dtdw,c);
 g2(5:6)=0;
 
    [g30,g3] = CalculateDriftbCoeffs08O4(wVol3dt,dwVol3dtdw,c);
 
 g3(5:6)=0;
 
 g40=0;
 g4(1:6)=0;

 a0=c0+b0;
 
 a(1:5)=c(1:5)+b(1:5);
    
     [a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A(a0,a,g10,g1,g20,g2,g30,g3,g40,g4,sigma0,ds(tt));
     %[a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A(a0,a,b0,ba,g0,g,h0,h,sigma0,ds(tt));
     
     
 c0=a0;
 c(1)=a1;
 c(2)=a2;
 c(3)=a3;
 c(4)=a4;
 c(5)=a5;
    
    w0=c0;
    
    end

end
wnStart=1;


yy0=((1-gamma).*w0).^(1/(1-gamma));

 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=ceil(Nn/2)-1:-1:1
     if((w(nn)<0)&&(Flag==0))
         wnStart=nn+1;
         Flag=1;
     end
 end
 
 
 %c0
 %c
 
 
 %str=input('Look at numbers')
 yy0(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));

 yy=yy0;
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

ItoHermiteMean=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %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


theta1=1;
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=30;
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(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');
 
title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
 
legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')
 
str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
 

end
.
.
2nd function.
%function [a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A(a0,a,c0,c,g0,g,h0,h,sigma0,dt)
function [a0,a1,a2,a3,a4,a5] = AdvanceSeriesCumulantSolution02NewO6Params5A(a0,a,g10,g1,g20,g2,g30,g3,g40,g4,sigma0,dt)


a1=a(1);
a2=a(2);
a3=a(3);
a4=a(4);
a5=a(5);
[F] = CalculateSixCumulantsParamsFive(a0,a1,a2,a3,a4,a5);

C1=F(1,1);
C2=F(2,1);
C3=F(3,1);
C4=F(4,1);
C5=F(5,1);
C6=F(6,1);

g10=g10+sigma0*sqrt(dt);
%[dC2,dC3,dC4] = CalculateCumulantsHighOrder(g0,g,h0,h);

 SeriesOrder=6;
      ab00=0;
      ab0(1:SeriesOrder)=0;
%      b30=0;
%      b3(1:SeriesOrder)=0;
%      b40=0;
%      b4(1:SeriesOrder)=0;
 [dC1,dC2,dC3,dC4,dC5,dC6] = CalculateDiffusionCumulants(ab00,ab0,g10,g1,g20,g2,g30,g3,g40,g4,4,6);
 %[dC1,dC2,dC3,dC4,dC5,dC6] = CalculateDiffusionCumulants(ab00,ab0,g0,g,h0,h,b30,b3,b40,b4,4,6);


%dC1=0;
C1=C1+dC1;
C2=C2+dC2;
C3=C3+dC3;
C4=C4+dC4;
C5=C5+dC5;
C6=C6+dC6;



a1=a(1);
a2=a(2);
a3=a(3);
a4=a(4);
a5=a(5);

DZa0=a(1);
DZa(1)=2*a(2);
DZa(2)=3*a(3);
DZa(3)=4*a(4);
DZa(4)=5*a(5);

[DZaI0,DZaI] = SeriesReciprocal(DZa0,DZa,4);


a1=a1+.5*g10^2*DZaI0;
a2=a2+.5*g10^2*DZaI(1)/1;
a3=a3+.5*g10^2*DZaI(2)/1;
a4=a4+.5*g10^2*DZaI(3)/1;
a5=a5+.5*g10^2*DZaI(4)/1;

% a0
% a1
% a2
% a3
% a4
% a5
%str=input('Look at coefficients');
da(1,1)=a0;
da(2,1)=a1;
da(3,1)=a2;
da(4,1)=a3;
da(5,1)=a4;
da(6,1)=a5;

a(1)=a1;
a(2)=a2;
a(3)=a3;
a(4)=a4;
a(5)=a5;

%a0=0;
%[Fa,dFa] = CalculateSixCumulantsParamsFiveAndDerivs(C1,C2,C3,C4,C5,C6,a0,a1,a2,a3,a4,a5);
%Fa
%[F,dF] = CalculateCumulantsAndDerivativesFromMoments(C1,C2,C3,C4,C5,C6,a0,a,5,6,6);
[F,dF] = CalculateCumulantsAndDerivativesFromMoments(C1,C2,C3,C4,C5,C6,a0,a1,a2,a3,a4,a5,5,6,6);
%F-Fa
%dF-dFa
%str=input('Look at coefficients');
%[F,dF] = CalculateCumulantsAndDerivativesFromMoments(C1,C2,C3,C4,C5,C6,a0,a,SeriesOrder,NZterms,NMoments)


ObjBest=F(1,1)^2+F(2,1)^2+F(3,1)^2+F(4,1)^2+F(5,1)^2+F(6,1)^2;

a0Best=a0;
a1Best=a1;
a2Best=a2;
a3Best=a3;
a4Best=a4;
a5Best=a5;

nn=0;
while((nn<5)&&((abs(F(1,1))>.000000000001) || (abs(F(2,1))>.000000000001) || (abs(F(3,1))>.000000000001) || (abs(F(4,1))>.00000000001)|| (abs(F(5,1))>.0000000001)|| (abs(F(6,1))>.000000001)))
nn=nn+1;
da=da-dF\F;

%inv(dF)

a0=da(1,1);
a1=da(2,1);
a2=da(3,1);
a3=da(4,1);
a4=da(5,1);
a5=da(6,1);

%[F,dF] = CalculateCumulantsAndDerivatives(dC1,dC2,dC3,dC4,a1,a2,a3,a4,da0,da1,da2,da3,da4);
%[F,dF] = CalculateSixCumulantsParamsFiveAndDerivs(C1,C2,C3,C4,C5,C6,a0,a1,a2,a3,a4,a5);
%function [F,dF] = CalculateCumulantsAndDerivativesFromMoments(C1,C2,C3,C4,C5,C6,a0,a1,a2,a3,a4,a5,SeriesOrder,NZterms,NMoments)
[F,dF] = CalculateCumulantsAndDerivativesFromMoments(C1,C2,C3,C4,C5,C6,a0,a1,a2,a3,a4,a5,5,6,6);

 ObjNew=F(1,1)^2+F(2,1)^2+F(3,1)^2+F(4,1)^2+F(5,1)^2+F(6,1)^2;
 
  if(ObjBest<ObjNew)
      a0=a0Best;
      a1=a1Best;
      a2=a2Best;
      a3=a3Best;
      a4=a4Best;
      a5=a5Best;
  else
      ObjBest=ObjNew;
      a0Best=a0;
      a1Best=a1;
      a2Best=a2;
      a3Best=a3;
      a4Best=a4;
      a5Best=a5;
  end

da(1,1)=a0;
da(2,1)=a1;
da(3,1)=a2;
da(4,1)=a3;
da(5,1)=a4;
da(6,1)=a5;



end
  a0=a0Best;   
  a1=a1Best;
  a2=a2Best;
  a3=a3Best;
  a4=a4Best;
  a5=a5Best;

  
% a0
% a1
% a2
% a3
% a4
nn
%str=input('Look at results');

end

.
.
.
function [F,dF] = CalculateCumulantsAndDerivativesFromMoments(C1,C2,C3,C4,C5,C6,a0,a1,a2,a3,a4,a5,SeriesOrder,NZterms,NMoments)


a(1)=a1;
a(2)=a2;
a(3)=a3;
a(4)=a4;
a(5)=a5;


aa0=a0;
a0=0;

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:50)=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:50)=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);
u5=EXZ(6,1);
u6=EXZ(7,1);


% u1
% u2
% u3
% u4
% u5
% u6

%str=input('Look at numbers')

du1(1)=0;
du2(1)=0;
du3(1)=0;
du4(1)=0;
du5(1)=0;
du6(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);
du5(mm)=5*EXZ(5,mm);
du6(mm)=6*EXZ(6,mm);
end

%du1(1)=1;



k1=aa0+a(2)+3*a(4);
%k1=u1;
k2=u2-u1^2;
k3=u3-3*u2*u1+2*u1^3;
k4=u4-4*u3*u1-3*u2^2+12*u2*u1^2-6*u1^4;
k5=u5-5*u4*u1-10*u3*u2+20*u3*u1^2+30*u2^2*u1-60*u2*u1^3+24*u1^5;
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;


% k1
% k2
% k3
% k4
% k5
% k6

% str=input('Look at numbers--2')
dk1(1)=1.0;
dk1(2)=0.0;
dk1(3)=1.0;
dk1(4)=0.0;
dk1(5)=3.0;
dk1(6)=0.0;

dk2(1)=0;
dk3(1)=0;
dk4(1)=0;
dk5(1)=0;
dk6(1)=0;

for mm=2:SeriesOrder+1
%dk1(mm)=du1(mm);
dk2(mm)=-2.* u1.*du1(mm)+du2(mm);
dk3(mm)=du3(mm)-3*du1(mm)*u2-3*u1*du2(mm)+2*3*u1^2*du1(mm);
dk4(mm)=du4(mm)-4*du3(mm)*u1-4*u3*du1(mm)-3*2*u2*du2(mm)+12*du2(mm)*u1^2+12*2*u2*u1*du1(mm)-6*4*u1^3*du1(mm);
dk5(mm)=du5(mm)-5*du4(mm)*u1-5*u4*du1(mm)-10*du3(mm)*u2-10*u3*du2(mm)+20*du3(mm)*u1^2+20*2*u3*u1*du1(mm)+ ...
    30*2*u2*du2(mm)*u1+30*u2^2*du1(mm)-60*du2(mm)*u1^3-60*3*u2*u1^2*du1(mm)+24*5*u1^4*du1(mm);
dk6(mm)=du6(mm)-6*du5(mm)*u1-6*u5*du1(mm)-15*du4(mm)*u2-15*u4*du2(mm)+30*du4(mm)*u1^2+30*2*u4*u1*du1(mm)-10*2*u3*du3(mm)+ ...
    120*du3(mm)*u2*u1-120*u3*du2(mm)*u1-120*u3*u2*du1(mm)-120*du3(mm)*u1^3-120*3*u3*u1^2*du1(mm)+30*3*u2^2*du2(mm) - ...
    270*2*u2*du2(mm)*u1^2+270*2*u2^2*u1*du1(mm)+360*du2(mm)*u1^4+360*4*u2*u1^3*du1(mm)-120*6*u1^5*du1(mm);

end





F(1,1)=k1-C1;
 F(2,1)=k2-C2;
 F(3,1)=k3-C3;
 F(4,1)=k4-C4;
 F(5,1)=k5-C5;
 F(6,1)=k6-C6;

 
     for mm=1:SeriesOrder+1
        dF(1,mm)=dk1(mm);
        dF(2,mm)=dk2(mm);
        dF(3,mm)=dk3(mm);
        dF(4,mm)=dk4(mm);
        dF(5,mm)=dk5(mm);
        dF(6,mm)=dk6(mm);
     end

end

.
.
.
function [dC1,dC2,dC3,dC4,dC5,dC6] = CalculateDiffusionCumulants(d00,d0,d10,d1,d20,d2,d30,d3,d40,d4,SeriesOrder,NMoments)



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;

%[dC2,dC3,dC4] = CalculateCumulantsHighOrder(g0,g,h0,h)

a0=d0-d2+3*d4;
a00=d00-d20+3*d40;
a0(SeriesOrder+1:SeriesOrder*NMoments)=0;

%size(d1)
%size(d3)

a1=d1-3*d3;
a10=d10-3*d30;
a1(SeriesOrder+1:SeriesOrder*NMoments)=0;

a2=d2-6.*d4;
a20=d20-6.*d40;
a2(SeriesOrder+1:SeriesOrder*NMoments)=0;

a3=d3;
a30=d30;
a3(SeriesOrder+1:SeriesOrder*NMoments)=0;

a4=d4;
a40=d40;
a4(SeriesOrder+1:SeriesOrder*NMoments)=0;







%                     str=input('Look at indices');
T0(1:NMoments,1:NMoments*SeriesOrder)=0;                    
T0(1,1:SeriesOrder)=a0(1:SeriesOrder);
T00(1)=a00;

T1(1:NMoments,1:NMoments*SeriesOrder)=0;                    
T1(1,1:SeriesOrder)=a1(1:SeriesOrder);
T10(1)=a10;

T2(1:NMoments,1:NMoments*SeriesOrder)=0; 
T2(1,1:SeriesOrder)=a2(1:SeriesOrder);
T20(1)=a20;

T3(1:NMoments,1:NMoments*SeriesOrder)=0; 
T3(1,1:SeriesOrder)=a3(1:SeriesOrder);
T30(1)=a30;

T4(1:NMoments,1:NMoments*SeriesOrder)=0; 
T4(1,1:SeriesOrder)=a4(1:SeriesOrder);
T40(1)=a40;
for k = 2 : (NMoments)                    
b=T0(k-1,1:SeriesOrder*k);
b0=T00(k-1);
[b0,b] =SeriesProduct(a00,a0,b0,b,SeriesOrder*k);                    
T0(k,1:SeriesOrder*k)=b(1:SeriesOrder*k);
T00(k)=b0;  
    
    
b=T1(k-1,1:SeriesOrder*k);
b0=T10(k-1);
[b0,b] =SeriesProduct(a10,a1,b0,b,SeriesOrder*k);                    
T1(k,1:SeriesOrder*k)=b(1:SeriesOrder*k);
T10(k)=b0;  

b=T2(k-1,1:SeriesOrder*k);
b0=T20(k-1);
[b0,b] =SeriesProduct(a20,a2,b0,b,SeriesOrder*k);                    
T2(k,1:SeriesOrder*k)=b(1:SeriesOrder*k);  
T20(k)=b0;

b=T3(k-1,1:SeriesOrder*k);
b0=T30(k-1);
[b0,b] =SeriesProduct(a30,a3,b0,b,SeriesOrder*k);                    
T3(k,1:SeriesOrder*k)=b(1:SeriesOrder*k);  
T30(k)=b0;

b=T4(k-1,1:SeriesOrder*k);
b0=T40(k-1);
[b0,b] =SeriesProduct(a40,a4,b0,b,SeriesOrder*k);                    
T4(k,1:SeriesOrder*k)=b(1:SeriesOrder*k);  
T40(k)=b0;

end

CC(1,1,1,1,1)=1;


Moment1(1:NMoments)=0;
for k = 0 : (NMoments)
    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;
                for p = 0:j
                    l1=j-p+1;
                    l_0=p+1;
                    if(k>0)
                        b(1:SeriesOrder*k)=0;
                        b0=1;
                        %CC(l_0,l1,l2,l3,l4)=1.0;
                        if(l_0>1)
                            a=T0(l_0-1,1:SeriesOrder*k);
                            a0=T00(l_0-1);
                            [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*k);
                            %CC(l_0,l1,l2,l3,l4)=CC(l_0,l1,l2,l3,l4).*4^(l_0-1);
                        end  
                        if(l1>1)
                            a=T1(l1-1,1:SeriesOrder*k);
                            a0=T10(l1-1);
                            [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*k);
                            %CC(l_0,l1,l2,l3,l4)=CC(l_0,l1,l2,l3,l4).*2^(l1-1);
                        end    
                        if(l2>1)
                            a=T2(l2-1,1:SeriesOrder*k);
                            a0=T20(l2-1);
                            [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*k);
                            %CC(l_0,l1,l2,l3,l4)=CC(l_0,l1,l2,l3,l4).*2^(l2-1);
                        end
                        if(l3>1)
                            a=T3(l3-1,1:SeriesOrder*k);
                            a0=T30(l3-1);
                            [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*k);
                            %CC(l_0,l1,l2,l3,l4)=CC(l_0,l1,l2,l3,l4).*0^(l3-1);
                        end
                        if(l4>1)
                            a=T4(l4-1,1:SeriesOrder*k);
                            a0=T40(l4-1);
                            [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*k);
                            %CC(l_0,l1,l2,l3,l4)=CC(l_0,l1,l2,l3,l4).*0^(l4-1);
                        end

                      %  CC(l_0,l1,l2,l3,l4)=CC(l_0,l1,l2,l3,l4).*4^(l_0-1).*2^(l1-1)*2^(l2-1)*0^(l3-1)*0^(l4-1);
                        CC(l_0,l1,l2,l3,l4)=b0;
                        for qq=1:SeriesOrder*k
                            CC(l_0,l1,l2,l3,l4)=CC(l_0,l1,l2,l3,l4)+b(qq)*EZ(qq);
                        end
                        
                        
%                     Moment1(k)=Moment1(k)+factorial(k)/(factorial(l_0-1).*factorial(l1-1).*factorial(l2-1).*factorial(l3-1).*factorial(l4-1)).*c0(l_0).*CC(l1,l2,l3,l4);%.*EZ((l1-1)+2*(l2-1)+3*(l3-1)+4*(l4-1));
%                     k
%                     l_0
%                     l1
%                     l2
%                     l3
%                     l4
%                     
%                     str=input('Look at indices');
                    end
                end
            end
        end
    end
end



Moment1(1:NMoments)=0;
for k = 0 : (NMoments)
    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;
                for p = 0:j
                    l1=j-p+1;
                    l_0=p+1;
                    if((k>0)&&((l_0+l1+l2+l3+l4)>5))
                        if((l1+l2+l3+l4==4))
                             Moment1(k)=Moment1(k)+factorial(k)/(factorial(l_0-1).*factorial(l1-1).*factorial(l2-1).*factorial(l3-1).*factorial(l4-1)).*CC(l_0,l1,l2,l3,l4);%.*EZ((l1-1)+2*(l2-1)+0*3*(l3-1)+0*4*(l4-1));
                        else
                             Moment1(k)=Moment1(k)+factorial(k)/(factorial(l_0-1).*factorial(l1-1).*factorial(l2-1).*factorial(l3-1).*factorial(l4-1)).*CC(l_0,l1,l2,l3,l4).*EZ((l1-1)+2*(l2-1)+3*(l3-1)+4*(l4-1));
                        end
                        %                     k
                    % l_0
                    % l1
                    % l2
                    % l3
                    % l4
                    %Moment1(k)=Moment1(k)+factorial(k)/(factorial(l_0-1).*factorial(l1-1).*factorial(l2-1).*factorial(l3-1).*factorial(l4-1)).*CC(l_0,l1,l2,l3,l4);%.*EZ((l1-1)+2*(l2-1)+0*3*(l3-1)+0*4*(l4-1));
%                     k
%                     l_0
%                     l1
%                     l2
%                     l3
%                     l4
%                     
%                     str=input('Look at indices');
                    end
                end
            end
        end
    end
end

%Moment1
%str=input('Look at Moments');

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



dC1=u1;
dC2=u2-u1^2;
dC3=u3-3*u1*u2+2*u1^3;
dC4=u4-4*u3*u1-3*u2^2+12*u2*u1^2-6*u1^4;
dC5=u5-5*u4*u1-10*u3*u2+20*u3*u1^2+30*u2^2*u1-60*u2*u1^3+24*u1^5;
dC6=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;
%dC5=0;
%dC6=0;




end

.
.
.
function [F] = CalculateSixCumulantsParamsFive(a0,a1,a2,a3,a4,a5)
%[F,dF] = CalculateCumulantsAndDerivatives(dC1,dC2,dC3,dC4,a1,a2,a3,a4,da0,da1,da2,da3,da4);

k1=a0+a2+3*a4;

k2=a1^2+2*a2^2+6*a1*a3+15*a3^2+24*a2*a4+96*a4^2+ ... 
30*(a1+7*a3)*a5+945*a5^2;
 
k3=8*a2^3+216*a2^2*a4+6*a1^2*(a2+6*a4)+72*a1*a3*(a2+8*a4)+ ...
540*a1*(a2+10*a4)*a5+ ...
18*a2*(15*a3^2+128*a4^2+280*a3*a5+1575*a5^2)+ ...
108*a4*(25*a3^2+88*a4^2+560*a3*a5+3675*a5^2);
 
k4=24*(2*a2^4+96*a2^3*a4+a1^3*(a3+10*a5)+ ...
30*a2^2*(6*a3^2+64*a4^2+140*a3*a5+945*a5^2)+ ...
a1^2*(2*(a2^2+9*a3^2+18*a2*a4+96*a4^2)+375*a3*a5+ ...
2250*a5^2)+ ...
36*a2*a4*(125*a3^2+528*a4^2+3360*a3*a5+25725*a5^2)+ ...
9*(45*a3^4+3560*a3^2*a4^2+8704*a4^4+ ...
35*a3*(69*a3^2+3104*a4^2)*a5+ ...
700*(79*a3^2+1332*a4^2)*a5^2+632625*a3*a5^3+ ...
3018750*a5^4)+ ...
3*a1*(12*a2^2*(a3+10*a5)+8*a2*a4*(32*a3+375*a5)+ ...
3*(15*a3^3+528*a3*a4^2+530*a3^2*a5+7120*a4^2*a5+ ...
7175*a3*a5^2+36750*a5^3)));

k5=24*(16*a2^5+5*a1^4*a4+1200*a2^4*a4+ ...
10*a1^3*(3*a2*a3+42*a3*a4+40*a2*a5+570*a4*a5)+ ...
300*a2^3*(10*a3^2+128*a4^2+280*a3*a5+2205*a5^2)+ ...
1080*a2^2*a4*(125*a3^2+616*a4^2+3920*a3*a5+34300*a5^2)+ ...
30*a1*(a3*(16*a2^3+225*a2*a3^2+640*a2^2*a4+3330*a3^2*a4+ ...
9504*a2*a4^2+52224*a4^3)+ ...
10*(20*a2^3+954*a2*a3^2+900*a2^2*a4+15735*a3^2*a4+ ...
14952*a2*a4^2+91296*a4^3)*a5+ ...
3675*a3*(41*a2+750*a4)*a5^2+15750*(56*a2+1129*a4)*a5^3)+ ...
27*a4*(14775*a3^4+507200*a3^2*a4^2+950784*a4^4+ ...
700*a3*(1469*a3^2+27904*a4^2)*a5+ ...
1750*(16973*a3^2+118128*a4^2)*a5^2+419002500*a3*a5^3+ ...
2419134375*a5^4)+ ...
10*a1^2*(2*a2^3+72*a2^2*a4+ ...
3*a2*(24*a3^2+320*a4^2+625*a3*a5+4500*a5^2)+ ...
9*a4*(109*a3^2+528*a4^2+3130*a3*a5+24975*a5^2))+ ...
90*a2*(270*a3^4+69632*a4^4+16905*a3^3*a5+8391600*a4^2*a5^2+ ...
30187500*a5^4+280*a3^2*(89*a4^2+1580*a5^2)+ ...
35*a3*(24832*a4^2*a5+162675*a5^3)));

k6=120*(32*a2^6+3456*a2^5*a4+6*a1^5*a5+ ...
3*a1^4*(9*a3^2+8*a4*(2*a2+21*a4)+330*a3*a5+2850*a5^2)+ ...
720*a2^4*(15*a3^2+224*a4^2+490*a3*a5+4410*a5^2)+ ...
6048*a2^3*a4*(125*a3^2+704*a4^2+4480*a3*a5+44100*a5^2)+ ...
12*a1^3*(12*a3*(a2^2+6*a3^2+35*a2*a4+298*a4^2)+ ...
10*(20*a2^2+351*a3^2+684*a2*a4+6132*a4^2)*a5+ ...
60450*a3*a5^2+374625*a5^3)+ ...
216*a2^2*(945*a3^4+99680*a3^2*a4^2+313344*a4^4+ ...
420*a3*(161*a3^2+9312*a4^2)*a5+ ...
25200*(79*a3^2+1665*a4^2)*a5^2+28468125*a3*a5^3+ ...
166031250*a5^4)+ ...
1296*a2*a4*(6*(985*a3^4+38040*a3^2*a4^2+79232*a4^4)+ ...
35*a3*(13221*a3^2+279040*a4^2)*a5+ ...
1925*(7715*a3^2+59064*a4^2)*a5^2+230451375*a3*a5^3+ ...
1451480625*a5^4)+ ...
18*a1*(4*a3*(20*a2^4+1280*a2^3*a4+ ...
27*a2^2*(25*a3^2+1232*a4^2)+ ...
6*a2*(3885*a3^2*a4+69632*a4^3)+ ...
9*(135*a3^4+24280*a3^2*a4^2+237696*a4^4))+ ...
15*(80*a2^4+5600*a2^3*a4+56*a2^2*(159*a3^2+2848*a4^2)+ ...
64*a2*(5245*a3^2*a4+34236*a4^3)+ ...
3*(9717*a3^4+1145552*a3^2*a4^2+4076032*a4^4))*a5+ ...
2100*a3*(1148*a2^2+8199*a3^2+47250*a2*a4+ ...
526740*a4^2)*a5^2+ ...
126000*(126*a2^2+2935*a3^2+5645*a2*a4+68241*a4^2)*a5^3+ ...
4314397500*a3*a5^4+21772209375*a5^5)+ ...
27*(9945*a3^6+2894760*a3^4*a4^2+61140480*a3^2*a4^4+ ...
92930048*a4^6+ ...
210*a3*(5553*a3^4+1177072*a3^2*a4^2+13526272*a4^4)*a5+ ...
2100*(29649*a3^4+4101044*a3^2*a4^2+17011776*a4^4)*a5^2+ ...
3500*a3*(551015*a3^2+41249952*a4^2)*a5^3+ ...
39375*(923683*a3^2+24833432*a4^2)*a5^4+ ...
394058306250*a3*a5^5+1907724656250*a5^6)+ ...
6*a1^2*(8*a2^4+480*a2^3*a4+ ...
180*a2^2*(4*a3^2+64*a4^2+125*a3*a5+1050*a5^2)+ ...
72*a2*a4*(327*a3^2+1848*a4^2+10955*a3*a5+99900*a5^2)+ ...
9*(225*a3^4+22824*a3^2*a4^2+69632*a4^4+15090*a3^3*a5+ ...
830240*a3*a4^2*a5+412925*a3^2*a5^2+8239800*a4^2*a5^2+ ...
5475750*a3*a5^3+29636250*a5^4)));




 F(1,1)=k1;
 F(2,1)=k2;
 F(3,1)=k3;
 F(4,1)=k4;
 F(5,1)=k5;
 F(6,1)=k6;
       
end

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


NoOfTerms=19;

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 [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 [wVol3dt,dwVol3dtdw] = BesselVolAndDerivativesH3(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder)


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

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



Fp2(1)=Fp1(1,1,2,4);
Fp2(2)=Fp1(1,2,1,4);
Fp2(3)=Fp1(2,1,1,4);


wVol3dt0=0;
dwVol3dt(1:ExpnOrder)=0.0;

wVol3dt=0;
dwVol3dtdw(1:ExpnOrder)=0.0;

for mm=1:NoOfTerms

    wVol3dt0=YqCoeffa(mm).*w0.^Fp2(mm);
    for nn=1:ExpnOrder
        if(nn==1)
            dwVol3dt(nn)=wVol3dt0*Fp2(mm)*1/w0;
        else
        dwVol3dt(nn)=dwVol3dt(nn-1)*(Fp2(mm)-(nn-1))/w0;
        end
    end
    wVol3dt=wVol3dt+wVol3dt0;
    for nn=1:ExpnOrder
        dwVol3dtdw(nn)=dwVol3dtdw(nn)+dwVol3dt(nn);
    end
        
    
end



end
.
.
.
function [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,ExpnOrder)



NoOfTerms=9;%excluding dt^3 terms

YqCoeff(1:NoOfTerms)=0.0;
Fp(1:NoOfTerms)=0.0;
a0=((beta1-gamma)/(1-gamma));
b0=((beta2-gamma)/(1-gamma));
A=mu1;
B=mu2;
C=-.5*sigma0^2*gamma;
dt2=dt^2/2;

YqCoeff(1:9)=0;

YqCoeff(1)=A*dt;
YqCoeff(2)=B*dt;
YqCoeff(3)=C*dt;

 YqCoeff(4)=A^2*a0*dt2*(1-gamma);
 YqCoeff(5)=B^2*b0*dt2*(1-gamma);
 YqCoeff(6)=1*1.5*.66*C^2*(-1)*dt2*(1-gamma)+1*1.5*.66* C*(-1)*(-2)*.5*sigma0^2 *dt2*(1-gamma)^2;
 YqCoeff(7)=A*B*(a0+b0)*dt2*(1-gamma);
 YqCoeff(8)=1*1.5*.66*A*C*(a0-1)*dt2*(1-gamma)^1+1*1.5*.66*.5*sigma0^2*A*a0*(a0-1)*dt2*(1-gamma)^2;
 YqCoeff(9)=1*1.5*.66*B*C*(b0-1)*dt2*(1-gamma)^1+1*1.5*.66*.5*sigma0^2*B*b0*(b0-1)*dt2*(1-gamma)^2;
 



Fp(1)=(beta1-gamma)/(1-gamma);
Fp(2)=(beta2-gamma)/(1-gamma);
Fp(3)=-1;

Fp(4)=2*(beta1-gamma)/(1-gamma)-1;
Fp(5)=2*(beta2-gamma)/(1-gamma)-1;
Fp(6)=-3;
Fp(7)=(beta1-gamma)/(1-gamma)+(beta2-gamma)/(1-gamma)-1;
Fp(8)=(beta1-gamma)/(1-gamma)-2;
Fp(9)=(beta2-gamma)/(1-gamma)-2;


 Ratio=2/dt*(YqCoeff(4)*((1-gamma)*w0).^Fp(4)+YqCoeff(5)*((1-gamma)*w0).^Fp(5)+YqCoeff(6)*((1-gamma)*w0).^Fp(6)+ ...
     YqCoeff(7)*((1-gamma)*w0).^Fp(7)+YqCoeff(8)*((1-gamma)*w0).^Fp(8)+YqCoeff(9)*((1-gamma)*w0).^Fp(9))/ ...
     (YqCoeff(1)*((1-gamma)*w0).^Fp(1)+YqCoeff(2)*((1-gamma)*w0).^Fp(2)+YqCoeff(3)*((1-gamma)*w0).^Fp(3));

Fp2(1:9)=Fp(1:9);%/(1-gamma);

wMu0dt0=0;
dwMu0dt(1:ExpnOrder)=0.0;

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

for mm=1:NoOfTerms

    wMu0dt0=YqCoeff(mm).*((1-gamma)*w0).^Fp2(mm);
    for nn=1:ExpnOrder
        if(nn==1)
            dwMu0dt(nn)=wMu0dt0*Fp2(mm)*1/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,gamma)



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


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


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




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

b(7:10)=0.0; 

Bound=abs(b(1)/b0)*gamma/(1-gamma);
%Bound=1/abs(b0/b(1))*gamma;%/(1-gamma);

for nn=1:9
   
    if(abs(b(nn+1))>abs(b(nn)*Bound))
  %      b(nn+1)=sign(b(nn+1))*abs(b(nn))*Bound;    
    end
end

%[b0,b] = ApplyBoundsOnSeriesCoeffs(b0,b);


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

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



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



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: 1839
Joined: July 14th, 2002, 3:00 am

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

March 22nd, 2022, 3:51 pm

Friends, mind control agents kept me awake almost all night (last night) and let me have a very poor sleep since I was supposed to travel from kot addu to Lahore early in the morning. When I travel for a few hundred kilometers, many times my journey becomes quite stimulating and I think of new ideas about my research and in order to thwart that, mind control agents did not let me sleep properly last night. We started very early in the morning and It was very difficult to drive in the early part of the journey. 
There is some possibility I might have corona. Day before yesterday, my father complained that he had dry cough with irritation in the throat and I brought medicine for him. I was thinking it was minor flu but today after coming to Lahore, I had similar dry cough with irritation in the throat. And my father already have 100 degrees farenheit fever. I think I have caught his germs and I am on his footsteps to get fever. I am not very worried since I (and also my father) had two doses of vaccine and also a third booster dose but I was feeling a bit weak during the day. I wanted to write explanation of major new functions in my program that I posted yesterday but I was feeling weak and could not get myself to do it. What is very disturbing that mind control agents continued to try to drug my food today and I had to be very careful. I was also injected with antipsychotic fluanxol injection around noon after arriving in Lahore. At the time of injection, I did not have any cough and if it were clear at that time that I was infected, I would have resisted the injection and might possibly have delayed it but I did not know. Mind control agents have no hesitation to drug my food and they want to use my corona infection as an opportunity to damage my brain and control the neurotransmitters they do not know. These agents never become better until they (or their backers) are reprimanded and scolded by good people. I will continue to post about their efforts to drug my food and not let me have peace (by keeping me awake etc.) and other attempts to damage my brain. 
I will also try to write explanation of the main functions of yesterday's program tomorrow if I could muster enough effort to do that.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
 
User avatar
Amin
Topic Author
Posts: 1839
Joined: July 14th, 2002, 3:00 am

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

March 22nd, 2022, 4:11 pm

Friends, after writing above post, I turned off the lights and wanted to try to sleep but mind control agents started giving me severe itch in my back. It is very disturbing and I have continued to beg people in the past to stop mind control agents from such lowly tactics but they continue to get bolder. I have quickly made this post and will try to sleep again and if they continued their tactics, I will make a post again in the morning.
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