Breakthrough in the theory of stochastic differential equations and their simulation

Amin · December 26th, 2021, 9:33 am

Friends I was not able to succesfully complete the project I mentioned. When I solved the linearized equations for change in cumulant values, I could not get the coefficients of Z any close to right. Apparently linearization was not a good idea but otherwise it might be too difficult to solve these equations. I have several other ideas that I want to develop fully to see if they can be useful in advancing coefficients of Z using the "monte carlo like equation" I will continue to work with these ideas and let the friends know of the results as they become available.
But I am writing this post to describe the tough torture I have faced in past few days. I had not written about it in past few days since all my energy was concentrated towards completing my work. Past few days have been a nightmare. I was so totally down because they keep taking neurotransmitters out of my brain. They have been using gases in my room, living room and in my car. They continue to use special gases in my car so that whenever I return home, I am totally energy-less and very tired. Mostly for past few days, I would simply sleep after coming back from driving outside in my car to get somewhat better.
A few days ago mind control agents asked my mother to give me drugged food since I had been having good food from outside. She especially cooked food (mutton chops) for me and I had to eat enough food. The food was badly drugged and I remained very down for several days after that. My walk is usually extremely easy for me but when I walked the next morning, I had to take several large breaks to sit next to the road because it was so difficult walking.
They continue to charge my eyes and many times it becomes very unnerving experience for me and I would have fits of extreme anxiety.
More than four months ago, I would sometimes work in drawing room of our house but they have consistently increased some heavy gas there that whenever I go inside the drawing room, my pants become charged and start to stick to my body and it feels very uncomfortable in the lower body.

Amin · December 30th, 2021, 10:47 am

I want to wish a wonderful new year to all friends and hope that this new year brings health and success to all of us.
I am very sure that on my thread we will remain more busy than ever during this new year. We have successfully developed several tools and it would be a pity if we do not find their intelligent applications in mathematical finance. Giving friends some idea, I would be working on the general stochastic volatility solution that I had started to work with a few months earlier but we were bogged down by problems. Then I want to find a simple semi-analytic type solution for stochastic volatility Cheyette model which is one of the most interesting interest rate models. In the meantime, I want to keep brainstorming how to easily solve for stochastic integrals of SDEs ( I have several ideas about it and I hope we will find a good solution to this problem pretty soon) So I hope that we will remain quite busy exploring new things during this coming year.

An important thing about previous work that I recall is that I have to upload a new version of SeriesReciprocal routine that I earlier distributed because the previous routine becomes unstable when the leading coefficient of the series is small. Frankly I am not an expert at series at all (though I dabbled with them before while doing time-dependent monte carlo that includes time exponentials). And I copied my series routines from math.stackexchange which is a wonderful website. In fact somebody there had mentioned that it would be helpful to convert the leading coefficient of series to unity when finding series reciprocal but when I started working with simple problems at start, my simple routine worked very well and all warnings went to back of my mind. I will be pasting my new routine for Series Reciprocal in the next post in a few hours. With this new routine, we would be able to find solution for larger range of parameters close to zero where series reciprocal was blowing with older routine. I am sure most experts in mathematics would already have realized my mistake but I am sure there are still many friends like me who were unaware of the problem so new routine would be helpful for them.

Lastly I want to tell friends that my persecution activity has remarkably decreased since my last post (even though mind control agents continue some bad things from time to time) and I am very thankful to good people who have taken a stand to force mind control agencies to end or decrease my persecution. If it were not for good people in US who have repeatedly tried to end my persecution in the past, I would never have been able to do anything meaningful. I was able to do interesting research since there were many times when my persecution activity remarkably decreased for long periods of time due to pressure from good friends in US. But I would still say that every time mind control agencies decreased my persecution activity, they increased it again after a few weeks or months to new levels since they did not like that I was doing interesting work. Please try to end my persecution forever. I am already 48 and it is very difficult for me to continue taking antipsychotics now as they take all energy out of my body so I will request friends to please force mind control agencies to stop giving me any antipsychotics. Other than trying to do good mathematics, my greatest wish is to live a life free of antipsychotics that are forced on me. But I am very thankful to friends who tried to help me in the past and tried to decrease mind control activity against me. I really owe them for their favors.

Amin · December 30th, 2021, 2:41 pm

Friends, this morning I thought of another way to solve for solution of densities of SDEs and stochastic integrals. I decided to share the basic details with friends now. Please pardon any errors.
To keep things simple, I suppose I am in Bessel SDE setup. A typical SDE in Bessel coordinates is

[$]dw(t) \, = \, \mu(w) dt \, + \, \sigma dZ(t)[$]

We can write the above SDE a bit differently as

[$]w(t_2) \, = \, w(t_1) \, + \mu(w) dt \, + \, \sigma dZ(t)[$]

where [$]t_2 = t_1 + \Delta t[$]

We want to find the Z-series representation for [$]w(t_2)[$] when we are given the Z-series representation of [$]w(t_1)[$]. Suppose Z-series representation of [$]w(t_1)[$] is known and is given as

[$]w(t_1)= a_0 \, + \, a_1 \, Z \,+ \, a_2 \, Z^2 \, + a_3 \, Z^3 + ...[$]

We can also find a Z-series representation of [$] \mu(w(t_1)) [$] which is given as

[$]\mu(w(t_1)) \Delta t= b_0 \, + \, b_1 \, Z \,+ \, b_2 \, Z^2 \, + b_3 \, Z^3 + ...[$]

We substitute the Z-series of [$]w(t_1)[$] and [$]\mu(w(t_1))[$] in SDE to find

[$]w(t_2) \, = \, \sum a_n Z^n + \, \sum b_n Z^n+ \, \sigma \sqrt{(\Delta t)} Z_1[$]

Please note that original gaussian at t_1 is represented as Z while innovation gausssian is represented as [$]Z_1[$]

Now using the fact the [$]Z[$] and [$] Z_1[$] are independent and suppose their standard normal densities are denoted respectively as [$]P(Z)[$] and [$]P(Z_1)[$], we can write the joint density by introducing repsective normal densities in the above expression for [$]w(t_2)[$] as

[$]P[w(t_2) | (Z,Z_1)]= \, \Big[ \, \sum a_n Z^n + \, \sum b_n Z^n+ \, \sigma \sqrt(\Delta t) Z_1 \Big] \, P(Z_1) \, P(Z) \,[$]

To represent the density of [$]w(t_2)[$] as a Z_series, we introduce a new gaussian variable such that [$]Z_2=Z+Z_1[$]
So [$]Z=Z_2 - Z_1[$] and we substitute it in the above equation to find as

[$]P[w(t_2) | (Z_2,Z_1)]= \, \Big[ \, \sum a_n {(Z_2-Z_1)}^n + \, \sum b_n {(Z_2-Z_1)}^n+ \, \sigma \sqrt(\Delta t) Z_1 \Big] \, P(Z_1) \, P(Z_2-Z_1) \,[$]

We integrate over [$]Z_1[$] in the above joint density to find only an expression of marginal density of [$]w(t_2)[$] in terms of a Z-series in [$]Z_2[$]. This is done as

[$]P[w(t_2) | (Z_2)]=\, \int_{-\infty}^{+\infty} \Big[ \, \sum a_n {(Z_2-Z_1)}^n + \, \sum b_n {(Z_2-Z_1)}^n+ \, \sigma \sqrt(\Delta t) Z_1 \Big] \, P(Z_1) \, P(Z_2-Z_1) \,dZ_1 \,[$]

After integrating over [$]Z_1[$], we can represent the above expression in the form (In a few hours, I will show a detailed derivation of different steps of this integration and conversion into form in next equation in a new post)

[$]P[w(t_2) | (Z_2)]=\, \Big[\, c_0 \, + \, c_1 \, Z_2 \,+ \, c_2 \, {Z_2}^2 \, + c_3 \, {Z_2}^3 + ...\Big] \, P(Z_2) \, [$],
we can easily identify the Z_series of [$]w(t_2)[$] in above expression given as

[$]\, w(t_2)=\, c_0 \, + \, c_1 \, Z_2 \,+ \, c_2 \, {Z_2}^2 \, + c_3 \, {Z_2}^3 + ...\,[$]

So we have a simple "monte carlo-like" algorithm to advance the density of an SDE and its stochastic integrals.
.
.
Copying the related equation from above, I think we can solve this equation by using additive property of central moments/cumulants and then mapping the resulting moments on our coefficients. I have not worked out everything but it should be relatively simple. I will come back with details in a day or two. Here is the equation on which we can use additive property of cumulants.

[$]w(t_2) \, = \, \sum a_n Z^n + \, \sum b_n Z^n+ \, \sigma \sqrt{(\Delta t)} Z_1[$]

briefly there are two independent gaussian random variables [$]Z[$] and [$]Z_1[$].
We can write the above equation as
[$]w(t_2) \, = \, X_0 \,+ \, X_1[$]

where
[$]X_0 \,=   \, \sum a_n Z^n + \, \sum b_n Z^n [$]
and
[$]X_1 \,=   \, \sigma \sqrt{(\Delta t)} Z_1 [$]

For example mth raw moments would be

[$]E[{(X_0)}^m]=E[{(\, \sum a_n Z^n + \, \sum b_n Z^n)}^m ][$] [$]Z[$] is a standard normal.
[$]E[{(X_1)}^m]=E[{( \, \sigma \sqrt{(\Delta t)} Z_1)}^m ][$] [$]Z_1[$] is a standard normal.

which can easily be calculated since nth coefficients and moments of standard normal are known to us. Only thing we need is an algorithm to be able to project the added cumulants back to coefficients of powers of Z which I think can be done iteratively starting from smaller moments/cumulants and moving to higher moments/cumulants. At every step [$]c_0[$] would give us value of median where we could do expansion for drift function (and volatility function if we are not in bessel coordinates) and continue to use the expansion to solve for value of coefficients on next step.
The above holds for SDEs in Bessel coordinates but finding a solution would be more difficult in original coordinates and for dependent functions but I am very sure it can be done.
.
Friends my analytic calculations for solution of densities of SDEs by cumulant matching is going well. I want to explain some simple details of my suggested method.
Suppose initial SDE variable at start of simulation time step is given as
[$]X_0 \,=   \, \sum a_n Z^n + \, \sum b_n Z^n = \, \sum c_n Z^n [$] where [$]c_n \,=\, a_n \,+ \,b_n[$]

Initially I am working with fourth order so our initial SDE variable is given as

[$]X_0 \,= c_0\, + \, c_ 1 \, Z + \, c_ 2 \, Z^2 + \, c_ 3 \, Z^3 + \, c_ 4 \, Z^4[$]

I suppose that after a small time step the coefficients have changed by a small value each

[$]X_2 \,= (c_0\, + dc_0) + \,( c_ 1 +dc_1) \, Z + \,( c_ 2 +dc_2) \, Z^2 + \, ( c_ 3 +dc_3) \, Z + \, (c_ 4+dc_4) \, Z^4[$]

I write the cumulants before the simulation step [$]C_n(X_0)[$] and after the simulation step [$]C_n(X_2)[$] from the coefficients of Z in above two respective equation and match them with change in cumulant values due to addition of [$]X_1[$] as

[$]C_n(X_2) \, = \, C_n(X_0) \, + C_n(X_1) \,[$]

I suppose that [$]dc_2[$], [$]dc_3[$] and [$]dc_4[$] are very small for a small time step and their squares and higher powers can be neglected. I retain the square of leading variance coefficient [$]dc_1[$] that is change in coefficient of Z. It gives me system of equations of sort

[$]dc_0 + dc_2 + 3dc_4=C_1(X_1)[$]
[$]A_1 \, {(dc_1)}^2 \, + \, B_1 \, dc_1 \, + \, C_1 \, dc_2 \,  + \, D_1 \, dc_3 \,  + \, E_1 \, dc_4 \,= \, C_2(X_1)[$]
[$]A_2 \, {(dc_1)}^2 \, + \, B_2 \, dc_1 \, + \, C_2 \, dc_2 \,  + \, D_2 \, dc_3 \,  + \, E_2 \, dc_4 \,= \, C_3(X_1)[$]
[$]A_3 \, {(dc_1)}^2 \, + \, B_3 \, dc_1 \, + \, C_3 \, dc_2 \,  + \, D_3 \, dc_3 \,  + \, E_3 \, dc_4 \,= \, C_4(X_1)[$]

In above syatem of equations, [$]A_n, \, B_n, \, C_n, D_n , \, E_n [$] are coefficients that are already known(These are (polynomial-like) expressions in terms of [$]c_1, \, c_2, \, c_3, \, c_4[$] which are known at start of simulation step and their values can be calculated at start of simulation. The above system of equations can easily be solved in mathematica and this solution can be imported into matlab or C++.
Just to elaborate a bit further, I am copying the second of the above system of equations based on difference in second cumulant from mathematica as
[$] ({dc_1})^2 +(2 c_1 + 6 c_3) dc_1 + (4 c_2 + 24 c_4) dc_2 + (6 c_1 + 30 c_3) dc_3 + (24 c_2 + 192 c_4) dc_4= \, C_2(X_1)[$]
meaning [$]A_1=1[$], [$]B_1=(2 c_1 + 6 c_3) [$], [$]C_1=(4 c_2 + 24 c_4)[$] , [$]D_1\,=\,(6 c_1 + 30 c_3) [$] and [$]E_1=\,  (24 c_2 + 192 c_4)[$]

However we have to exogeneously specify the value of [$]dc_0[$] that is change in median. I could not find a simple analytic way to solve for median and I used leading order of fokker-planck equation solution for that and it is a very simple equation in leading coefficients. I will try to find some other easy way to solve for median and continue to think about it later. I have still not get the analytic solution to numbers but I hope to provide more details and post a working program in another two days.

.
https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1305#p868904

May be friends who are experts at neural networks might want to find the solution by root-finding(of system of equations) and then train neural networks on it. Can be a very interesting problem to work.

Amin · December 30th, 2021, 4:45 pm

https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1305#p868904

Friends, I have decided to give above problem a try again by inverting for parameters by matrix version of Newton iteration to find roots at every step. I am using matlab so it should be helpful to play with matrices. For the problem formulation I had, the condition number might be a problem in matrix inversion but I want to see how it goes. I will keep friends updated.

Amin · January 1st, 2022, 11:45 am

Friends, this morning I thought of another way to solve for solution of densities of SDEs and stochastic integrals. I decided to share the basic details with friends now. Please pardon any errors.
To keep things simple, I suppose I am in Bessel SDE setup. A typical SDE in Bessel coordinates is

[$]dw(t) \, = \, \mu(w) dt \, + \, \sigma dZ(t)[$]

We can write the above SDE a bit differently as

[$]w(t_2) \, = \, w(t_1) \, + \mu(w) dt \, + \, \sigma dZ(t)[$]

where [$]t_2 = t_1 + \Delta t[$]

We want to find the Z-series representation for [$]w(t_2)[$] when we are given the Z-series representation of [$]w(t_1)[$]. Suppose Z-series representation of [$]w(t_1)[$] is known and is given as

[$]w(t_1)= a_0 \, + \, a_1 \, Z \,+ \, a_2 \, Z^2 \, + a_3 \, Z^3 + ...[$]

We can also find a Z-series representation of [$] \mu(w(t_1)) [$] which is given as

[$]\mu(w(t_1)) \Delta t= b_0 \, + \, b_1 \, Z \,+ \, b_2 \, Z^2 \, + b_3 \, Z^3 + ...[$]

We substitute the Z-series of [$]w(t_1)[$] and [$]\mu(w(t_1))[$] in SDE to find

[$]w(t_2) \, = \, \sum a_n Z^n + \, \sum b_n Z^n+ \, \sigma \sqrt{(\Delta t)} Z_1[$]

Please note that original gaussian at t_1 is represented as Z while innovation gausssian is represented as [$]Z_1[$]

Now using the fact the [$]Z[$] and [$] Z_1[$] are independent and suppose their standard normal densities are denoted respectively as [$]P(Z)[$] and [$]P(Z_1)[$], we can write the joint density by introducing repsective normal densities in the above expression for [$]w(t_2)[$] as

[$]P[w(t_2) | (Z,Z_1)]= \, \Big[ \, \sum a_n Z^n + \, \sum b_n Z^n+ \, \sigma \sqrt(\Delta t) Z_1 \Big] \, P(Z_1) \, P(Z) \,[$]

To represent the density of [$]w(t_2)[$] as a Z_series, we introduce a new gaussian variable such that [$]Z_2=Z+Z_1[$]
So [$]Z=Z_2 - Z_1[$] and we substitute it in the above equation to find as

[$]P[w(t_2) | (Z_2,Z_1)]= \, \Big[ \, \sum a_n {(Z_2-Z_1)}^n + \, \sum b_n {(Z_2-Z_1)}^n+ \, \sigma \sqrt(\Delta t) Z_1 \Big] \, P(Z_1) \, P(Z_2-Z_1) \,[$]

We integrate over [$]Z_1[$] in the above joint density to find only an expression of marginal density of [$]w(t_2)[$] in terms of a Z-series in [$]Z_2[$]. This is done as

[$]P[w(t_2) | (Z_2)]=\, \int_{-\infty}^{+\infty} \Big[ \, \sum a_n {(Z_2-Z_1)}^n + \, \sum b_n {(Z_2-Z_1)}^n+ \, \sigma \sqrt(\Delta t) Z_1 \Big] \, P(Z_1) \, P(Z_2-Z_1) \,dZ_1 \,[$]

After integrating over [$]Z_1[$], we can represent the above expression in the form (In a few hours, I will show a detailed derivation of different steps of this integration and conversion into form in next equation in a new post)

[$]P[w(t_2) | (Z_2)]=\, \Big[\, c_0 \, + \, c_1 \, Z_2 \,+ \, c_2 \, {Z_2}^2 \, + c_3 \, {Z_2}^3 + ...\Big] \, P(Z_2) \, [$],
we can easily identify the Z_series of [$]w(t_2)[$] in above expression given as

[$]\, w(t_2)=\, c_0 \, + \, c_1 \, Z_2 \,+ \, c_2 \, {Z_2}^2 \, + c_3 \, {Z_2}^3 + ...\,[$]

So we have a simple "monte carlo-like" algorithm to advance the density of an SDE and its stochastic integrals.
.
.
Copying the related equation from above, I think we can solve this equation by using additive property of central moments/cumulants and then mapping the resulting moments on our coefficients. I have not worked out everything but it should be relatively simple. I will come back with details in a day or two. Here is the equation on which we can use additive property of cumulants.

[$]w(t_2) \, = \, \sum a_n Z^n + \, \sum b_n Z^n+ \, \sigma \sqrt{(\Delta t)} Z_1[$]

briefly there are two independent gaussian random variables [$]Z[$] and [$]Z_1[$].
We can write the above equation as
[$]w(t_2) \, = \, X_0 \,+ \, X_1[$]

where
[$]X_0 \,=   \, \sum a_n Z^n + \, \sum b_n Z^n [$]
and
[$]X_1 \,=   \, \sigma \sqrt{(\Delta t)} Z_1 [$]

For example mth raw moments would be

[$]E[{(X_0)}^m]=E[{(\, \sum a_n Z^n + \, \sum b_n Z^n)}^m ][$] [$]Z[$] is a standard normal.
[$]E[{(X_1)}^m]=E[{( \, \sigma \sqrt{(\Delta t)} Z_1)}^m ][$] [$]Z_1[$] is a standard normal.

which can easily be calculated since nth coefficients and moments of standard normal are known to us. Only thing we need is an algorithm to be able to project the added cumulants back to coefficients of powers of Z which I think can be done iteratively starting from smaller moments/cumulants and moving to higher moments/cumulants. At every step [$]c_0[$] would give us value of median where we could do expansion for drift function (and volatility function if we are not in bessel coordinates) and continue to use the expansion to solve for value of coefficients on next step.
The above holds for SDEs in Bessel coordinates but finding a solution would be more difficult in original coordinates and for dependent functions but I am very sure it can be done.
.
Friends my analytic calculations for solution of densities of SDEs by cumulant matching is going well. I want to explain some simple details of my suggested method.
Suppose initial SDE variable at start of simulation time step is given as
[$]X_0 \,=   \, \sum a_n Z^n + \, \sum b_n Z^n = \, \sum c_n Z^n [$] where [$]c_n \,=\, a_n \,+ \,b_n[$]

Initially I am working with fourth order so our initial SDE variable is given as

[$]X_0 \,= c_0\, + \, c_ 1 \, Z + \, c_ 2 \, Z^2 + \, c_ 3 \, Z^3 + \, c_ 4 \, Z^4[$]

I suppose that after a small time step the coefficients have changed by a small value each

[$]X_2 \,= (c_0\, + dc_0) + \,( c_ 1 +dc_1) \, Z + \,( c_ 2 +dc_2) \, Z^2 + \, ( c_ 3 +dc_3) \, Z + \, (c_ 4+dc_4) \, Z^4[$]

I write the cumulants before the simulation step [$]C_n(X_0)[$] and after the simulation step [$]C_n(X_2)[$] from the coefficients of Z in above two respective equation and match them with change in cumulant values due to addition of [$]X_1[$] as

[$]C_n(X_2) \, = \, C_n(X_0) \, + C_n(X_1) \,[$]

I suppose that [$]dc_2[$], [$]dc_3[$] and [$]dc_4[$] are very small for a small time step and their squares and higher powers can be neglected. I retain the square of leading variance coefficient [$]dc_1[$] that is change in coefficient of Z. It gives me system of equations of sort

[$]dc_0 + dc_2 + 3dc_4=C_1(X_1)[$]
[$]A_1 \, {(dc_1)}^2 \, + \, B_1 \, dc_1 \, + \, C_1 \, dc_2 \,  + \, D_1 \, dc_3 \,  + \, E_1 \, dc_4 \,= \, C_2(X_1)[$]
[$]A_2 \, {(dc_1)}^2 \, + \, B_2 \, dc_1 \, + \, C_2 \, dc_2 \,  + \, D_2 \, dc_3 \,  + \, E_2 \, dc_4 \,= \, C_3(X_1)[$]
[$]A_3 \, {(dc_1)}^2 \, + \, B_3 \, dc_1 \, + \, C_3 \, dc_2 \,  + \, D_3 \, dc_3 \,  + \, E_3 \, dc_4 \,= \, C_4(X_1)[$]

In above syatem of equations, [$]A_n, \, B_n, \, C_n, D_n , \, E_n [$] are coefficients that are already known(These are (polynomial-like) expressions in terms of [$]c_1, \, c_2, \, c_3, \, c_4[$] which are known at start of simulation step and their values can be calculated at start of simulation. The above system of equations can easily be solved in mathematica and this solution can be imported into matlab or C++.
Just to elaborate a bit further, I am copying the second of the above system of equations based on difference in second cumulant from mathematica as
[$] ({dc_1})^2 +(2 c_1 + 6 c_3) dc_1 + (4 c_2 + 24 c_4) dc_2 + (6 c_1 + 30 c_3) dc_3 + (24 c_2 + 192 c_4) dc_4= \, C_2(X_1)[$]
meaning [$]A_1=1[$], [$]B_1=(2 c_1 + 6 c_3) [$], [$]C_1=(4 c_2 + 24 c_4)[$] , [$]D_1\,=\,(6 c_1 + 30 c_3) [$] and [$]E_1=\,  (24 c_2 + 192 c_4)[$]

However we have to exogeneously specify the value of [$]dc_0[$] that is change in median. I could not find a simple analytic way to solve for median and I used leading order of fokker-planck equation solution for that and it is a very simple equation in leading coefficients. I will try to find some other easy way to solve for median and continue to think about it later. I have still not get the analytic solution to numbers but I hope to provide more details and post a working program in another two days.

.
Friends, I have been able to find a solution to simulating the probability density of diffusions by matching first four cumulants as I mentioned in the copied post. It works very well and is very stable as well. I modelled series solution of SDE plus its drift as one random variable and SDE volatility term as second random variable and was able to retrieve coefficients of SDE solution series expansion at next step by matching first four cumulants using additive property of cumulants. It works extremely well and matrix inversion in Newton root-finding is also stable. I will be sharing my program in another day. I have only solved the problem in bessel coordinates only but I expect that original coordinates would be just as easy. I took the series to fourth power of Z only and it can easily be extended to six powers (seven terms in the series) but it would be a lot more tedious to deal with first seven terms as cumulant expansions become a bit cumbersome but it is certainly doable. I hope to post my program tomorrow or on Monday morning.

Amin · January 2nd, 2022, 12:59 pm

Friends, I am posting some graphs from SDE density simulation program based on cumulant matching.
I will be copying the prototype matlab program in next post. I want to work on this problem for a few weeks and prepare a formal paper on it and submit it to Wilmott Journal. However, I will still be explaining most details in a post tomorrow. I will also upload my mathematica notebook used to calculate cumulants and derivatives of the cumulant equation that is used for Newton root finding of system of equations.
I want to warn friends that it is just a prototype quickly prepared over a few days and it is not the final program. It has to be improved in many aspects.
Current solution has many minor problems like.
1. There is some slack in calculation of drift.
2. I have just assumed that volatility term is perfectly normal which is true over very small periods of time and other volatility terms in powers of Z have to be added.
3. Median calculation has to be improved. It works well but can break down under stress when volatility is very high or in some other cases
4. Since volatility accuracy is only of order of sqrt(dt), step size is very small.
5. Mean can become too off for extreme parameters.

Here are the graphs. Parameters and all other information is on the title.

Amin · January 2nd, 2022, 2:00 pm

I am copying the matlab program in this post. I will come back with a detailed post explaining both analytics and my matlab program. I will also upload my mathamtica notebook used to calculate different cumulants and their derivatives.

Here is the main routine.
,

function [] = FPESeriesCoeffsFromCumulantsWilmott()

%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/16*1/2;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.
Tt=128*2*2;%128*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=.1/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=46*2;  % No of normal density subdivisions

NnMid=4.0;

x0=1.0;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.75;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=1.0;%.950;   %mean reversion parameter.
theta=1.0;%mean reversion target
sigma0=1.0;%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)*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=4;
SeriesOrder=4;

wnStart=1;
tic

for tt=1:Ts
   % t1=(tt-1)*dt;
   % t2=(tt)*dt;
    if(tt==1)
        %dt=ds(1);
        %[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder);
        [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),ExpnOrder);
        c0=wMu0dt+w0;
        c(1)=sigma0*sqrt(ds(tt))
        c(2:10)=0.0;
        w0=c0;
    
    elseif(tt==2) 

    %dt=ds(tt);    
        
    [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),8);
     [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);
            
            
    else
        
    %dt=ds(tt);    
    %below wMu0dta is drift and dwMu0dtdwa are its derivatives.
    %wMu1dt is its first derivative and dwMu1dtdw are derivatives of first
    %derivative.


     %[wMu0dta,dwMu0dtdwa,wMu1dt,dwMu1dtdw] = BesselDriftAndDerivatives08(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,10);
     
     [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),8);
     [b0,b] = CalculateDriftbCoeffs08(wMu0dt,dwMu0dtdw,c,gamma);
%b0
%b
%str=input('Look at drift values');
     
%Below bm0 and bm10 are used for claulations of median.    
     
      [wMu0dta,dwMu0dtdwa,wMu1dt,dwMu1dtdw] = BesselDriftAndDerivatives08Bounded(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),10);
      bm0=wMu0dta;
      bm10=wMu1dt;
      
      [NewMedian] = CalculateTotalChangeInMedian(c0,c,bm0,bm10,sigma0,ds(tt));
     
   %  da0=NewMedian-w0-bm0*dt*1.0-0*bm10*bm0*dt^2;
     
   %  da0=.5*sigma0^2*(2*c(2))/c(1)^2*dt;%+(.5*sigma0^2*(2*c(2))/c(1)^2).*( .5*sigma0^2*(6*c(3))/c(1)^3- sigma0^2*(2*c(2))^2/c(1)^4)*dt^2/2;
     da0=NewMedian-w0-b0;
    
     a0=NewMedian;
     
     %Initial random variable is determined by coefficients given in a
     %below.
     a(1:4)=c(1:4)+b(1:4);
    
    [da1,da2,da3,da4] = AdvanceSeriesCumulantSolution(da0,a,sigma0,ds(tt));

    c0=a0;
    c(1)=a(1)+da1;
    c(2)=a(2)+da2;
    c(3)=a(3)+da3;
    c(4)=a(4)+da4;
    
    
    w0=c0;
    
    end

end
wnStart=1;

w(1:Nn)=c0;
for nn=1:ExpnOrder
    w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^nn;
end
Flag=0;
for nn=Nn-1:-1:1
    if((w(nn)>w(nn+1))&&(Flag==0))
        wnStart=nn+1;
        Flag=1;
    end
end


c0
c


%str=input('Look at numbers')
yy(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
yy0=((1-gamma).*w0).^(1/(1-gamma));


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

.

.

function [F,dF] = CalculateCumulantsAndDerivatives(dC1,dC2,dC3,dC4,a1,a2,a3,a4,da0,da1,da2,da3,da4)

F1=da0+da2+3*da4-dC1;

dF1da1=0;
dF1da2=1;
dF1da3=0;
dF1da4=3;

F2=da1^2 + 2* da2^2 + (6* a1 + 30* a3)* da3 + 15 *da3^2 +  ...
 da1* (2 *a1 + 6* a3 + 6 *da3) + (24 *a2 + 192 *a4) *da4 + 96 *da4^2 + ...
 da2* (4 *a2 + 24 *a4 + 24 *da4)-dC2;
 
dF2da1= 2 *a1 + 6 *a3 + 2 *da1 + 6 *da3;
dF2da2= 4 *a2 + 24 *a4 + 4 *da2 + 24 *da4;
dF2da3=6 *a1 + 30 *a3 + 6 *da1 + 30 *da3;
dF2da4=24 *a2 + 192 *a4 + 24 *da2 + 192 *da4;


F3=8 *da2^3 + (36 *a1^2 + 216 *a2^2 + 576 *a1 *a3 + 2700 *a3^2 + 4608 *a2 *a4 + ...
    28512 *a4^2) *da4 + (2304 *a2 + 28512 *a4) *da4^2 + 9504 *da4^3 + ...
 da1^2 *(6 *a2 + 36 *a4 + 6 *da2 + 36 *da4) + ...
 da2^2* (24 *a2 + 216 *a4 + 216 *da4) + ...
 da3^2 *(270 *a2 + 2700 *a4 + 2700 *da4) + ...
 da3* (72 *a1 *a2 + 540 *a2 *a3 + 576 *a1 *a4 + ...
    5400 *a3 *a4 + (576 *a1 + 5400 *a3) *da4) + ...
 da2* (6 *a1^2 + 24 *a2^2 + 72 *a1 *a3 + 270 *a3^2 + 432 *a2 *a4 + ... 
    2304 *a4^2 + (72 *a1 + 540 *a3) *da3 + ...
    270 *da3^2 + (432 *a2 + 4608 *a4) *da4 + 2304 *da4^2) + ...
 da1 *(12 *a1 *a2 + 72 *a2 *a3 + 72 *a1 *a4 + 576 *a3 *a4 + ...
    da2 *(12 *a1 + 72 *a3 + 72 *da3) + (72 *a1 + 576 *a3) *da4 + ... 
    da3* (72 *a2 + 576 *a4 + 576 *da4))-dC3;
    
 dF3da1=12 *a1 *a2 + 72 *a2 *a3 + 72 *a1 *a4 + 576 *a3 *a4 + ...
 da2* (12 *a1 + 72 *a3 + 72 *da3) + (72 *a1 + 576 *a3) *da4 + ...
 2 *da1* (6 *a2 + 36 *a4 + 6 *da2 + 36 *da4) + ...
 da3* (72 *a2 + 576 *a4 + 576 *da4);
 
 
 dF3da2=6 *a1^2 + 24 *a2^2 + 72 *a1 *a3 + 270 *a3^2 + 432 *a2 *a4 + 2304 *a4^2 + ...
 6 *da1^2 + 24 *da2^2 + (72 *a1 + 540 *a3) *da3 + 270 *da3^2 + ...
 da1* (12 *a1 + 72 *a3 + 72 *da3) + (432 *a2 + 4608 *a4) *da4 + 2304 *da4^2 + ...
 2 *da2* (24 *a2 + 216 *a4 + 216 *da4);
 
 dF3da3=72 *a1 *a2 + 540 *a2 *a3 + 576 *a1 *a4 + 5400 *a3 *a4 + ...
 da2* (72 *a1 + 540 *a3 + 540 *da3) + (576 *a1 + 5400 *a3) *da4 + ...
 da1* (72 *a2 + 576 *a4 + 72 *da2 + 576 *da4) + ...
 2 *da3* (270 *a2 + 2700 *a4 + 2700 *da4);
 
 
 dF3da4=36 *a1^2 + 216 *a2^2 + 576 *a1 *a3 + 2700 *a3^2 + 4608 *a2 *a4 + ... 
 28512 *a4^2 + 36 *da1^2 + 216 *da2^2 + (576 *a1 + 5400 *a3) *da3 + ...
 2700 *da3^2 + da1* (72 *a1 + 576 *a3 + 576 *da3) + ...
 2* (2304 *a2 + 28512 *a4) *da4 + 28512 *da4^2 + ...
 da2 *(432 *a2 + 4608 *a4 + 4608 *da4);
 


F4=-3 *a1^4 - 12 *a1^2 *a2^2 - 12 *a2^4 - 36 *a1^3 *a3 - 72 *a1 *a2^2 *a3 - ...
 198 *a1^2 *a3^2 - 180 *a2^2 *a3^2 - 540 *a1 *a3^3 - 675 *a3^4 - ...
 144 *a1^2 *a2 *a4 - 288 *a2^3 *a4 - 864 *a1 *a2 *a3 *a4 - 2160 *a2 *a3^2 *a4 - ...
 576 *a1^2 *a4^2 - 2880 *a2^2 *a4^2 - 3456 *a1 *a3 *a4^2 - 8640 *a3^2 *a4^2 - ...
 13824 *a2 *a4^3 - 27648 *a4^4 + ...
 3 *(a1^2 - a2^2 + 114 *a4^2 + 15 *(a3^2 + 2 *a2 *a4) + ...
    3 *(a2^2 + 2 *a1 *a3 - 2 *a2 *a4 - 6 *a4^2))^2 + ...
 48 *da2^4 + (3240 *a1 + 38880 *a3) *da3^3 + 9720 *da3^4 + ...
 da1^3* (24 *a3 + 24 *da3) + (864 *a1^2 *a2 + 2304 *a2^3 + 18432 *a1 *a2 *a3 + ...
    108000 *a2 *a3^2 + 9216 *a1^2 *a4 + 92160 *a2^2 *a4 + 228096 *a1 *a3 *a4 + ...
    1537920 *a3^2 *a4 + 1368576 *a2 *a4^2 + ...
    7520256 *a4^3) *da4 + (4608 *a1^2 + 46080 *a2^2 + 114048 *a1 *a3 + ...
    768960 *a3^2 + 1368576 *a2 *a4 + 11280384 *a4^2) *da4^2 + (456192 *a2 + ...
    7520256 *a4) *da4^3 + 1880064 *da4^4 + ...
 da2^3* (192 *a2 + 2304 *a4 + 2304 *da4) + ...
 da2^2* (48 *a1^2 + 288 *a2^2 + 864 *a1 *a3 + 4320 *a3^2 + 6912 *a2 *a4 + ...
    46080 *a4^2 + (864 *a1 + 8640 *a3) *da3 + ...
    4320 *da3^2 + (6912 *a2 + 92160 *a4) *da4 + 46080 *da4^2) + ...
 da3^2* (432 *a1^2 + 4320 *a2^2 + 9720 *a1 *a3 + 58320 *a3^2 + ...
    108000 *a2 *a4 + 768960 *a4^2 + (108000 *a2 + 1537920 *a4) *da4 + ...
    768960 *da4^2) + ...
 da3* (24 *a1^3 + 864 *a1 *a2^2 + 864 *a1^2 *a3 + 8640 *a2^2 *a3 + ...
    9720 *a1 *a3^2 + 38880 *a3^3 + 18432 *a1 *a2 *a4 + 216000 *a2 *a3 *a4 + ...
    114048 *a1 *a4^2 + ...
    1537920 *a3 *a4^2 + (18432 *a1 *a2 + 216000 *a2 *a3 + 228096 *a1 *a4 + ...
       3075840 *a3 *a4) *da4 + (114048 *a1 + 1537920 *a3) *da4^2) + ...
 da1^2* (48 *a2^2 + 72 *a1 *a3 + 432 *a3^2 + 864 *a2 *a4 + 4608 *a4^2 + ...
    48 *da2^2 + (72 *a1 + 864 *a3) *da3 + ...
    432 *da3^2 + (864 *a2 + 9216 *a4) *da4 + 4608 *da4^2 + ...
    da2* (96 *a2 + 864 *a4 + 864 *da4)) + ...
 da2* (96 *a1^2 *a2 + 192 *a2^3 + 1728 *a1 *a2 *a3 + 8640 *a2 *a3^2 + ...
    864 *a1^2 *a4 + 6912 *a2^2 *a4 + 18432 *a1 *a3 *a4 + 108000 *a3^2 *a4 + ...
    92160 *a2 *a4^2 + ...
    456192 *a4^3 + (864 *a1^2 + 6912 *a2^2 + 18432 *a1 *a3 + 108000 *a3^2 + ...
       184320 *a2 *a4 + 1368576 *a4^2) *da4 + (92160 *a2 + ...
       1368576 *a4) *da4^2 + 456192 *da4^3 + ...
    da3^2* (8640 *a2 + 108000 *a4 + 108000 *da4) +... 
    da3* (1728 *a1 *a2 + 17280 *a2 *a3 + 18432 *a1 *a4 + ...
       216000 *a3 *a4 + (18432 *a1 + 216000 *a3) *da4)) + ...
 da1* (96 *a1 *a2^2 + 72 *a1^2 *a3 + 864 *a2^2 *a3 + 864 *a1 *a3^2 + ...
    3240 *a3^3 + 1728 *a1 *a2 *a4 + 18432 *a2 *a3 *a4 + 9216 *a1 *a4^2 + ...
    114048 *a3 *a4^2 + (864 *a1 + 9720 *a3) *da3^2 + 3240 *da3^3 + ...
    da2^2* (96 *a1 + 864 *a3 + 864 *da3) + (1728 *a1 *a2 + 18432 *a2 *a3 + ...
       18432 *a1 *a4 + 228096 *a3 *a4) *da4 + (9216 *a1 + 114048 *a3) *da4^2 + ...
     da3* (72 *a1^2 + 864 *a2^2 + 1728 *a1 *a3 + 9720 *a3^2 + 18432 *a2 *a4 + ...
       114048 *a4^2 + (18432 *a2 + 228096 *a4) *da4 + 114048 *da4^2) + ...
    da2* (192 *a1 *a2 + 1728 *a2 *a3 + 1728 *a1 *a4 + ...
       18432 *a3 *a4 + (1728 *a1 + 18432 *a3) *da4 + ...
       da3* (1728 *a2 + 18432 *a4 + 18432 *da4)))-dC4;
       
  dF4da1=     96 *a1 *a2^2 + 72 *a1^2 *a3 + 864 *a2^2 *a3 + 864 *a1 *a3^2 + 3240 *a3^3 + ...
 1728 *a1 *a2 *a4 + 18432 *a2 *a3 *a4 + 9216 *a1 *a4^2 + ...
 114048 *a3 *a4^2 + (864 *a1 + 9720 *a3) *da3^2 + 3240 *da3^3 + ...
 3 *da1^2 *(24 *a3 + 24 *da3) + ...
 da2^2* (96 *a1 + 864 *a3 + 864 *da3) + (1728 *a1 *a2 + 18432 *a2 *a3 + ...
    18432 *a1 *a4 + 228096 *a3 *a4) *da4 + (9216 *a1 + 114048 *a3) *da4^2 + ... 
 da3* (72 *a1^2 + 864 *a2^2 + 1728 *a1 *a3 + 9720 *a3^2 + 18432 *a2 *a4 + ...
    114048 *a4^2 + (18432 *a2 + 228096 *a4) *da4 + 114048 *da4^2) + ...
 2 *da1* (48 *a2^2 + 72 *a1 *a3 + 432 *a3^2 + 864 *a2 *a4 + 4608 *a4^2 + ... 
    48 *da2^2 + (72 *a1 + 864 *a3) *da3 + ...
    432 *da3^2 + (864 *a2 + 9216 *a4) *da4 + 4608 *da4^2 + ...
    da2 *(96 *a2 + 864 *a4 + 864 *da4)) + ...
 da2* (192 *a1 *a2 + 1728 *a2 *a3 + 1728 *a1 *a4 + ...
    18432 *a3 *a4 + (1728 *a1 + 18432 *a3) *da4 + ...
    da3* (1728 *a2 + 18432 *a4 + 18432 *da4));
    
    
 dF4da2=96 *a1^2 *a2 + 192 *a2^3 + 1728 *a1 *a2 *a3 + 8640 *a2 *a3^2 + 864 *a1^2 *a4 + ...
 6912 *a2^2 *a4 + 18432 *a1 *a3 *a4 + 108000 *a3^2 *a4 + 92160 *a2 *a4^2 + ...
 456192 *a4^3 + ...
 192 *da2^3 + (864 *a1^2 + 6912 *a2^2 + 18432 *a1 *a3 + 108000 *a3^2 + ...
    184320 *a2 *a4 + 1368576 *a4^2) *da4 + (92160 *a2 + ...
    1368576 *a4) *da4^2 + 456192 *da4^3 + ...
 da1^2* (96 *a2 + 864 *a4 + 96 *da2 + 864 *da4) + ...
 3 *da2^2* (192 *a2 + 2304 *a4 + 2304 *da4) + ...
 da3^2* (8640 *a2 + 108000 *a4 + 108000 *da4) + ...
 da3* (1728 *a1 *a2 + 17280 *a2 *a3 + 18432 *a1 *a4 + ...
    216000 *a3 *a4 + (18432 *a1 + 216000 *a3) *da4) + ...
 2 *da2 *(48 *a1^2 + 288 *a2^2 + 864 *a1 *a3 + 4320 *a3^2 + 6912 *a2 *a4 + ...
    46080 *a4^2 + (864 *a1 + 8640 *a3) *da3 + ...
    4320 *da3^2 + (6912 *a2 + 92160 *a4) *da4 + 46080 *da4^2) + ...
 da1* (192 *a1 *a2 + 1728 *a2 *a3 + 1728 *a1 *a4 + 18432 *a3 *a4 + ...
    2 *da2 *(96 *a1 + 864 *a3 + 864 *da3) + (1728 *a1 + 18432 *a3) *da4 + ...
    da3* (1728 *a2 + 18432 *a4 + 18432 *da4));
    
    
    
dF4da3=24 *a1^3 + 864 *a1 *a2^2 + 864 *a1^2 *a3 + 8640 *a2^2 *a3 + 9720 *a1 *a3^2 + ...
 38880 *a3^3 + 18432 *a1 *a2 *a4 + 216000 *a2 *a3 *a4 + 114048 *a1 *a4^2 + ...
 1537920 *a3 *a4^2 + 24 *da1^3 + 3 *(3240 *a1 + 38880 *a3) *da3^2 + ...
 38880 *da3^3 + da1^2 *(72 *a1 + 864 *a3 + 864 *da3) + ...
 da2^2* (864 *a1 + 8640 *a3 + 8640 *da3) + (18432 *a1 *a2 + 216000 *a2 *a3 + ...
    228096 *a1 *a4 + 3075840 *a3 *a4) *da4 + (114048 *a1 + ...
    1537920 *a3) *da4^2 + ...
 2 *da3* (432 *a1^2 + 4320 *a2^2 + 9720 *a1 *a3 + 58320 *a3^2 + ...
    108000 *a2 *a4 + 768960 *a4^2 + (108000 *a2 + 1537920 *a4) *da4 + ...
    768960 *da4^2) + ...
 da1* (72 *a1^2 + 864 *a2^2 + 1728 *a1 *a3 + 9720 *a3^2 + 18432 *a2 *a4 + ...
    114048 *a4^2 + 864 *da2^2 + 2 *(864 *a1 + 9720 *a3) *da3 + ...
    9720 *da3^2 + (18432 *a2 + 228096 *a4) *da4 + 114048 *da4^2 + ... 
    da2 *(1728 *a2 + 18432 *a4 + 18432 *da4)) + ...
 da2* (1728 *a1 *a2 + 17280 *a2 *a3 + 18432 *a1 *a4 + ...
    216000 *a3 *a4 + (18432 *a1 + 216000 *a3) *da4 + ...
    2 *da3* (8640 *a2 + 108000 *a4 + 108000 *da4));
    
    
 dF4da4=   864 *a1^2 *a2 + 2304 *a2^3 + 18432 *a1 *a2 *a3 + 108000 *a2 *a3^2 + ...
 9216 *a1^2 *a4 + 92160 *a2^2 *a4 + 228096 *a1 *a3 *a4 + 1537920 *a3^2 *a4 + ...
 1368576 *a2 *a4^2 + 7520256 *a4^3 + 2304 *da2^3 + ...
 2* (4608 *a1^2 + 46080 *a2^2 + 114048 *a1 *a3 + 768960 *a3^2 + ...
    1368576 *a2 *a4 + 11280384 *a4^2) *da4 + ...
 3 *(456192 *a2 + 7520256 *a4) *da4^2 + 7520256 *da4^3 + ...
 da1^2* (864 *a2 + 9216 *a4 + 864 *da2 + 9216 *da4) + ...
 da2^2* (6912 *a2 + 92160 *a4 + 92160 *da4) + ...
 da3^2* (108000 *a2 + 1537920 *a4 + 1537920 *da4) + ...
 da3* (18432 *a1 *a2 + 216000 *a2 *a3 + 228096 *a1 *a4 + 3075840 *a3 *a4 + ...
    2* (114048 *a1 + 1537920 *a3) *da4) + ...
 da2* (864 *a1^2 + 6912 *a2^2 + 18432 *a1 *a3 + 108000 *a3^2 + ...
    184320 *a2 *a4 + 1368576 *a4^2 + (18432 *a1 + 216000 *a3) *da3 + ...
    108000 *da3^2 + 2 *(92160 *a2 + 1368576 *a4) *da4 + 1368576 *da4^2) + ... 
 da1* (1728 *a1 *a2 + 18432 *a2 *a3 + 18432 *a1 *a4 + 228096 *a3 *a4 + ...
    da2* (1728 *a1 + 18432 *a3 + 18432 *da3) + ...
    2* (9216 *a1 + 114048 *a3) *da4 + ...
    da3 *(18432 *a2 + 228096 *a4 + 228096 *da4));
    
    
 F(1,1)=F1;
 F(2,1)=F2;
 F(3,1)=F3;
 F(4,1)=F4;
       
dF(1,1)=dF1da1;
dF(1,2)=dF1da2;
dF(1,3)=dF1da3;
dF(1,4)=dF1da4;

dF(2,1)=dF2da1;
dF(2,2)=dF2da2;
dF(2,3)=dF2da3;
dF(2,4)=dF2da4;

dF(3,1)=dF3da1;
dF(3,2)=dF3da2;
dF(3,3)=dF3da3;
dF(3,4)=dF3da4;

dF(4,1)=dF4da1;
dF(4,2)=dF4da2;
dF(4,3)=dF4da3;
dF(4,4)=dF4da4;



end

.

.

function [da1,da2,da3,da4] = AdvanceSeriesCumulantSolution(da0,a,sigma0,dt)
%da0=0;


dC1=0;
dC2=sigma0^2*dt;
dC3=0;
dC4=0;

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

da1=sqrt(a1^2+sigma0^2*dt)-a1;
da2=0;
da3=0;
da4=0;

da(1,1)=da1;
da(2,1)=da2;
da(3,1)=da3;
da(4,1)=da4;

[F,dF] = CalculateCumulantsAndDerivatives(dC1,dC2,dC3,dC4,a1,a2,a3,a4,da0,da1,da2,da3,da4);

nn=0;
while((abs(F(1,1))>.00000000001) || (abs(F(2,1))>.00000000001) || (abs(F(3,1))>.00000000001) || (abs(F(4,1))>.00000000001))
nn=nn+1;
da=da-dF\F;

%inv(dF)

da1=da(1,1);
da2=da(2,1);
da3=da(3,1);
da4=da(4,1);

[F,dF] = CalculateCumulantsAndDerivatives(dC1,dC2,dC3,dC4,a1,a2,a3,a4,da0,da1,da2,da3,da4);

end

%a(1)=a1+da1;
%a(2)=a2+da2;
%a(3)=a3+da3;
%a(4)=a4+da4;



end

.
.
.

function [NewMedian] = CalculateTotalChangeInMedian(a0,a,b0,b10,sigma0,dt)

NewMedian=a0+(b0+.5*sigma0^2*(2*a(2))/a(1)^2)*dt + ...
    (b0+.5*sigma0^2*(2*a(2))/a(1)^2)*(b10*0+.5*sigma0^2*1/a(1)^2 + .5*sigma0^2*(6*a(3))/a(1)^3- sigma0^2*(2*a(2))^2/a(1)^4)*dt^2/2;


end

.
.
.

function [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02(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;


YqCoeff(1)=mu1*dt;
YqCoeff(2)=mu2*dt;
YqCoeff(3)=-.5*sigma0^2*gamma*dt;

YqCoeff(4)=mu1^2*(beta1-gamma)*dt^2;
YqCoeff(5)=mu2^2*(beta2-gamma)*dt^2;
YqCoeff(6)=-.25*sigma0^4*gamma^2*(1-gamma)*dt^2;
YqCoeff(7)=(mu1*mu2*(beta1-gamma)+mu1*mu2*(beta2-gamma))*dt^2;
YqCoeff(8)=.5*mu1*sigma0^2*gamma*(1-gamma)*dt^2;
YqCoeff(9)=.5*mu2*sigma0^2*gamma*(1-gamma)*dt^2;


% YqCoeff(10)=YqCoeff0(1,1,4,1)*dt^3;
% YqCoeff(11)=YqCoeff0(1,2,3,1)*dt^3;
% YqCoeff(12)=YqCoeff0(2,1,3,1)*dt^3;
% YqCoeff(13)=YqCoeff0(1,3,2,1)*dt^3;
% YqCoeff(14)=YqCoeff0(2,2,2,1)*dt^3;
% YqCoeff(15)=YqCoeff0(3,1,2,1)*dt^3;
% YqCoeff(16)=YqCoeff0(1,4,1,1)*dt^3;
% YqCoeff(17)=YqCoeff0(2,3,1,1)*dt^3;
% YqCoeff(18)=YqCoeff0(3,2,1,1)*dt^3;
% YqCoeff(19)=YqCoeff0(4,1,1,1)*dt^3;

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;

% Fp(10)=Fp1(1,1,4,1);
% Fp(11)=Fp1(1,2,3,1);
% Fp(12)=Fp1(2,1,3,1);
% Fp(13)=Fp1(1,3,2,1);
% Fp(14)=Fp1(2,2,2,1);
% Fp(15)=Fp1(3,1,2,1);
% Fp(16)=Fp1(1,4,1,1);
% Fp(17)=Fp1(2,3,1,1);
% Fp(18)=Fp1(3,2,1,1);
% Fp(19)=Fp1(4,1,1,1);

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-gamma)/((1-gamma)*w0);
        else
        dwMu0dt(nn)=dwMu0dt(nn-1)*(Fp2(mm)-(nn-1))*(1-gamma)/((1-gamma)*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(9: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

.
.

Amin · January 3rd, 2022, 2:07 pm

Friends, there was a potential error in my calculations of objective function i.e the system of cumulant equations. I had taken one expectation out of order in calculations of cumulant equations. I re-did the calculations and found out that even with perfectly correct order of taking expectations, I get exactly the same result as in my program I posted yesterday. However, I am posting the corrected matlab function which is a more computationally intensive expression but it gives precisely the same results as my previous function. Here is the new function that has to replace the earlier function with similar name. You can remove the suffix 02 in the function below and make it work with the old program. One warning is that previous function is more efficient but gives exactly the same results as this new function. But I am posting corrected function for the sake of accuracy.
.

function [F,dF] = CalculateCumulantsAndDerivatives02(dC1,dC2,dC3,dC4,a1,a2,a3,a4,da0,da1,da2,da3,da4)

F1=da0+da2+3*da4-dC1;

dF1da1=0;
dF1da2=1;
dF1da3=0;
dF1da4=3;


F2=3* a2^2 + 6* a1* a3 + 15* a3^2 + 24* a2* a4 - 18* a4^2 - ... 
 15 *(a3^2 + 2* a2* a4) - ...
 3 *(a2^2 + 2* a1* a3 - 2* a2* a4 - 6* a4^2) + ... 
 da1^2 + 2* da2^2 + (6* a1 + 30* a3)* da3 + 15 *da3^2 +  ...
 da1* (2 *a1 + 6* a3 + 6 *da3) + (24 *a2 + 192 *a4) *da4 + 96 *da4^2 + ...
 da2* (4 *a2 + 24 *a4 + 24 *da4)-dC2;

 
dF2da1= 2 *a1 + 6 *a3 + 2 *da1 + 6 *da3;
dF2da2= 4 *a2 + 24 *a4 + 4 *da2 + 24 *da4;
dF2da3=6 *a1 + 30 *a3 + 6 *da1 + 30 *da3;
dF2da4=24 *a2 + 192 *a4 + 24 *da2 + 192 *da4;


    
    
F3=9* a1^2 *a2 + 6 *a2^3 + 72* a1* a2* a3 + 270* a2* a3^2 + 45* a1^2* a4 + ... 
 207* a2^2* a4 + 576* a1* a3* a4 + 2700* a3^2* a4 + 2304* a2* a4^2 - ... 
 864* a4^3 - ... 
 15* (a2^3 + 6* a1* a2* a3 - 3* a2* a3^2 + 3* a1^2* a4 - 6* a2^2 *a4 - ... 
    9* a3^2* a4 - 18* a2* a4^2) - 945* (3* a3^2* a4 + 3* a2* a4^2) - ... 
 105 *(3* a2* a3^2 + 3* a2^2* a4 + 6* a1* a3* a4 - 3 *a2* a4^2 - 9* a4^3) - ... 
 3* (3* a1^2* a2 - 3* a2^3 - 6* a1* a2* a3 - 6 *a2^2* a4 - 18* a1* a3 *a4 + ... 
    18* a2* a4^2 + 27* a4^3) + ... 
    8 *da2^3 + (36 *a1^2 + 216 *a2^2 + 576 *a1 *a3 + 2700 *a3^2 + 4608 *a2 *a4 + ...
    28512 *a4^2) *da4 + (2304 *a2 + 28512 *a4) *da4^2 + 9504 *da4^3 + ...
 da1^2 *(6 *a2 + 36 *a4 + 6 *da2 + 36 *da4) + ...
 da2^2* (24 *a2 + 216 *a4 + 216 *da4) + ...
 da3^2 *(270 *a2 + 2700 *a4 + 2700 *da4) + ...
 da3* (72 *a1 *a2 + 540 *a2 *a3 + 576 *a1 *a4 + ...
    5400 *a3 *a4 + (576 *a1 + 5400 *a3) *da4) + ...
 da2* (6 *a1^2 + 24 *a2^2 + 72 *a1 *a3 + 270 *a3^2 + 432 *a2 *a4 + ... 
    2304 *a4^2 + (72 *a1 + 540 *a3) *da3 + ...
    270 *da3^2 + (432 *a2 + 4608 *a4) *da4 + 2304 *da4^2) + ...
 da1 *(12 *a1 *a2 + 72 *a2 *a3 + 72 *a1 *a4 + 576 *a3 *a4 + ...
    da2 *(12 *a1 + 72 *a3 + 72 *da3) + (72 *a1 + 576 *a3) *da4 + ... 
    da3* (72 *a2 + 576 *a4 + 576 *da4))-dC3;
    
 dF3da1=12 *a1 *a2 + 72 *a2 *a3 + 72 *a1 *a4 + 576 *a3 *a4 + ...
 da2* (12 *a1 + 72 *a3 + 72 *da3) + (72 *a1 + 576 *a3) *da4 + ...
 2 *da1* (6 *a2 + 36 *a4 + 6 *da2 + 36 *da4) + ...
 da3* (72 *a2 + 576 *a4 + 576 *da4);
 
 
 dF3da2=6 *a1^2 + 24 *a2^2 + 72 *a1 *a3 + 270 *a3^2 + 432 *a2 *a4 + 2304 *a4^2 + ...
 6 *da1^2 + 24 *da2^2 + (72 *a1 + 540 *a3) *da3 + 270 *da3^2 + ...
 da1* (12 *a1 + 72 *a3 + 72 *da3) + (432 *a2 + 4608 *a4) *da4 + 2304 *da4^2 + ...
 2 *da2* (24 *a2 + 216 *a4 + 216 *da4);
 
 dF3da3=72 *a1 *a2 + 540 *a2 *a3 + 576 *a1 *a4 + 5400 *a3 *a4 + ...
 da2* (72 *a1 + 540 *a3 + 540 *da3) + (576 *a1 + 5400 *a3) *da4 + ...
 da1* (72 *a2 + 576 *a4 + 72 *da2 + 576 *da4) + ...
 2 *da3* (270 *a2 + 2700 *a4 + 2700 *da4);
 
 
 dF3da4=36 *a1^2 + 216 *a2^2 + 576 *a1 *a3 + 2700 *a3^2 + 4608 *a2 *a4 + ... 
 28512 *a4^2 + 36 *da1^2 + 216 *da2^2 + (576 *a1 + 5400 *a3) *da3 + ...
 2700 *da3^2 + da1* (72 *a1 + 576 *a3 + 576 *da3) + ...
 2* (2304 *a2 + 28512 *a4) *da4 + 28512 *da4^2 + ...
 da2 *(432 *a2 + 4608 *a4 + 4608 *da4);
 



F4=42* a1^2* a2^2 + 51* a2^4 + 24* a1^3 *a3 + 864* a1* a2^2* a3 + ... 
 432* a1^2* a3^2 + 4320* a2^2* a3^2 + 3240* a1* a3^3 + 9720* a3^4 + ... 
 828* a1^2* a2* a4 + 2328* a2^3* a4 + 18432* a1 *a2* a3* a4 + ... 
 108000* a2* a3^2* a4 + 4554* a1^2* a4^2 + 46134 *a2^2 *a4^2 + ... 
 114048 *a1 *a3 *a4^2 + 768960 *a3^2* a4^2 + 456192* a2* a4^3 - ... 
 147042* a4^4 - ...
 945* (6 *a2^2 *a3^2 + 4* a1* a3^3 + 4* a2^3 *a4 + 24 *a1 *a2 *a3 *a4 - ...
    12* a2* a3^2* a4 + 6* a1^2* a4^2 - 12* a2^2* a4^2 - 36* a3^2* a4^2 - ... 
    36* a2* a4^3) - 135135* (6* a3^2* a4^2 + 4 *a2* a4^3) - ... 
 15* (6* a1^2* a2^2 - 4 *a2^4 + 4* a1^3* a3 - 24* a1* a2^2 *a3 + 6 *a2^2* a3^2 - ... 
    12* a1^2* a2* a4 - 72* a1* a2* a3* a4 + 36* a2* a3^2* a4 - 36* a1^2* a4^2 + ... 
    72* a2^2* a4^2 + 54* a3^2* a4^2 + 108* a2* a4^3) - ... 
 3* (a1^4 - 12* a1^2* a2^2 + 6* a2^4 + 12* a1* a2^2* a3 - 36* a1^2* a2* a4 + ... 
    32* a2^3* a4 + 72* a1* a2* a3* a4 + 18* a2^2* a4^2 + 108* a1* a3* a4^2 - ... 
    108* a2* a4^3 - 108* a4^4) - ... 
 10395* (a3^4 + 12* a2* a3^2* a4 + 6* a2^2 *a4^2 + 12* a1* a3* a4^2 - ... 
    4* a2* a4^3 - 12* a4^4) - ... 
 105* (a2^4 + 12* a1* a2^2* a3 + 6* a1^2* a3^2 - 12* a2^2 *a3^2 + ... 
    12* a1^2* a2* a4 - 12* a2^3* a4 - 24* a1* a2* a3* a4 - 36* a2* a3^2* a4 - ... 
    30* a2^2 *a4^2 - 72* a1* a3* a4^2 + 36* a2* a4^3 + 54* a4^4) + ... 
 3* (a1^2 - a2^2 + 114* a4^2 + 15* (a3^2 + 2* a2* a4) + ... 
    3* (a2^2 + 2* a1 *a3 - 2* a2* a4 - 6* a4^2))^2 + ... 
 48* da2^4 + (3240* a1 + 38880* a3)* da3^3 + 9720* da3^4 + ... 
 da1^3* (24* a3 + 24* da3) + (864* a1^2* a2 + 2304 *a2^3 + 18432* a1* a2* a3 + ...
    108000* a2* a3^2 + 9216* a1^2* a4 + 92160 *a2^2* a4 + 228096* a1* a3 *a4 + ...
    1537920* a3^2* a4 + 1368576* a2* a4^2 + ... 
    7520256* a4^3)* da4 + (4608* a1^2 + 46080* a2^2 + 114048* a1* a3 + ... 
    768960* a3^2 + 1368576* a2* a4 + 11280384* a4^2)* da4^2 + (456192* a2 + ... 
    7520256* a4)* da4^3 + 1880064* da4^4 + ... 
 da2^3 *(192* a2 + 2304 *a4 + 2304* da4) + ...
 da2^2 *(48 *a1^2 + 288* a2^2 + 864* a1* a3 + 4320* a3^2 + 6912 *a2* a4 + ... 
    46080 *a4^2 + (864* a1 + 8640* a3)* da3 + ... 
    4320* da3^2 + (6912* a2 + 92160 *a4)* da4 + 46080* da4^2) +  ...
 da3^2 *(432* a1^2 + 4320 *a2^2 + 9720* a1* a3 + 58320* a3^2 + ...
    108000 *a2 *a4 + 768960 *a4^2 + (108000* a2 + 1537920* a4)* da4 + ... 
    768960* da4^2) + ... 
 da3* (24* a1^3 + 864 *a1* a2^2 + 864* a1^2* a3 + 8640 *a2^2 *a3 + ...
    9720* a1* a3^2 + 38880* a3^3 + 18432* a1* a2* a4 + 216000* a2* a3* a4 + ... 
    114048* a1 *a4^2 + ... 
    1537920* a3* a4^2 + (18432* a1* a2 + 216000 *a2* a3 + 228096* a1* a4 + ... 
       3075840* a3* a4)* da4 + (114048* a1 + 1537920* a3)* da4^2) + ... 
 da1^2* (48* a2^2 + 72* a1* a3 + 432 *a3^2 + 864 *a2* a4 + 4608* a4^2 + ... 
    48* da2^2 + (72* a1 + 864* a3)* da3 + ... 
    432* da3^2 + (864* a2 + 9216 *a4)* da4 + 4608* da4^2 + ... 
    da2* (96* a2 + 864* a4 + 864* da4)) + ... 
 da2 *(96 *a1^2* a2 + 192* a2^3 + 1728* a1* a2* a3 + 8640* a2* a3^2 + ... 
    864* a1^2* a4 + 6912* a2^2* a4 + 18432* a1* a3* a4 + 108000* a3^2* a4 + ... 
    92160* a2* a4^2 + ... 
    456192* a4^3 + (864* a1^2 + 6912* a2^2 + 18432 *a1* a3 + 108000* a3^2 + ... 
       184320* a2* a4 + 1368576* a4^2)* da4 + (92160 *a2 + ...
       1368576* a4)* da4^2 + 456192* da4^3 + ... 
    da3^2 *(8640* a2 + 108000* a4 + 108000* da4) + ... 
    da3* (1728* a1* a2 + 17280* a2* a3 + 18432* a1* a4 + ... 
       216000* a3* a4 + (18432* a1 + 216000 *a3)* da4)) + ... 
 da1* (96* a1* a2^2 + 72* a1^2* a3 + 864* a2^2 *a3 + 864* a1* a3^2 + ... 
    3240* a3^3 + 1728* a1* a2* a4 + 18432* a2* a3* a4 + 9216* a1* a4^2 + ... 
    114048* a3 *a4^2 + (864* a1 + 9720* a3)* da3^2 + 3240* da3^3 + ... 
    da2^2* (96* a1 + 864* a3 + 864* da3) + (1728* a1 *a2 + 18432* a2 *a3 + ...
       18432* a1* a4 + 228096 *a3 *a4)* da4 + (9216 *a1 + 114048* a3)* da4^2 + ...
     da3* (72* a1^2 + 864* a2^2 + 1728* a1* a3 + 9720* a3^2 + 18432 *a2* a4 + ... 
       114048* a4^2 + (18432* a2 + 228096* a4)* da4 + 114048* da4^2) + ... 
    da2* (192* a1* a2 + 1728* a2* a3 + 1728* a1* a4 + ... 
       18432* a3* a4 + (1728* a1 + 18432* a3)* da4 + ... 
       da3* (1728* a2 + 18432* a4 + 18432* da4)))-dC4;

 
  dF4da1=     96 *a1 *a2^2 + 72 *a1^2 *a3 + 864 *a2^2 *a3 + 864 *a1 *a3^2 + 3240 *a3^3 + ...
 1728 *a1 *a2 *a4 + 18432 *a2 *a3 *a4 + 9216 *a1 *a4^2 + ...
 114048 *a3 *a4^2 + (864 *a1 + 9720 *a3) *da3^2 + 3240 *da3^3 + ...
 3 *da1^2 *(24 *a3 + 24 *da3) + ...
 da2^2* (96 *a1 + 864 *a3 + 864 *da3) + (1728 *a1 *a2 + 18432 *a2 *a3 + ...
    18432 *a1 *a4 + 228096 *a3 *a4) *da4 + (9216 *a1 + 114048 *a3) *da4^2 + ... 
 da3* (72 *a1^2 + 864 *a2^2 + 1728 *a1 *a3 + 9720 *a3^2 + 18432 *a2 *a4 + ...
    114048 *a4^2 + (18432 *a2 + 228096 *a4) *da4 + 114048 *da4^2) + ...
 2 *da1* (48 *a2^2 + 72 *a1 *a3 + 432 *a3^2 + 864 *a2 *a4 + 4608 *a4^2 + ... 
    48 *da2^2 + (72 *a1 + 864 *a3) *da3 + ...
    432 *da3^2 + (864 *a2 + 9216 *a4) *da4 + 4608 *da4^2 + ...
    da2 *(96 *a2 + 864 *a4 + 864 *da4)) + ...
 da2* (192 *a1 *a2 + 1728 *a2 *a3 + 1728 *a1 *a4 + ...
    18432 *a3 *a4 + (1728 *a1 + 18432 *a3) *da4 + ...
    da3* (1728 *a2 + 18432 *a4 + 18432 *da4));
    
    
 dF4da2=96 *a1^2 *a2 + 192 *a2^3 + 1728 *a1 *a2 *a3 + 8640 *a2 *a3^2 + 864 *a1^2 *a4 + ...
 6912 *a2^2 *a4 + 18432 *a1 *a3 *a4 + 108000 *a3^2 *a4 + 92160 *a2 *a4^2 + ...
 456192 *a4^3 + ...
 192 *da2^3 + (864 *a1^2 + 6912 *a2^2 + 18432 *a1 *a3 + 108000 *a3^2 + ...
    184320 *a2 *a4 + 1368576 *a4^2) *da4 + (92160 *a2 + ...
    1368576 *a4) *da4^2 + 456192 *da4^3 + ...
 da1^2* (96 *a2 + 864 *a4 + 96 *da2 + 864 *da4) + ...
 3 *da2^2* (192 *a2 + 2304 *a4 + 2304 *da4) + ...
 da3^2* (8640 *a2 + 108000 *a4 + 108000 *da4) + ...
 da3* (1728 *a1 *a2 + 17280 *a2 *a3 + 18432 *a1 *a4 + ...
    216000 *a3 *a4 + (18432 *a1 + 216000 *a3) *da4) + ...
 2 *da2 *(48 *a1^2 + 288 *a2^2 + 864 *a1 *a3 + 4320 *a3^2 + 6912 *a2 *a4 + ...
    46080 *a4^2 + (864 *a1 + 8640 *a3) *da3 + ...
    4320 *da3^2 + (6912 *a2 + 92160 *a4) *da4 + 46080 *da4^2) + ...
 da1* (192 *a1 *a2 + 1728 *a2 *a3 + 1728 *a1 *a4 + 18432 *a3 *a4 + ...
    2 *da2 *(96 *a1 + 864 *a3 + 864 *da3) + (1728 *a1 + 18432 *a3) *da4 + ...
    da3* (1728 *a2 + 18432 *a4 + 18432 *da4));
    
    
    
dF4da3=24 *a1^3 + 864 *a1 *a2^2 + 864 *a1^2 *a3 + 8640 *a2^2 *a3 + 9720 *a1 *a3^2 + ...
 38880 *a3^3 + 18432 *a1 *a2 *a4 + 216000 *a2 *a3 *a4 + 114048 *a1 *a4^2 + ...
 1537920 *a3 *a4^2 + 24 *da1^3 + 3 *(3240 *a1 + 38880 *a3) *da3^2 + ...
 38880 *da3^3 + da1^2 *(72 *a1 + 864 *a3 + 864 *da3) + ...
 da2^2* (864 *a1 + 8640 *a3 + 8640 *da3) + (18432 *a1 *a2 + 216000 *a2 *a3 + ...
    228096 *a1 *a4 + 3075840 *a3 *a4) *da4 + (114048 *a1 + ...
    1537920 *a3) *da4^2 + ...
 2 *da3* (432 *a1^2 + 4320 *a2^2 + 9720 *a1 *a3 + 58320 *a3^2 + ...
    108000 *a2 *a4 + 768960 *a4^2 + (108000 *a2 + 1537920 *a4) *da4 + ...
    768960 *da4^2) + ...
 da1* (72 *a1^2 + 864 *a2^2 + 1728 *a1 *a3 + 9720 *a3^2 + 18432 *a2 *a4 + ...
    114048 *a4^2 + 864 *da2^2 + 2 *(864 *a1 + 9720 *a3) *da3 + ...
    9720 *da3^2 + (18432 *a2 + 228096 *a4) *da4 + 114048 *da4^2 + ... 
    da2 *(1728 *a2 + 18432 *a4 + 18432 *da4)) + ...
 da2* (1728 *a1 *a2 + 17280 *a2 *a3 + 18432 *a1 *a4 + ...
    216000 *a3 *a4 + (18432 *a1 + 216000 *a3) *da4 + ...
    2 *da3* (8640 *a2 + 108000 *a4 + 108000 *da4));
    
    
 dF4da4=   864 *a1^2 *a2 + 2304 *a2^3 + 18432 *a1 *a2 *a3 + 108000 *a2 *a3^2 + ...
 9216 *a1^2 *a4 + 92160 *a2^2 *a4 + 228096 *a1 *a3 *a4 + 1537920 *a3^2 *a4 + ...
 1368576 *a2 *a4^2 + 7520256 *a4^3 + 2304 *da2^3 + ...
 2* (4608 *a1^2 + 46080 *a2^2 + 114048 *a1 *a3 + 768960 *a3^2 + ...
    1368576 *a2 *a4 + 11280384 *a4^2) *da4 + ...
 3 *(456192 *a2 + 7520256 *a4) *da4^2 + 7520256 *da4^3 + ...
 da1^2* (864 *a2 + 9216 *a4 + 864 *da2 + 9216 *da4) + ...
 da2^2* (6912 *a2 + 92160 *a4 + 92160 *da4) + ...
 da3^2* (108000 *a2 + 1537920 *a4 + 1537920 *da4) + ...
 da3* (18432 *a1 *a2 + 216000 *a2 *a3 + 228096 *a1 *a4 + 3075840 *a3 *a4 + ...
    2* (114048 *a1 + 1537920 *a3) *da4) + ...
 da2* (864 *a1^2 + 6912 *a2^2 + 18432 *a1 *a3 + 108000 *a3^2 + ...
    184320 *a2 *a4 + 1368576 *a4^2 + (18432 *a1 + 216000 *a3) *da3 + ...
    108000 *da3^2 + 2 *(92160 *a2 + 1368576 *a4) *da4 + 1368576 *da4^2) + ... 
 da1* (1728 *a1 *a2 + 18432 *a2 *a3 + 18432 *a1 *a4 + 228096 *a3 *a4 + ...
    da2* (1728 *a1 + 18432 *a3 + 18432 *da3) + ...
    2* (9216 *a1 + 114048 *a3) *da4 + ...
    da3 *(18432 *a2 + 228096 *a4 + 228096 *da4));
    
    
 F(1,1)=F1;
 F(2,1)=F2;
 F(3,1)=F3;
 F(4,1)=F4;
       
dF(1,1)=dF1da1;
dF(1,2)=dF1da2;
dF(1,3)=dF1da3;
dF(1,4)=dF1da4;

dF(2,1)=dF2da1;
dF(2,2)=dF2da2;
dF(2,3)=dF2da3;
dF(2,4)=dF2da4;

dF(3,1)=dF3da1;
dF(3,2)=dF3da2;
dF(3,3)=dF3da3;
dF(3,4)=dF3da4;

dF(4,1)=dF4da1;
dF(4,2)=dF4da2;
dF(4,3)=dF4da3;
dF(4,4)=dF4da4;



end

Amin · January 3rd, 2022, 2:26 pm

Friends, just to clarify again that even though I apparently made an error in order of subtraction and expectation, the extra terms in corrected expression evaluate to zero. I checked their expectation and found that these extra terms consistently result in a value smaller than 10^-20 so you can safely continue to use the older function.

Amin · January 3rd, 2022, 4:10 pm

Friends, I am attaching my mathematica notebook for calculations of cumulants and their derivatives with respect to "change in coefficients" of Z-series representation of SDE variable.

SeriesFromCumulants02.zip: (17.62 KiB) Downloaded 33 times

Here I will explain some simple commands with comments
.

W1 := a0 + a1 Z + a2 Z^2 + a3 Z^3 + a4 Z^4

EW1 := a0 + a2 + 3 a4

W2 := (a0 + da0) + (a1 + da1) Z + (a2 + da2) Z^2 + (a3 + da3) Z^3 + (a4 + da4) Z^4

EW2 := (a0 + da0) + (a2 + da2) + 3 (a4 + da4)

.
In the above quoted part of the notebook, W1 is series representation of the "SDE variable + drift" and W2 is series representation of SDE variable after Bessel noise term has been added in the SDE. EW1 and EW2 are their respective expectations or means that follow from taking expectation on powers of standard normal.
.

Collect[(W1 - EW1)^2, Z]

a2^2 + 6 a2 a4 + 9 a4^2 +
(-2 a1 a2 - 6 a1 a4) Z +
(a1^2 - 2 a2^2 - 6 a2 a4) Z^2 +
(2 a1 a2 - 2 a2 a3 - 6 a3 a4) Z^3 +
(a2^2 + 2 a1 a3 - 2 a2 a4 - 6 a4^2) Z^4 +
(2 a2 a3 + 2 a1 a4) Z^5 +
(a3^2 + 2 a2 a4) Z^6 +
2 a3 a4 Z^7 +
a4^2 Z^8

a2^2 + 6 a2 a4 + 9 a4^2 +
(-2 a1 a2 - 6 a1 a4) 0 +
(a1^2 - 2 a2^2 - 6 a2 a4) 1 +
(2 a1 a2 - 2 a2 a3 - 6 a3 a4) 0 +
(a2^2 + 2 a1 a3 - 2 a2 a4 - 6 a4^2) 3 +
(2 a2 a3 + 2 a1 a4) 0 +
(a3^2 + 2 a2 a4) 15 +
(2 a3 a4) 0 +
a4^2 105

a1^2 - a2^2 + 114 a4^2 + 15 (a3^2 + 2 a2 a4) + 3 (a2^2 + 2 a1 a3 - 2 a2 a4 - 6 a4^2)

Ep2W1nEW1 := Out[28]

.
In the above quoted part, I calculate the value of (W1 - EW1)^2 and collect the terms with various powers of Z. I copy the answer and replace powers of standard normal Z with their respective expectations. This is especially convenient to collect terms of each powers of Z for higher cumulants where we have cumbersome expressions to just replace the expectation of each power of Z just once. I then re-calculate the value of expected expression and assign it to 2nd central moment or 2nd cumulant.
The nomenclature in notebook is that
Ep2W1nEW1 means second central moment of W1.
Ep3W1nEW1 means third central moment of W1.
Ep4W1nEW1 means fourth central moment of W1. (used to calculate fourth cumulant)
Ep2W2nEW2 means second central moment of W2.
Ep3W2nEW2 means third central moment of W2.
Ep4W2nEW2 means fourth central moment of W2. (used to calculate fourth cumulant)

I am sure with this orientation friends would easily understand the mathematica notebook.
I would explain rest of the simple things with equations tomorrow.

Amin · January 4th, 2022, 3:38 pm

Friends, I am posting corrected version of the drift calculation function. There were several minor errors in this function and I have fixed them now. Please use this new function in your program.
I hope to work tonight and be able to add functionality for taking into account higher order volatility terms of order dt^1.5. I will be getting antipsychotic injection tomorrow and do not know if I would be able to work properly after that so I want to try to complete it tonight.

Here is the corrected function you have to use.

function [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives02(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*(1-gamma)^a0;
B=mu2*(1-gamma)^b0;
C=-.5*sigma0^2*gamma/(1-gamma);
dt2=dt^2/2;

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

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


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;


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).*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

Amin · January 4th, 2022, 7:38 pm

Friends, I have been able to improve the cumulant matching density simulation program by adding volatility to order dt^1.5. I can easily take a step size of 1/64 and retain accuracy. Tomorrow I will try to further increase the order in volatility to dt^2. I hope to post the new program tomorrow.

Amin · January 6th, 2022, 12:30 am

I have been able to add second hermite polynomial in volatility and I also increased order everywhere and with all of this addition, I can now easily take a step size of 1/16 with very good accuracy. I have still not tried step sizes larger than 1/16 but I suspect that for easy diffusions, it would easily be possible to take a step size of 1/8 or 1/4 as we can take in monte carlo. Even closer to zero, it is a much better program than what we have previously been able to do. I have also increased accuracy in drift to order dt^3 for larger step sizes.
I will post the new program later today.

Amin · January 6th, 2022, 10:36 am

Friends, mind control agencies are drugging water(with mind control chemicals) in Lahore city again. I had not mentioned it on this forum but they had actively been drugging specific things I like in my food for past few weeks. I go once out of my home and have a quick meal on the fly somewhere. And bring some water and some food (like bread and some spread/jam etc) that I can eat rest of the day when I would be hungry. But mind control agencies have targeted bread and many spreads and distributed them in the markets in surrounding areas. I had told friends that it was impossible to get good coffee that would not have mind control chemicals in it. Pakistan army generals are given tens of billions of rupees every year and they ensure that every coffee manufacturer or importer drugs the coffee with mind control chemicals before distributing it in the market. I switched to tea to stimulate myself and think about mathematics but mind control crooks have started to drug most brands of tea and all the tea manufactured in past two months is drugged with mind control chemicals. Ground water was not getting drugged until recently but mind control agencies have started drugging groundwater and water I got today was completely drugged with mind control chemicals.
I want to request good Americans and especially good Jews to please stop this activity of retarding talented Muslims by some jewish actors who are behind this.
My original post was very different but I took a deep sigh and deleted most of it but I want to request people to not test my patience after a brutal and inhuman persecution spread over more than two decades. I just cannot allow my inhuman rape to continue forever at the hand of racist American defense.

Amin · January 7th, 2022, 7:37 am

https://forum.wilmott.com/viewtopic.php?f=15&t=94796&start=705#p821204

I want to ask a lot of people this philosophical question. Do most of us have naiive expectations of goodness from our societies? Most of us probably believe in these expectations of goodness until the society is put to a test when we realize that our expectations of goodness are probably a mirage. I still believe when I was in Hong Kong for three months at the end of 1999. I went to the largest newspaper group's offices to tell them about my persecution by mind control weapons. I was so thoroughly naiive that the staff at the newspaper group will be horrified at the sheer human rights abuses but once I left the group after a long interview with them, it was my turn to become aghast when I realized that they had already prepared the guy who was supposed to interview me at the Newspaper office and preparations had already been made by them how to respond to my accusations.
I have over the years seen again and again how mind control agents use influential, rich and connected people again and again to further their agenda of my persecution and when these influential people enthusiastically argue that I must be persecuted, most good natured people become quiet just because they do not want to get into problems for someone else and are afraid to fall out of favor with the influential people for the somebody else's cause. I keep thinking again and again how our societies continue to do something that is totally wrong and evil but still most people remain quiet because somebody rich or influential wants them to not make a big deal about these wrongdoings and cruelties of our society.
Again and again there were times when I was aghast how many people I looked up to with a sense of great uprightness failed me and wanted me to be persecuted when I came across them. When I was in London in 2010, I wrote about my persecution on a lot of LinkedIn groups. Though most of the linkedin groups allowed my posts to remain active and people had a god time making fun of me, Amnesty international of UK removed my posts about my human rights abuses from their group and banned me from making any further posts. I was horrified and wrote several emails to request the group owner to restore my posts on Amnesty international group. Apparently the group owner was told that I was some kind of a staunch muslim and he kept ranting against my religion. Mind control agents knew how to deal with different people and told people different stories that suited the nature of people so as to make them feel antagonistic towards me.
Again and again, I noticed that the society and the very people who are supposed to be flag-bearers of society's goodness and her values allow the society to perpetrate the crimes when influential people ask them and convince them to allow the wrong and bad things and human rights abuses to happen. Some people want to curry favor with rich and influential and rest of the people do not want to get into any trouble. And evil continues to get bigger and bigger in the society when rich and influential people learn that they can easily shape the behavior of the society in wrong ways without any resistance and they continue to become bolder every time they cause somebody's persecution successfully.
I still recall when I was in Karachi and I wrote email letters to my undergrad university. I had the same perception of the university that most other people had that mine was a very good school. I told the professors how mind control agents were given me some chemicals in water that my saliva in mouth would become dry and it would be an extremely painful situation. With great naivette, I thought that my problem would become over after the professors at my prestigious undergrad school would make sure that my persecution stops. Alas, because of the effect of some influential people, my university started convincing people that my persecution must continue and increase.
I had so many such experiences when very gentle and learned people who I thoroughly respected made fun of me. I knew they had been told lies about me and I had no way of telling them that these lies were not true. I just had to bear the ridicule and become quiet. These mind control people had a very good idea of psyche of people and they know how to deal with different types of people. In front of good and gentle people, they would tell lies so as to appeal to the goodness of people and convince them that I must be stopped because I am some kind of evil.
There are an estimated ten thousand humans who are a target of constant mind control abuse.
So I would request all good people to please understand the true nature of this mind control abuse. You would always notice someone influential who would try to make sure that innocent people must be made a target of mind control abuse and try to appeal to your goodness that this mind control abuse is indeed a good thing. If good people continue to yield, evil people would continue to get bolder. I would wish that all good people take a stand against such evils in our society so that our innocent perceptions of the goodness of the society in fact becomes a truth. Everybody is good when there is no incentive for being bad or evil, but true goodness is that you remain good when put to test and you have to become good against the wind. If all of us could do that, our societies would become a far better place.

Breakthrough in the theory of stochastic differential equations and their simulation

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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