Breakthrough in the theory of stochastic differential equations and their simulation

Amin · June 7th, 2022, 1:46 pm

For the series solution of SV equation, we have viewtopic.php?f=4&t=99702&start=1200#p867718

From the derivation above, we find that constant CDF lines of bessel form of the SV Fokker-planck equation are given as
[$]\frac{dw}{dt}\, =\, f(w)\, =\, \mu(w) \, + \, .5 {\sigma}^2 Z \frac{\partial Z}{\partial w} - \, .5 {\sigma}^2 \frac{\partial^2 Z}{\partial w^2} [\frac{\partial Z}{\partial w}]^{-1}[$]

Our second order solution to differential equation of constant CDF lines of SV SDE is given as

[$]w(t_1) \, =\, w(t_0)+ f(w(t_0)) \, \Delta t+\, f'(w(t_0)) f(w(t_0)) \,\frac{ {\Delta t}^2}{2} \,[$]

Coming to the asset equation that has dependence both on q itself and stochastic volatility V. We write the evolution equation for asset along constant CDF lines as

[$]\frac{dq}{dt}\, =\, g(q,V)\, =\, \mu(q,V) \, + \, .5 {\sigma_x}^2 {V(t)}^{2\, \gamma_v} \Big[ Z_2 \frac{\partial Z_2}{\partial q} - \, \frac{\partial^2 Z_2}{\partial q^2} [\frac{\partial Z_2}{\partial q}]^{-1} \Big][$]

Our second order solution to differential equation of constant CDF lines of Asset SDE is given as

[$]q(t_1) \, =\, q(t_0)+ g(q(t_0),V(t_0)) \, \Delta t+\, g_q(q(t_0),V(t_0)) g(q(t_0),V(t_0)) \,\frac{ {\Delta t}^2}{2} \,+\, g_V(q(t_0),V(t_0)) h(V(t_0)) \,\frac{ {\Delta t}^2}{2} \,[$]

I want to tell friends that we have still not included last second order term that depends on derivative of volatility in our matlab program and this remains to be done.
In the above equation
[$]h(V(t))\,= \, \frac{dV}{dt}\,[$]
which is the derivative of volatility (in order to remain on constant CDF lines) in original stochastic volatility coordinates. We only calculated [$]\frac{dw}{dt}\,=f(w)\,[$] in Bessel coordinates. For a proper second order solution, we will have to calculate [$]h(V(t))\,= \, \frac{dV}{dt}\,[$] in original stochastic volatility coordinates and then use it in our analysis. I plan on doing that with second version of the program due in a few weeks.

For reference of friends
[$] \frac{\partial V}{\partial t} = \, h(v(t))\, =\mu(V) - .5 \, \frac{\partial \big[{\sigma}^2 {V(t)}^{2 \gamma} \big]}{\partial V}+ .5\,{\sigma}^2 {V(t)}^{2 \gamma} \Big[ Z \frac{\partial Z}{\partial V} - \, \frac{\partial^2 Z}{\partial V^2} [\frac{\partial Z}{\partial V}]^{-1} \Big][$]

But on a second thought, it seems that we do not need to calculate [$]h(V(t))\,= \, \frac{dV}{dt}\,[$] in original coordinates and we can replace it with
[$]h(V(t))\,= \, \frac{dw}{dt}\,\, \frac{dV}{dw}\, =f(w(t)) \, \frac{dV}{dw}\,[$], where [$]\, \frac{dV}{dw}\,[$] can be calculated from bessel transformation relationship. So our asset equation to second order accuracy would be

[$]q(t_1) \, =\, q(t_0)+ g(q(t_0),V(t_0)) \, \Delta t+\, g_q(q(t_0),V(t_0)) g(q(t_0),V(t_0)) \,\frac{ {\Delta t}^2}{2} \,+\, g_V(q(t_0),V(t_0)) f(w(t_0)) \frac{dV}{dw} \, \,\frac{ {\Delta t}^2}{2} \,[$]

Amin · June 8th, 2022, 10:15 am

In this post I want to share my analytics for stochastic dt-integral of mean-reverting SDEs. A typical mean reverting SDE is of the form.

[$]dV(t) \, = \kappa \, (\theta \, - \, V(t) \,) dt + \, \sigma \, {V(t)}^{\gamma} \, dZ(t)\,[$]

We want to find the value of [$]\int_0^t \, V(s) \,ds[$] and later find step values of that integral over each time period given as [$]\int_{t_0}^{t_1} \, V(s) \,ds[$] and [$]\int_{t_1}^{t_2} \, V(s) \,ds[$] and [$]\int_{t_2}^{t_3} \, V(s) \,ds[$] and so on.

We first solve for the SDE(for different times) like we have usual treatment of classic mean-reverting SDEs. This gives us

[$]V(t_1)\, = \exp(- \, \kappa t_1) V(t_0) + \theta \, (1-\exp(- \, \kappa t_1)) +\, \exp(- \, \kappa t_1) \int_{t_0}^{t_1} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,[$]

for [$]V(t_2)\,[$] we have

[$]V(t_2)\, = \exp(- \, \kappa t_2) V(t_0) + \theta \, (1-\exp(- \, \kappa t_2)) +\, \exp(- \, \kappa t_2) \int_{t_0}^{t_2} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,[$]
[$] = \exp(- \, \kappa t_2) V(t_0) + \theta \, (1-\exp(- \, \kappa t_2)) +\, \exp(- \, \kappa t_2) \int_{t_0}^{t_1} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,+\, \exp(- \, \kappa t_2) \int_{t_1}^{t_2} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,[$]
In above, We have split the integral into orthogonal parts as
[$] \int_{t_0}^{t_2} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,=\,\int_{t_0}^{t_1} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,+ \int_{t_1}^{t_2} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,[$]

Similarly at date [$]t_3[$] our solution to V(t) is given as
[$]V(t_3) = \exp(- \, \kappa t_3) V(t_0) + \theta \, (1-\exp(- \, \kappa t_3)) +\, \exp(- \, \kappa t_3) \int_{t_0}^{t_1} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,[$]
[$]+\, \exp(- \, \kappa t_3) \int_{t_1}^{t_2} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,+\, \exp(- \, \kappa t_3) \int_{t_2}^{t_3} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,[$]

Now trying to approximate the dt-integral over first three time steps, we get
[$]\int_{t_0}^{t_3} \, V(t) \, dt\, = \big[ \, V(t_1) \, + \, V(t_2) \,+\, V(t_3) \, \big] \Delta t[$]
[$]=\Big[ \big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] V(t_0) +\big[ (1-\exp(- \, \kappa t_1)) + (1-\exp(- \, \kappa t_2)) + (1-\exp(- \, \kappa t_3)) \big] \theta [$]
[$]+\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] \,\int_{t_0}^{t_1} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,[$]
[$]+\big[ (\exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] \,\int_{t_1}^{t_2} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,[$]
[$]+\big[ (\exp(- \, \kappa t_3) ) \big] \,\int_{t_2}^{t_3} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\, \Big] \, \Delta t[$]
We can represent the stochastic integrals in above equation in hermite polynomial basis since they are already in the form of Z-series (We just take first three hermite polynomials in our analysis for brevity) , this gives us
[$]\,\int_{t_0}^{t_1} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,=a_{11} H_1(t_1) \, +a_{12} H_2(t_1) +\,a_{13} H_3(t_1) \, + ...[$]
[$]\,\int_{t_1}^{t_2} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,=a_{21} H_1(t_2) \, +a_{22} H_2(t_2) +\,a_{23} H_3(t_2) \, + ...[$]
[$]\,\int_{t_2}^{t_3} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,=a_{31} H_1(t_3) \, +a_{32} H_2(t_3) +\,a_{33} H_3(t_3) \, + ...[$]
Inserting above hermite representation in the approximate dt-integral, we get
[$]\int_{t_0}^{t_3} \, V(t) \, dt\, = \big[ \, V(t_1) \, + \, V(t_2) \,+\, V(t_3) \, \big] \Delta t[$]
[$]=\Big[ \big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] V(t_0) +\big[ (1-\exp(- \, \kappa t_1)) + (1-\exp(- \, \kappa t_2)) + (1-\exp(- \, \kappa t_3)) \big] \theta [$]
[$]+\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] \,\big[a_{11} H_1(t_1) \, +a_{12} H_2(t_1) +\,a_{13} H_3(t_1) \, + ...\big]\,[$]
[$]+\big[ (\exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] \,\,\big[a_{21} H_1(t_2) \, +a_{22} H_2(t_2) +\,a_{23} H_3(t_2) \, + ...\big]\,[$]
[$]+\big[ (\exp(- \, \kappa t_3) ) \big] \,\,\big[a_{31} H_1(t_3) \, +a_{32} H_2(t_3) +\,a_{33} H_3(t_3) \, + ...\big]\, \Big] \, \Delta t[$]

This gives us variance of dt-integral for first three time steps as
[$]Var \Big[\int_{t_0}^{t_3} \, V(t) \, dt\, \Big][$]
[$]=\Big[{\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {a_{11}}^2 Var\big[H_1(t_1)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {a_{12}}^2 Var\big[H_2(t_1)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {a_{13}}^2 Var\big[H_3(t_1)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {a_{21}}^2 Var\big[H_1(t_2)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {a_{22}}^2 Var\big[H_2(t_2)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {a_{23}}^2 Var\big[H_3(t_2)\big][$]
[$]+{\big[ \exp(- \, \kappa t_3) \big] }^2 \, {a_{31}}^2 Var\big[H_1(t_3)\big][$]
[$]+{\big[ \exp(- \, \kappa t_3) \big] }^2 \, {a_{32}}^2 Var\big[H_2(t_3)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_3) \big] }^2 \, {a_{33}}^2 Var\big[H_3(t_3)\big]\, \Big] {\Delta t}^2[$]

Since we want to match variances on every step, we have to calculate the above variances on every time step. On second time step the variances would be

[$]Var \Big[\int_{t_0}^{t_2} \, V(t) \, dt\, \Big][$]
[$]=\Big[{\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) ) \big] }^2 \, {a_{11}}^2 Var\big[H_1(t_1)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) ) \big] }^2 \, {a_{12}}^2 Var\big[H_2(t_1)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) ) \big] }^2 \, {a_{13}}^2 Var\big[H_3(t_1)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_2) ) \big] }^2 \, {a_{21}}^2 Var\big[H_1(t_2)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_2) ) \big] }^2 \, {a_{22}}^2 Var\big[H_2(t_2)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_2) ) \big] }^2 \, {a_{23}}^2 Var\big[H_3(t_2)\big] \Big] {\Delta t}^2[$]

Similarly variance of stochastic integral for first time step is given as

[$]Var \Big[\int_{t_0}^{t_1} \, V(t) \, dt\, \Big][$]
[$]=\Big[{\big[ (\exp(- \, \kappa t_1) ) \big] }^2 \, {a_{11}}^2 Var\big[H_1(t_1)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_1) ) \big] }^2 \, {a_{12}}^2 Var\big[H_2(t_1)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_1) ) \big] }^2 \, {a_{13}}^2 Var\big[H_3(t_1)\big] \Big] {\Delta t}^2[$]

Now we want to find coefficients of hermite representation so that on each step, total variance of each hermite polynomial matches with the variance of stochastic dt-integral at that step.
On the first step of the integral from [$]t_0[$] to [$] t_1[$] , there are no calculations and
[$] \Big[\int_{t_0}^{t_1} \, V(t) \, dt\, \Big][$]
[$]=c_1(t_1) \, H_1\,+ \,c_2(t_1) \, H_2\, + \, c_3(t_1) \, H_3\, [$]
where
[$]c_1(t_1) \,= \exp(- \, \kappa t_1) \, a_{11} \, \Delta t[$]
[$]c_2(t_1) \,= \exp(- \, \kappa t_1) \, a_{12} \, \Delta t[$]
[$]c_3(t_1) \,= \exp(- \, \kappa t_1) \, a_{13} \, \Delta t[$]

On second step, we have for one step addition that hermite coefficients are chosen so that total variance of dt-integral after second step given our hermite coefficients on previous step match with total variance after second step. We match coefficients of first hermite polynomial with cumulative variance of first hermite polynomial in the equations below. For that we have
[$]{\big[ c_1(t_1) \, + \, c_1(t_2) \, \big] }^2 = \Big[ \, {\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) ) \big] }^2 \, {a_{11}}^2+{\big[ (\exp(- \, \kappa t_2) ) \big] }^2 \, {a_{21}}^2 \Big] \, {\Delta t}^2[$]
(The reason to add [$]c_1(t_1)[$] and [$]c_1(t_2)[$] linearly inside the square is that our hermite polynomials are perfectly correlated across time and there is no way to maintain orthogonality over time as in monte carlo.)
which gives us
[$]\, c_1(t_2) \, = \Big[\,\sqrt{ {\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) ) \big] }^2 \, {a_{11}}^2+{\big[ (\exp(- \, \kappa t_2) ) \big] }^2 \, {a_{21}}^2}\,- \exp(- \, \kappa t_1) \, a_{11} \Big] \, \Delta t[$]
Similarly value of step coefficient of second hermite would be
[$]\, c_2(t_2) \, = \Big[ \,\sqrt{ {\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) ) \big] }^2 \, {a_{12}}^2+{\big[ (\exp(- \, \kappa t_2) ) \big] }^2 \, {a_{22}}^2}\,- \exp(- \, \kappa t_1) \, a_{12}\Big] \, \Delta t [$]
The value of step coefficients of third hermite would be
[$]\, c_3(t_2) \, = \Big[ \,\sqrt{ {\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) ) \big] }^2 \, {a_{13}}^2+{\big[ (\exp(- \, \kappa t_2) ) \big] }^2 \, {a_{23}}^2}\,- \exp(- \, \kappa t_1) \, a_{13}\Big] \, \Delta t [$]
With above values our hermite representation of second step of the integral would be
[$] \Big[\int_{t_1}^{t_2} \, V(t) \, dt\, \Big][$]
[$]=c_1(t_2) \, H_1\,+ \,c_2(t_2) \, H_2\, + \, c_3(t_2) \, H_3\, [$]

Similarly hermite coefficient for third time step [$]\, c_1(t_3) \,[$] should be that when it is added to coefficients of first two steps, the total variance equals the variance of dt-integral after third step. We match cumulative variance of first hermite polynomial over first three time steps with coefficients of first hermite polynomial over three time steps.
[$]{\big[ c_1(t_1) \, + \, c_1(t_2) \,\, + \, c_1(t_3) \, \big] }^2 = \Big[\, {\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2)+ \exp(- \, \kappa t_3) ) \big] }^2 \, {a_{11}}^2[$]
[$]+{\big[ (\exp(- \, \kappa t_2) + \exp(- \, \kappa t_3)) \big] }^2 \, {a_{21}}^2 + \big[\exp(- \, \kappa t_3) \, \big] \, {a_{31}}^2 \Big] \, {\Delta t}^2[$]
Solving for [$]\, c_1(t_3) \,[$], we get one step value of the first hermite coefficient as
[$]\, c_1(t_3) \, = \Big[ \, \sqrt{\big[ {\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2)+ \exp(- \, \kappa t_3) ) \big] }^2 \, {a_{11}}^2+{\big[ (\exp(- \, \kappa t_2) + \exp(- \, \kappa t_3)) \big] }^2 \, {a_{21}}^2 + \big[\exp(- \, \kappa t_3) \, \big] \, {a_{31}}^2 \big] } \, [$]
[$]- \,\sqrt{ {\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) ) \big] }^2 \, {a_{11}}^2+{\big[ (\exp(- \, \kappa t_2) ) \big] }^2 \, {a_{21}}^2}\, \Big] \, \Delta t[$]
Similarly we can calculate value of step coefficients of second hermite as
[$]\, c_2(t_3) \, =\Big[ \, \sqrt{\big[ {\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2)+ \exp(- \, \kappa t_3) ) \big] }^2 \, {a_{12}}^2+{\big[ (\exp(- \, \kappa t_2) + \exp(- \, \kappa t_3)) \big] }^2 \, {a_{22}}^2 + \big[\exp(- \, \kappa t_3) \, \big] \, {a_{32}}^2 \big] } \, [$]
[$]- \,\sqrt{ {\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) ) \big] }^2 \, {a_{12}}^2+{\big[ (\exp(- \, \kappa t_2) ) \big] }^2 \, {a_{22}}^2}\, \Big] \, \Delta t[$]
and
[$]\, c_3(t_3) \, =\Big[\, \sqrt{\big[ {\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2)+ \exp(- \, \kappa t_3) ) \big] }^2 \, {a_{13}}^2+{\big[ (\exp(- \, \kappa t_2) + \exp(- \, \kappa t_3)) \big] }^2 \, {a_{23}}^2 + \big[\exp(- \, \kappa t_3) \, \big] \, {a_{33}}^2 \big] } \, [$]
[$]- \,\sqrt{ {\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) ) \big] }^2 \, {a_{13}}^2+{\big[ (\exp(- \, \kappa t_2) ) \big] }^2 \, {a_{23}}^2}\, \, \Big] \, \Delta t[$]
With above values our hermite representation of third step of the integral would be
[$] \Big[\int_{t_2}^{t_3} \, V(t) \, dt\, \Big][$]
[$]=c_1(t_3) \, H_1\,+ \,c_2(t_3) \, H_2\, + \, c_3(t_3) \, H_3\, [$]

It turns out that coefficients on each hermite polynomial for one particular step of the stochastic time integral would be square root of the total variance of that particular hermite polynomial at the end of the time interval minus square root of the total variance of that particular hermite at the start of the time interval for a given dt-integral.

This is how I did the calculation of value of stochastic time integral over each step in the program.
I still have to discuss how to calculate the values of orthogonal hermite coefficients on each time step like [$] a_{11} \, , a_{21} \, ,a_{31} [$] for first hermite polynomial and [$] a_{12} \, , a_{22} \, ,a_{32} [$] for second hermite polynomial and so on.

Amin · June 8th, 2022, 11:57 am

To clarify, if there is only one stochastic variable(or dt-integral) in the future, its value can easily be represented by adding variances through addition of squared hermite coefficients. But if dt-integral is represented with hermites and there are large number of dependent parts of the integral, then we have to solve it in the way I wrote in the previous post.

Amin · June 10th, 2022, 12:57 pm

Friends, very sorry for this delay. I had an antipsychotic injection on Tuesday and it was getting difficult to concentrate on work.
Now coming to the topic of this post which is about calculation of Z-series and hermite representation of stochastic integrals in the Stochastic volatility SDE.
As I mentioned in yesterday's post that we want to find a hermite representation of the following integrals
[$]\,\int_{t_0}^{t_1} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,=ah_1(t_1) H_1(t_1) \, +ah_2(t_1) H_2(t_1) +\,ah_3(t_1) H_3(t_1) \, + ...[$]
[$]\,\int_{t_1}^{t_2} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,=ah_1(t_2) H_1(t_2) \, +ah_2(t_2) H_2(t_2) +\,ah_3(t_2) H_3(t_2) \, + ...[$]
[$]\,\int_{t_2}^{t_3} \, \exp( \, \kappa t)\sigma {V(t)}^{\gamma} \, dZ(t)\,=ah_1(t_3) H_1(t_3) \, +ah_2(t_3) H_2(t_3) +\,ah_3(t_3) H_3(t_3) \, + ...[$]

What we are given are solution of the stochastic volatility SDE in bessel coordinates. This solution in bessel coordinates is given as [$] w(t_1) [$], [$] w(t_2) [$], [$]w(t_3) [$] and so on. These solution are in the form of Z-series and can be represented as

[$]w(t_1)\, = \, c_0(t_1) \, + \, c_1(t_1) \, Z \, + \, c_2(t_1) \, Z^2 \, + \, c_3(t_1) \, Z^3 \, + \, c_4(t_1) \, Z^4 + ...[$]

and similarly for [$] t_2 [$], [$]t_3 [$], .. [$]t_n[$]

First we convert the above Z-series of SDE in bessel coordinates to Z-series of the SDE in original coordinates. We had learnt earlier how to convert a Z-series of a variable to Z-series of its functions and that is what we will use here. here
[$]V(t) \, = {\big[(1-\, \gamma) \, w(t) \big]}^{(1/(1- \, \gamma))}[$]

Then we will get a Z-series representation of V(t) in original coordinates. This would be in the form for example as

[$]V(t)\, = \, b_0(t) \, + \, b_1(t) \, Z \, + \, b_2(t) \, Z^2 \, + \, b_3(t) \, Z^3 \, + \, b_4(t) \, Z^4 + ...[$]

We can further find an equivalent hermite polynomial representation of the above Z-series as

[$]V(t)\, = \, bh_0(t) \, + \, bh_1(t) \, Z \, + \, bh_2(t) \, (Z^2 \, -1)\, + \, bh_3(t) \, (Z^3\, -\, 3 Z) \, + \, bh_4(t) \, (Z^4\, -6 Z^2 \, + \, 3) + ...[$]

We would have the above representation of V(t) on every time stage on the time grid.

From the solution of mean-reversion form of the stochastic volatility SDE, we have

[$]V(t_n) \, = V(t_{n-1}) \, \exp(- \kappa (t_n \, - t_{n -1} )) \, + \theta \, \big[ 1 \, - \, \exp(- \kappa (t_n \, - t_{n -1} ))[$]
[$]+ \,  \exp(- \kappa t_n \,) \, \int_{t_{n-1}}^{t_n} \,  \exp( \kappa t \,) \, dZ(t) \, [$]

We see in above that theta term does not have a coefficient in Z and it is not affecting the noises. So to gain intuition we write the above equation as

[$]\,  \exp(- \kappa t_n \,) \, \int_{t_{n-1}}^{t_n} \,  \exp( \kappa t \,) \, dZ(t) \, =V(t_n) \, - V(t_{n-1}) \, \exp(- \kappa (t_n \, - t_{n -1} )) \, [$]
which means
[$] \int_{t_{n-1}}^{t_n} \,  \exp( \kappa t \,) \, dZ(t) \, =\,  \exp(\kappa t_n \,) \,V(t_n) \, -\, \exp( \kappa \, t_{n -1} ) \, V(t_{n-1}) [$]

Let us suppose, we represent
[$] \int_{t_{n-1}}^{t_n} \,  \exp( \kappa t \,) \, dZ(t) \, = \, ah_0(t_n) \, + \, ah_1(t_n) \, Z \, + \, ah_2(t_n) \, (Z^2 \, -1)\, + \, ah_3(t_n) \, (Z^3\, -\, 3 Z) \, + \, ah_4(t_n) \, (Z^4\, -6 Z^2 \, + \, 3) + ...[$]

Then
[$]ah_1(t_n) \, = \, \sqrt{ \,{\big[ \exp(\kappa t_n \,) \, bh_1(t_n)\, \big]}^2 \, - \, \,{\big[ \exp(\kappa t_{n-1} \,) \, bh_1(t_{n-1})\, \big]}^2} [$]

[$]ah_2(t_n) \, = \, \sqrt{ \,{\big[ \exp(\kappa t_n \,) \, bh_2(t_n)\, \big]}^2 \, - \, \,{\big[ \exp(\kappa t_{n-1} \,) \, bh_2(t_{n-1})\, \big]}^2} [$]

[$]ah_3(t_n) \, = \, \sqrt{ \,{\big[ \exp(\kappa t_n \,) \, bh_3(t_n)\, \big]}^2 \, - \, \,{\big[ \exp(\kappa t_{n-1} \,) \, bh_3(t_{n-1})\, \big]}^2} [$]

Please note that above three equations are variance solution version of the equation below(that we proved earlier)

[$] \int_{t_{n-1}}^{t_n} \,  \exp( \kappa t \,) \, dZ(t) \, =\,  \exp(\kappa t_n \,) \,V(t_n) \, -\, \exp( \kappa \, t_{n -1} ) \, V(t_{n-1}) [$]

After we have found a hermite representation of the above integrals on every time step of the simulation, we can easily calculate the variance of the dt-integral as we did in a previous post as
For example the variance of dt-integral over first three dates would be
[$]Var \Big[\int_{t_0}^{t_3} \, V(t) \, dt\, \Big][$]
[$]=\Big[{\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, Var \big[ \int_{t_{0}}^{t_1} \,  \exp( \kappa t \,) \, dZ(t) \,\big][$]
[$]+{\big[ (\exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, Var \big[ \int_{t_{1}}^{t_2} \,  \exp( \kappa t \,) \, dZ(t) \,\big][$]
[$]+{\big[ \exp(- \, \kappa t_3) \big] }^2 \,   Var \big[ \int_{t_{2}}^{t_3} \, \exp( \kappa t \,) \, dZ(t) \,\big] \Big] \, {\Delta t}^2[$]

which could further be written as

[$]Var \Big[\int_{t_0}^{t_3} \, V(t) \, dt\, \Big][$]
[$]=\Big[{\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {ah_1(t_1)}^2 Var\big[H_1(t_1)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {ah_2(t_1)}^2 Var\big[H_2(t_1)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {ah_3(t_1)}^2 Var\big[H_3(t_1)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {ah_1(t_2)}^2 Var\big[H_1(t_2)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {ah_2(t_2)}^2 Var\big[H_2(t_2)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {ah_3(t_2)}^2 Var\big[H_3(t_2)\big][$]
[$]+{\big[ \exp(- \, \kappa t_3) \big] }^2 \, {ah_1(t_3)}^2 Var\big[H_1(t_3)\big][$]
[$]+{\big[ \exp(- \, \kappa t_3) \big] }^2 \, {ah_2(t_3)}^2 Var\big[H_2(t_3)\big][$]
[$]+{\big[ (\exp(- \, \kappa t_3) \big] }^2 \, {ah_3(t_3)}^2 Var\big[H_3(t_3)\big]\, \Big] {\Delta t}^2[$]

Suppose hermite representation of the integral  [$]\int_{t_0}^{t_3} \, V(t) \, dt\,[$] is given as
[$]\int_{t_0}^{t_3} \, V(t) \, dt\, [$]
[$]=vh0\, + \, vh_1 \, Z \, + \, vh_2 \, (Z^2-1) \,\, + \, vh_3 \, (Z^3-3 Z) \,\, + \, vh_4 \, (Z^4-6Z^2 +3) \,[$]

Then
[$]vh_1 \, = \, \sqrt{{\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {ah_1(t_1)}^2 +{\big[ (\exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {ah_1(t_2)}^2 + {\big[ \exp(- \, \kappa t_3) \big] }^2 \, {ah_1(t_3)}^2 }[$]
[$]vh_2 \, = \, \sqrt{{\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {ah_2(t_1)}^2 +{\big[ (\exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {ah_2(t_2)}^2 + {\big[ \exp(- \, \kappa t_3) \big] }^2 \, {ah_2(t_3)}^2 }[$]
[$]vh_3 \, = \, \sqrt{{\big[ (\exp(- \, \kappa t_1) + \exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {ah_3(t_1)}^2 +{\big[ (\exp(- \, \kappa t_2) +\exp(- \, \kappa t_3) ) \big] }^2 \, {ah_3(t_2)}^2 + {\big[ \exp(- \, \kappa t_3) \big] }^2 \, {ah_3(t_3)}^2 }[$]
and so on
so when above coefficients have been calculated, we know the values of coefficients of the series below
[$]\int_{t_0}^{t_3} \, V(t) \, dt\, [$]
[$]=vh0\, + \, vh_1 \, Z \, + \, vh_2 \, (Z^2-1) \,\, + \, vh_3 \, (Z^3-3 Z) \,\, + \, vh_4 \, (Z^4-6Z^2 +3) \,[$]

The above series in form of hermite polynomials and can easily be converted into a simpler Z-series that can be used for graphing of density of dt-integral and other purposes etc.

Amin · June 15th, 2022, 9:49 am

Friends, today, I was able to come back to work that had stopped due to antipsychotic injection. I have made several changes to the program. I saw there was still some discrepancy between monte carlo and our analytic density. I was able to determine the cause of small mismatch between monte carlo and analytic density and I would be posting my new matlab program in a day or two. Briefly I was able to improve the fit to program by adding second order correction term due to stochastic volatility. Soon I will be sharing my new program .

Amin · June 15th, 2022, 2:04 pm

For the series solution of SV equation, we have viewtopic.php?f=4&t=99702&start=1200#p867718

From the derivation above, we find that constant CDF lines of bessel form of the SV Fokker-planck equation are given as
[$]\frac{dw}{dt}\, =\, f(w)\, =\, \mu(w) \, + \, .5 {\sigma}^2 Z \frac{\partial Z}{\partial w} - \, .5 {\sigma}^2 \frac{\partial^2 Z}{\partial w^2} [\frac{\partial Z}{\partial w}]^{-1}[$]

Our second order solution to differential equation of constant CDF lines of SV SDE is given as

[$]w(t_1) \, =\, w(t_0)+ f(w(t_0)) \, \Delta t+\, f'(w(t_0)) f(w(t_0)) \,\frac{ {\Delta t}^2}{2} \,[$]

Coming to the asset equation that has dependence both on q itself and stochastic volatility V. We write the evolution equation for asset along constant CDF lines as

[$]\frac{dq}{dt}\, =\, g(q,V)\, =\, \mu(q,V) \, + \, .5 {\sigma_x}^2 {V(t)}^{2\, \gamma_v} \Big[ Z_2 \frac{\partial Z_2}{\partial q} - \, \frac{\partial^2 Z_2}{\partial q^2} [\frac{\partial Z_2}{\partial q}]^{-1} \Big][$]

Our second order solution to differential equation of constant CDF lines of Asset SDE is given as

[$]q(t_1) \, =\, q(t_0)+ g(q(t_0),V(t_0)) \, \Delta t+\, g_q(q(t_0),V(t_0)) g(q(t_0),V(t_0)) \,\frac{ {\Delta t}^2}{2} \,+\, g_V(q(t_0),V(t_0)) h(V(t_0)) \,\frac{ {\Delta t}^2}{2} \,[$]

I want to tell friends that we have still not included last second order term that depends on derivative of volatility in our matlab program and this remains to be done.
In the above equation
[$]h(V(t))\,= \, \frac{dV}{dt}\,[$]
which is the derivative of volatility (in order to remain on constant CDF lines) in original stochastic volatility coordinates. We only calculated [$]\frac{dw}{dt}\,=f(w)\,[$] in Bessel coordinates. For a proper second order solution, we will have to calculate [$]h(V(t))\,= \, \frac{dV}{dt}\,[$] in original stochastic volatility coordinates and then use it in our analysis. I plan on doing that with second version of the program due in a few weeks.

For reference of friends
[$] \frac{\partial V}{\partial t} = \, h(v(t))\, =\mu(V) - .5 \, \frac{\partial \big[{\sigma}^2 {V(t)}^{2 \gamma} \big]}{\partial V}+ .5\,{\sigma}^2 {V(t)}^{2 \gamma} \Big[ Z \frac{\partial Z}{\partial V} - \, \frac{\partial^2 Z}{\partial V^2} [\frac{\partial Z}{\partial V}]^{-1} \Big][$]

But on a second thought, it seems that we do not need to calculate [$]h(V(t))\,= \, \frac{dV}{dt}\,[$] in original coordinates and we can replace it with
[$]h(V(t))\,= \, \frac{dw}{dt}\,\, \frac{dV}{dw}\, =f(w(t)) \, \frac{dV}{dw}\,[$], where [$]\, \frac{dV}{dw}\,[$] can be calculated from bessel transformation relationship. So our asset equation to second order accuracy would be

[$]q(t_1) \, =\, q(t_0)+ g(q(t_0),V(t_0)) \, \Delta t+\, g_q(q(t_0),V(t_0)) g(q(t_0),V(t_0)) \,\frac{ {\Delta t}^2}{2} \,+\, g_V(q(t_0),V(t_0)) f(w(t_0)) \frac{dV}{dw} \, \,\frac{ {\Delta t}^2}{2} \,[$]

.
.
[$]q(t_1) \, =\, q(t_0)+ g(q(t_0),V(t_0)) \, \Delta t+\, g_q(q(t_0),V(t_0)) g(q(t_0),V(t_0)) \,\frac{ {\Delta t}^2}{2} \,+\, g_V(q(t_0),V(t_0)) f(w(t_0)) \frac{dV}{dw} \, \,\frac{ {\Delta t}^2}{2} \,[$]
Friends I think that last second order term in the above formula is incorrect.
What is the correct term, I still do not know.
Try
[$]q(t_1) \, =\, q(t_0)+ g(q(t_0),V(t_0)) \, \Delta t+\, g_q(q(t_0),V(t_0)) g(q(t_0),V(t_0)) \,\frac{ {\Delta t}^2}{2} \,+\, g_V(q(t_0),V(t_0)) f(w(t_0)) \, \,\frac{ {\Delta t}^2}{2} \,[$]
or
[$]q(t_1) \, =\, q(t_0)+ g(q(t_0),V(t_0)) \, \Delta t+\, g_q(q(t_0),V(t_0)) g(q(t_0),V(t_0)) \,\frac{ {\Delta t}^2}{2} \,+\, g_V(q(t_0),V(t_0)) f(w(t_0)) \, \frac{dw}{dV} \, \,\frac{ {\Delta t}^2}{2} \,[$]

I think this term somewhere takes a change of probability derivative and gets multiplied by an extra [$]\, \frac{dw}{dV} \,[$]. I still do not know exactly why. But above alternatives are very close to exact.

Amin · June 17th, 2022, 6:49 am

Friends, I have been begging the scum in American army to end their racist practices but it seems that racist scumbags in American army consider it my weakness and want to somehow continue my persecution forcefully through the rest of my life. Though many pseudo-patriotic(Patriotism is true and intelligent love for your country and people of your country and please do not equate it for love for scum armed forces and their weaponry). Americans may not like my speaking the truth about scum American army, matter of the fact is that American army is a large group of racist scumbags. Basically most large armies if they are not kept properly in check become scum. Is it not interesting that many Americans love their army but think of Pakistani army as full of scum bags. Pakistanis love their army and call American army full of scumbags. No, all large armies, yours and mine are equally full of scumbags. Basically scumbags in American army suck the dick of rich and racist jews so as to gain jobs paying millions of dollars in their companies after their retirement. Friends believe it or not it is really all about Benjamins. When tens of millions of Benjamins are given by corrupt racist crooks in American army to scumbag generals in Pakistani army, they start retarding hundreds of talented people in my country and they love to do it because they can never make such a large number of Benjamins any other way. So this chain of scumbags starting from racist scumbags in American army to scumbags in Pakistani army(other countries army) all wanting big millions of Benjamins. destroys talent all over the world. Please thwart the scum in American army and ask them to respect human beings.

Amin · June 17th, 2022, 4:28 pm

Friends, I reduced the step size in my simulation by 1/4 i.e. from 1/64 to 1/256 and made some other minor changes and now the program works perfectly well. We need very small steps at start of simulation and we should be able to take large steps later but that would require a time dependent step size and would take some effort.
Here are some new graphs. I would be posting the program used to make these graphs in a few minutes.
The fit to monte carlo is very close to perfect.
.
.

Amin · June 17th, 2022, 4:45 pm

Here is the program used to create above graphs (I have rescaled the graphs though to show the main body)
.
.
.

function [] = FPESeriesCoeffsFromGaussiansWilmottNew02A_Lite_2DFPEA4AB02()

%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334 
%or skype ahsan.amin2999

%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/2/4;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.
Tt=64*4*3;%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

Ts=Tt;
ds(1:Tt)=dt;

T
sum(ds(1:Ts))
Ts

%Ts=4;
str=input('Check if time is right');

OrderA=4;  %
OrderM=4;  %

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

%dtM=dt;
%TtM=Tt;

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

NnMid=4.0;

v00=1.0;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.950;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=1.0;%.950;   %mean reversion parameter.
theta=1.0;%mean reversion target
sigma0=.9500;%Volatility value
NoOfCumulants=6;
SeriesOrder=NoOfCumulants-1;

%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=-1*kappa;

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=v00^(1-gamma)/(1-gamma);
%y0=x0;


x0=1.00;
gammaX=.5;
gammaV=.5;
sigmaX=.50;
thetaX=.50;
kappaX=0.0;
mu1X=thetaX*kappaX;
mu2X=-1*kappaX;
beta1X=0.0;
beta2X=1.0;
q0=x0^(1-gammaX)/(1-gammaX);
alpha1X=1-gammaX;



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


sigma11X(1:OrderA+1)=0;
mu11X(1:OrderA+1)=0;
mu22X(1:OrderA+1)=0;
sigma22X(1:OrderA+1)=0;

sigma11X(1)=1;
mu11X(1)=1;
mu22X(1)=1;
sigma22X(1)=1;

for k=1:(OrderA+1)
    if sigmaX~=0
        sigma11X(k)=sigmaX^(k-1);
    end
    if mu1X ~= 0
        mu11X(k)=mu1X^(k-1);
    end
    if mu2X ~= 0
        mu22X(k)=mu2X^(k-1);
    end
    if sigmaX~=0
        sigma22X(k)=sigmaX^(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.

Fp1X(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
%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;
YqCoeff0X(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)
YqCoeffX = ItoTaylorCoeffsNew(alpha1X,beta1X,beta2X,gammaX);%expand y^alpha1 where alpha1=(1-gamma)
YqCoeffX=YqCoeffX/(1-gammaX); 
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));
                
                Fp1X(l1,l2,l3,l4) = (alpha1X + (l1-1) * beta1X + (l2-1) * beta2X + (l3-1) * 2* gammaX + (l4-1) * gammaX ...
                    - (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);
                YqCoeff0X(l1,l2,l3,l4) =YqCoeffX(l1,l2,l3,l4).*mu11X(l1).*mu22X(l2).*sigma22X(l3).*sigma11X(l4);
            end
        end
    end
end


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

ExpnOrder=6;
SeriesOrder=6;
c(1:SeriesOrder)=0;
c(1)=w0;
wnStart=1;
yydt(wnStart:Nn)=0.0;
tic

for tt=1:Ts
    tt

    if(tt<=1)
        [wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
    cprev(1:6)=0;
    cprev(1)=w0;
    c(1)=wMu0dt+w0;
    c(2)=sigma0*sqrt(ds(tt));
    c(3:6)=0.0;
    w0prev=w0;
    w0=c(1);
        
    elseif(tt<=2) 

    [wMu0dt,dwMu0dtdw,Ratio] = BesselDriftAndDerivatives04A(w0,mu1,mu2,beta1,beta2,sigma0,gamma,ds(tt),8);
     [b] = CalculateDriftbCoeffs08AB(wMu0dt,dwMu0dtdw,c,SeriesOrder);
     
        if(tt==2)
            
            cprev2(1)=v00^(1-gamma)/(1-gamma);
            cprev2(2:6)=0;
            cprev(1:6)=c(1:6);
        end
        if(tt==3)
            cprev2(1:6)=cprev(1:6);
            cprev(1:6)=c(1:6);
        end
    
     c(1)=c(1)+b(1);
     c(2)=b(2)+sqrt(c(2)^2+sigma0^2*ds(tt));
     c(3:SeriesOrder)=b(3:SeriesOrder);
     w0prev=w0;
     w0=c(1);
    else
        
       [wMu0dta,dwMu0dtdwa,wMu1dt,dwMu1dtdw] = BesselDriftAndDerivatives08(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,ExpnOrder+2);
     
     
     %The program CalculateDriftbCoeffs04 takes value of drift and its
     %derivatives wrt w and also takes Z-series of w and converts them into
     %Z-series of drift. It is used to calculate Z-series of drift and Then
     %to calculate Z_series of first derivative of drift wrt w.
   %  [b0,b] = CalculateDriftbCoeffs0(wMu0dta,dwMu0dtdwa,a);
   %  [b10,b1] = CalculateDriftbCoeffs04(wMu1dt,dwMu1dtdw,a);
     
     [b] = CalculateDriftbCoeffs08AB(wMu0dta,dwMu0dtdwa,c,SeriesOrder);
     [b1] = CalculateDriftbCoeffs08AB(wMu1dt,dwMu1dtdw,c,SeriesOrder);
        
        
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%555
     
     
     %SeriesOrder=ExpnOrder;
d1(1:SeriesOrder)=0;

for nn=1:ExpnOrder-1
    d1(nn)=(nn)*c(nn+1);
end
d2(1:ExpnOrder)=0;
for nn=1:ExpnOrder-2
    d2(nn)=(nn)*(nn+1)*c(nn+2);
end
d3(1:ExpnOrder)=0;
for nn=1:ExpnOrder-3
    d3(nn)=(nn)*(nn+1)*(nn+2)*c(nn+3);
end

[d1Sq] =SeriesProductAB(d1,d1,SeriesOrder);
[d1Cube] =SeriesProductAB(d1Sq,d1,SeriesOrder);
[d1Quad] =SeriesProductAB(d1Cube,d1,SeriesOrder);
[d1Inv] = SeriesReciprocalAB(d1,SeriesOrder);
[d1SqInv] = SeriesReciprocalAB(d1Sq,SeriesOrder);
[d1CubeInv] = SeriesReciprocalAB(d1Cube,SeriesOrder);
[d1QuadInv] = SeriesReciprocalAB(d1Quad,SeriesOrder);
[d2Sq] =SeriesProductAB(d2,d2,SeriesOrder);

Z_d1Inv(1:SeriesOrder)=0;

Z_d1Inv(2:SeriesOrder)=d1Inv(1:SeriesOrder-1);


T_Odt(1:SeriesOrder)=b(1:SeriesOrder)+.5*sigma0^2*(Z_d1Inv(1:SeriesOrder)+SeriesProductAB(d2,d1SqInv,SeriesOrder));

T_Odt
b
b1
c(1)
c
%str=input('Look at Torder dt')


dc1=T_Odt;

[d2_d1CubeInv] =SeriesProductAB(d2,d1CubeInv,SeriesOrder);

Z_d2_d1CubeInv(1:SeriesOrder)=0;
Z_d2_d1CubeInv(2:SeriesOrder)=d2_d1CubeInv(1:SeriesOrder-1);

T_Odt2(1:SeriesOrder)=b1(1:SeriesOrder)+.5*sigma0^2*(d1SqInv(1:SeriesOrder)-Z_d2_d1CubeInv(1:SeriesOrder)+ ...
    SeriesProductAB(d3,d1CubeInv,SeriesOrder)-2*SeriesProductAB(d2Sq,d1QuadInv,SeriesOrder));

cprev2(1:SeriesOrder)=cprev(1:SeriesOrder);
cprev(1:SeriesOrder)=c(1:SeriesOrder);

c(1:SeriesOrder)=c(1:SeriesOrder)+T_Odt(1:SeriesOrder)*dt+SeriesProductAB(T_Odt,T_Odt2,SeriesOrder)*dt^2/2;

     
     
     
   
    w0prev2=w0prev;
    w0prev=w0;
    w0=c(1);
    end
    
     if(tt==1)
      %  cmid0=c0;
      %  cmid=c;
        
     else
      %   cmid0prev=cmid0;
      %   cmidprev=cmid;
     %  cmid0=(c0+c0prev)*.5;
     %   cmid=(c+cprev)*.5;  
        %cmid0=c0;
        %cmid=c;  
        
     end
    
 
     
      gammaV2=2*gammaV;
      [vw0,dvw0] = CalculatevgammaV2AndDerivatives(w0prev,gamma,gammaV2,SeriesOrder);
%      %[V2gamma0,V2gamma1] = CalculateDriftbCoeffs08A(vw0,dvw0,cmid,SeriesOrder);
      [V2gamma] = CalculateDriftbCoeffs08AB(vw0,dvw0,cprev,SeriesOrder);
      
      gammaV2=2*gammaV;
      [vw0,dvw0] = CalculatevgammaV2AndDerivatives(w0,gamma,gammaV2,SeriesOrder);
%      %[V2gamma0,V2gamma1] = CalculateDriftbCoeffs08A(vw0,dvw0,cmid,SeriesOrder);
      [V2gammaNew] = CalculateDriftbCoeffs08AB(vw0,dvw0,c,SeriesOrder);
      
     % V2gamma(1)=V2gamma0;
     % V2gamma(2:6)=V2gamma1(1:5);
%      
     
     
    if(tt==1)

        cq(1:6,1:6)=0;    
        cq(1,1)=x0^(1-gammaX)/(1-gammaX);
        
        
    end
    %    cprev=c;
    
    q0=cq(1,1);

    if(tt<=1)
     
%     %gammaV2=gammaV;
%     [vw0,dvw0] = CalculatevgammaV2AndDerivatives(w0prev,gamma,gamma,SeriesOrder);
%     [bv] = CalculateDriftbCoeffs08AB(vw0,dvw0,cprev,SeriesOrder);
        
     gammaV2=2*gammaV;
     [vw0,dvw0] = CalculatevgammaV2AndDerivatives(w0prev,gamma,gammaV2,SeriesOrder);
     [ch] = CalculateDriftbCoeffs08AB(vw0,dvw0,cprev,SeriesOrder);
     
     V2gammaSum(1:5)=ch(2:6)*dt;
     
     [ch] = ConvertZCoeffsToHCoeffsAB(ch,SeriesOrder);    
    
 %    [ch_g] = ConvertZCoeffsToHCoeffsAB(bv,SeriesOrder);    
    
     
     if(tt==1)
     %cH0(tt)=ch(1);
     cH0(tt)=theta*dt+(v00-theta)/kappa*(exp(-kappa*(tt-1)*dt)-(exp(-kappa*tt*dt)));;
     
     cH(tt,1:5)=exp(kappa*dt).*ch(2:6);
     
     %cVdt(1)=ch(1);
     cVdt(1)=cH0(tt);
     cVdt(2:6)=exp(-1*kappa*dt).*ch(2:6);
     
     end
    else
        
    if(tt>=2)    
        
     [cH0,cH,cVdt] = CalculateVPathStepIntegral04(cprev,c,v00,kappa,theta,gamma,dt,tt,cH0,cH,SeriesOrder);
   
 %   cVdt
 %   tt
 %   str=input('Look at cVdt');
  
  
    end
    end
    
    
    

    DoubleExpnOrder=ExpnOrder*2;
    
    [qMu0dt,dqMu0dtdq,qMu1dt,dqMu1dtdq] = BesselDriftAndDerivatives08A0(q0,mu1X,mu2X,beta1X,beta2X,gammaX,DoubleExpnOrder);
    [q2Mu0dt,dq2Mu0dtdq,q2Mu1dt,dq2Mu1dtdq] = BesselDriftAndDerivatives08Aq(q0,sigmaX,gammaX,DoubleExpnOrder);
 
    [Muq0] = CalculateDriftbCoeffs08A2Dim02(qMu0dt,dqMu0dtdq,cq,SeriesOrder);
    [Muq2] = CalculateDriftbCoeffs08A2Dim02(q2Mu0dt,dq2Mu0dtdq,cq,SeriesOrder);
    [DMuq0] = CalculateDriftbCoeffs08A2Dim02(qMu1dt,dqMu1dtdq,cq,SeriesOrder);
    [DMuq2] = CalculateDriftbCoeffs08A2Dim02(q2Mu1dt,dq2Mu1dtdq,cq,SeriesOrder);
    
    DvMuq2=Muq2*dt;
 
    
    
    if(tt<=1)
    [Muq2] =SeriesProduct2D1D(Muq2,V2gamma)*dt;
    [DMuq2] =SeriesProduct2D1D(DMuq2,V2gamma)*dt;
    %[Tvol] =SeriesProduct2D1D(Tvol,V2gamma)*dt;
    else
    [Muq2] =SeriesProduct2D1D(Muq2,cVdt);
    [DMuq2] =SeriesProduct2D1D(DMuq2,cVdt);
    cVdt
    
    end
    
    
   % tt
   % str=input('Look at numbers');
    
    
    
    Muq=Muq0*dt+Muq2;%+Muq01*dt/2;%-tt*Muq02*dt/2+tt*Muq01*dt/2;
    DMuq=DMuq0*dt+DMuq2;
    
    if(tt==1)
        cq=cq+Muq+SeriesProduct2D(Muq,DMuq)*1/2;
        cq(2,1)=cq(2,1)+sigmaX*v00^gammaV*sqrt(dt);
    elseif(tt<=1)
%        cq
        [D1cq] = Series2DZ1Derivative(cq);
%        D1cq
        [D1cqInv] = Series2DReciprocal02(D1cq);
%        D1cqInv
        [D1cqInvZ] = MultiplySeries2DWithZ1(D1cqInv);
%        D1cqInvZ
        fq=Muq+.5*sigmaX^2*SeriesProduct2D1D(D1cqInvZ,V2gamma)*dt;
%        fq*dt
        Dfq=DMuq;
        cq=cq+fq+SeriesProduct2D(fq,Dfq)*dt^0/2;
        
        
        SeriesProduct2D1D(D1cqInvZ,V2gamma)
    else
        [D1cq] = Series2DZ1Derivative(cq);
        [D1cqSq] =SeriesProduct2D(D1cq,D1cq);
        %[d1Quad] =SeriesProductAB(d1Cube,d1,SeriesOrder);
        [D1cqInv] = Series2DReciprocal02(D1cq);
        [D1cqInvZ] = MultiplySeries2DWithZ1(D1cqInv);
        [D2cq] = Series2DZ1Derivative(D1cq);
        [D1cqInvSqr]=SeriesProduct2D(D1cqInv,D1cqInv);
        [D1cqInvCub]=SeriesProduct2D(D1cqInvSqr,D1cqInv);
        [D2cqD1cqInvSqr]=SeriesProduct2D(D2cq,D1cqInvSqr);
        [D2cqD1cqInvCub]=SeriesProduct2D(D2cq,D1cqInvCub);
        
        
        D2cqT2=D1cqInvSqr-MultiplySeries2DWithZ1(D2cqD1cqInvCub);%+SeriesProduct2D(D3cq,D1cqInvCub);
        Dfq=DMuq+.5*sigmaX^2*SeriesProduct2D1D(D2cqT2,cVdt);
        %Dfq=DMuq+.5*sigmaX^2*SeriesProduct2D(D2cqT2,Tvol);
        
        D2cqT1=D1cqInvZ+D2cqD1cqInvSqr;
        
        
        Dvfq=DvMuq2+.5*sigmaX^2*D2cqT1*dt;

        fq=Muq+.5*sigmaX^2*SeriesProduct2D1D(D2cqT1,cVdt);    
        %fq=Muq+.5*sigmaX^2*SeriesProduct2D(D2cqT1,Tvol);    
 
        if(tt>3)
        Dvfqdc1=SeriesProduct2D1D(Dvfq,dc1)*dt;
        end
        %Dvfqdc1=SeriesProduct2D1D(Dvfqdc1,bv)
        if(tt>3)
        cq=cq+(fq+SeriesProduct2D(fq,Dfq)*dt^0/2+sqrt(1).*Dvfqdc1.*dt^0/2);
        V2gammaNew
        else
        cq=cq+fq;    
        %cq=cq+(fq+SeriesProduct2D(fq,Dfq)*dt^0/2);%+sqrt(1).*Dvfqdc1.*dt^0/2);
        end
            
    end
 
 


end

Z1=Z;


[b] = EvaluateSeriesForZ2(cq,Z1,Nn);

Mm=1201;
MmMid=601;
dMm=.05/2;
Q(MmMid+1:Mm)=cq(1,1)+cq(2,1).*dMm.*(1:MmMid-1);
Q(MmMid)=cq(1,1);
%Q(1:70)=cq(1,1)-cq(2,1).*.25.*(70:-1:1);
for mm=1:MmMid-1

    Q(MmMid-mm)=cq(1,1)-cq(2,1).*dMm.*mm;
end
Q(Q<0)=0;

Qa=Q-cq(2,1).*dMm/2;
Qb=Q+cq(2,1).*dMm/2;

Qa(Qa<0.000000001)=0;





Z1Prob(1:Mm)=0;
for nn=1:Nn
    bq(1:6)=b(1:6,nn);

    
    [Z0Prob] = CalculateProbabilitySeriesInv(Qa,Qb,bq);
    
%    Z0Prob
    
%str=input('Look at numbers-tttt')
    
    
    Z1Prob(1:Mm)=Z1Prob(1:Mm)+Z0Prob(1:Mm).*ZProb(nn);
    
end

pQ=Z1Prob/(cq(2,1).*dMm);
    
xx(1:Mm)=((1-gammaX).*Q(1:Mm)).^(1/(1-gammaX));

xx(xx<.00000001)=0;

pxx(1:Mm)=pQ(1:Mm).*xx(1:Mm).^(-gammaX);
    
sum(Z1Prob(1:Mm))    
Mean=sum(xx(1:Mm).*Z1Prob(1:Mm));
Mean
  w(1:Nn)=c(1);
  for nn=2:SeriesOrder-1
      w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^(nn-1);
  end
 
  
  w
  
  str=input('Look at w---2');
  
%  %str=input('Look at numbers')
  yy0(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
 yy=yy0;
% xx=xx0;
Dfyy(wnStart:Nn)=0;
%Dfyy1(wnStart:Nn)=0;
%Dfyy2(wnStart:Nn)=0;
%Dfxx(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    
    Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%    Dfyy1(nn) = (yy1(nn + 1) - yy1(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%    Dfyy2(nn) = (yy2(nn + 1) - yy2(nn - 1))/(Z(nn + 1) - Z(nn - 1));
 %   Dfxx(nn) = (xx(nn + 1) - xx(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    
    %Change of variable derivative for densities
end
%Dfyy(Nn)=Dfyy(Nn-1);
%Dfyy(1)=Dfyy(2);
%Dfyydt(Nn)=Dfyydt(Nn-1);
%Dfyydt(1)=Dfyydt(2);


pyy(1:Nn)=0;
%pyy1(1:Nn)=0;
%pyy2(1:Nn)=0;
%pxx(1:Nn)=0;
for nn = wnStart+1:Nn-1
    
    pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
%    pyy1(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy1(nn));
%    pyy2(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy2(nn));
 %   pxx(nn) = (normpdf(Z(nn),0, 1))/abs(Dfxx(nn));
end


toc

ItoHermiteMeanVar=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
ItoHermiteMeanAsset=Mean
disp('true Mean only applicable to standard SV mean reverting type models otherwise disregard');
%TrueMean=theta+(v0-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average
%yy0


theta1=1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%rng(29079137, 'twister')
rng(15898837, 'twister')
paths=100000;
V(1:paths)=v00;  %Original process monte carlo.
X=0.0;
X(1:paths)=x0;
alpha1=0;
alpha2=1;
a=mu1X;
b=mu2X;
rho=0;
sigma1=sigmaX;
gammaV=.5;
Random1(1:paths)=0;
Random2(1:paths)=0;

for ttM=1:TtM
Random1=randn(size(Random1));
Random2=randn(size(Random2));



X(1:paths)=X(1:paths)+ ...
    (a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dtM + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) * sqrt(dtM) + ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*sqrt(dtM) + ...
    (a*alpha1* X(1:paths).^(alpha1-1)+b*alpha2* X(1:paths).^(alpha2-1)).* ...
    (((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dtM^2/2)+ ...
    (sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) *(1-1/sqrt(3)).*dtM^1.5+ ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*(1-1/sqrt(3)).*dtM^1.5))+ ...
    .5*(a*alpha1*(alpha1-1)* X(1:paths).^(alpha1-2)+b*alpha2*(alpha2-1).* X(1:paths).^(alpha2-2)).* ...
    ( sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX)) *dtM^2/2 + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
    ((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random1(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*(Random1(1:paths).^2-1) * dtM/2 + ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths)*dtM/2)+ ...
    .5*sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
    (sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random1(1:paths).*1/sqrt(3).*dtM^1.5)+ ...
    sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
    ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random1(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
    sigma0*V(1:paths).^gamma.*Random1(1:paths).*Random2(1:paths)*dtM/2)+ ...
    .5*sqrt(1-rho^2)* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
    (sigma0^2*V(1:paths).^(2*gamma).*Random1(1:paths)*1/sqrt(3)*dtM^1.5)+ ...
    sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
    rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random1(1:paths)*1/sqrt(3)*dtM^1.5+ ...
    rho* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
    ((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random2(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
    sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths) * dtM/2 + ...
    rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*(Random2(1:paths).^2-1)*dtM/2)+ ...
    .5*rho* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
    (sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random2(1:paths).*1/sqrt(3).*dtM^1.5)+ ...
    rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
    ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths) * 1/sqrt(3).* dtM^1.5 + ...
    sigma0*V(1:paths).^gamma.*(Random2(1:paths).^2-1)*dtM/2)+ ...
    .5*rho* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
    sigma0^2.*V(1:paths).^(2*gamma).*Random2(1:paths) * 1/sqrt(3).* dtM^1.5+ ...
    rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
    rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random2(1:paths)*1/sqrt(3)*dtM^1.5;
    
VBefore=V;

    V(1:paths)=V(1:paths)+ ...
        (mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dtM + ...
        sigma0*V(1:paths).^gamma .*Random2(1:paths)*sqrt(dtM) + ...
        (mu1.*beta1*V(1:paths).^(beta1-1) + mu2.*beta2.*V(1:paths).^(beta2-1)).* ...
        ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dtM^2/2 + ...
        sigma0*V(1:paths).^gamma .*Random2(1:paths)*(1-1/sqrt(3))*dtM^1.5) + ...
        .5*(mu1.*beta1.*(beta1-1).*V(1:paths).^(beta1-2) + mu2.*beta2.*(beta2-1).*V(1:paths).^(beta2-2)).* ...    
        sigma0^2.*V(1:paths).^(2*gamma).*dtM^2/2 + ...
        sigma0*gamma*V(1:paths).^(gamma-1) .* ...
        ((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths).*1/sqrt(3)*dtM^1.5 + ...
        sigma0.*V(1:paths).^gamma .*(Random2(1:paths).^2-1)*dtM/2) + ...
        .5*sigma0*gamma*(gamma-1).*V(1:paths).^(gamma-2) .* ...
        sigma0^2.*V(1:paths).^(2*gamma) .*Random2(1:paths).*1/sqrt(3)*dtM^1.5;

%Vm=.5*VBefore+.5*V;
Vm=V;
    
%    X(1:paths)=X(1:paths)+ ...
%    sigma1* Vm(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) * sqrt(dtM) + ...
%    sigma1* Vm(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
%    (sigma1* Vm(1:paths).^gammaV.* X(1:paths).^gammaX .*(Random1(1:paths).^2-1) * dtM/2)+ ...
%    .5*sigma1* Vm(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
%    (sigma1^2* Vm(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random1(1:paths).*1/sqrt(3).*dtM^1.5);%+ ...

    
    
    
end

%SVolMeanAnalytic=thetaV+(V0-thetaV)*exp(-kappaV*dt*Tt)
SVolMeanAnalytic=theta+(v00-theta)*exp(-kappa*dt*Tt)
SVolMeanMC=sum(V(1:paths))/paths
AssetMeanAnalytic=x0
AssetMeanMC=sum(X(1:paths))/paths

% theta
% v00
% kappa
% dt
% Tt
% T


%MeanX

disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(v00-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average

cq

MaxCutOff=30;
NoOfBins=.25*300*ceil(10*gammaX)*T;%round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too 
[XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti_NEW(X,paths,NoOfBins,MaxCutOff );
plot(xx(wnStart+1:Mm-1),pxx(wnStart+1:Mm-1),'r',IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'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 = %.2f,thetaX=%.2f,kappaX=%.2f,gammaX=%.2f,sigmaX=%.2f,v0 =%.2f,kappa=%.2f,theta=%.2f,gamma=%.2f,sigma0=%.2f,T=%.2f,dt=%.3f,M=%.4f',x0,thetaX,kappaX,gammaX,sigmaX,v00,kappa,theta,gamma,sigma0,T,dt,ItoHermiteMeanAsset));
%,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.');


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

NoOfBins=round(2*500*gamma^2*4*sigma0/sqrt(SVolMeanMC)/(1+kappa))/10;
[VDensity,IndexOutV,IndexMaxV] = MakeDensityFromSimulation_Infiniti_NEW(V,paths,NoOfBins,MaxCutOff );

%[YdtDensity,IndexOutYdt,IndexMaxYdt] = MakeDensityFromSimulation_Infiniti_NEW(YYdt,paths,NoOfBins,MaxCutOff );
plot(yy(wnStart+1:Nn-1),pyy(wnStart+1:Nn-1),'r',IndexOutV(1:IndexMaxV),VDensity(1:IndexMaxV),'g');
title(sprintf('v0 =%.2f,kappa=%.2f,theta=%.2f,gamma=%.2f,sigma0=%.2f,T=%.2f,dt=%.3f,M=%.4f',v00,kappa,theta,gamma,sigma0,T,dt,ItoHermiteMeanVar));



end

.
.
.

function [cH0,cH,cVdt] = CalculateVPathStepIntegral04(cprev,c,V0,kappa,theta,gamma,dt,Ts,cH0,cH,SeriesOrder)


    gammaV2New=1;
    [vw0,dvw0] = CalculatevgammaV2AndDerivatives(c(1),gamma,gammaV2New,SeriesOrder);
    %[V2gamma0,V2gamma1] = CalculateDriftbCoeffs08A(vw0,dvw0,cmid,SeriesOrder);
    [cva] = CalculateDriftbCoeffs08AB(vw0,dvw0,c,SeriesOrder);

    cv0=cva(1);
    cv(1:SeriesOrder-1)=cva(2:SeriesOrder);
    
    gammaV2New=1;
    [vw0,dvw0] = CalculatevgammaV2AndDerivatives(cprev(1),gamma,gammaV2New,SeriesOrder);
    %[V2gamma0,V2gamma1] = CalculateDriftbCoeffs08A(vw0,dvw0,cmid,SeriesOrder);
    [cvaprev] = CalculateDriftbCoeffs08AB(vw0,dvw0,cprev,SeriesOrder);
    
    cv0prev=cvaprev(1);
    cvprev(1:SeriesOrder-1)=cvaprev(2:SeriesOrder);
    
    
    
    SeriesOrder=SeriesOrder-1;
    [fH0,fH] = ConvertZCoeffsToHCoeffs(cv0,cv,SeriesOrder);   
    
    
    [fH0prev,fHprev] = ConvertZCoeffsToHCoeffs(cv0prev,cvprev,SeriesOrder);    
 
    %dH0=(V0-theta)*exp(-kappa*Ts*dt)-(V0-theta)*exp(-kappa*(Ts-1)*dt)
    dH0=theta*dt+(V0-theta)/kappa*(exp(-kappa*(Ts-1)*dt)-(exp(-kappa*Ts*dt)));

    dH(1:SeriesOrder)=sign(exp(0*kappa*dt).*fH(1:SeriesOrder)-exp(-1*kappa*dt).*(fHprev(1:SeriesOrder))).*sqrt(abs(sign(fH(1:SeriesOrder)).*(exp(0*kappa*dt).*fH(1:SeriesOrder)).^2-sign(exp(-kappa*dt).*(fHprev(1:SeriesOrder))).*(exp(-2*1*kappa*dt).*(fHprev(1:SeriesOrder)).^2)));    
 
    cH0(Ts)=dH0;
    cH(Ts,1:SeriesOrder)=exp(kappa*dt*Ts).*dH(1:SeriesOrder);
 
 
     Coeff1(1:Ts)=0;
    Coeff2(1:Ts)=0;
    

 for tt=1:Ts
%     Coeff02=Coeff02+(theta+(V0-theta).*exp(-1*kappa*dt*(tt)))*dt;
     for ss=tt:Ts
         
         Coeff2(tt)=Coeff2(tt)+(exp(-1*kappa*dt*(ss))).*dt;
         %Coeff02=Coeff02+theta*(1-exp(-1*kappa*dt*(ss)));

     end
 end
 
 Coeff2(:)=Coeff2(:).^2;
 
  for tt=1:Ts-1
 %     Coeff01=Coeff01+(theta+(V0-theta).*exp(-1*kappa*dt*(tt)))*dt;
     for ss=tt:Ts-1
         
         Coeff1(tt)=Coeff1(tt)+(exp(-1*kappa*dt*(ss))).*dt;
         %Coeff01=Coeff01+theta*(1-exp(-1*kappa*dt*(ss)));
         

     end
 end
 
 Coeff1(:)=Coeff1(:).^2;
    
    eH1(1:SeriesOrder)=0;
    eH2(1:SeriesOrder)=0;
    for tt=1:Ts-1
  %  ts=Ts-1-tt+1;
    eH1(1:SeriesOrder)=eH1(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*Coeff1(tt).*cH(tt,1:SeriesOrder).^2;   
    end 
eH0=0;
eH2(1:SeriesOrder)=0;
for tt=1:Ts
    %*dt;
    
    eH2(1:SeriesOrder)=eH2(1:SeriesOrder)+sign(cH(tt,1:SeriesOrder)).*Coeff2(tt).*cH(tt,1:SeriesOrder).^2;   
end

eH0=cH0(Ts);


eH(1:SeriesOrder)=0.0;

eH(1:SeriesOrder)=sign(eH2(1:SeriesOrder)).*sqrt(abs(sign(eH2(1:SeriesOrder)).*eH2(1:SeriesOrder)))-sign(eH1(1:SeriesOrder)).*sqrt(abs(sign(eH1(1:SeriesOrder)).*eH1(1:SeriesOrder)));
[c0vdt,cvdt] = ConvertHCoeffsToZCoeffs(eH0,eH,SeriesOrder);% 




 cVdt(1)=c0vdt;
 cVdt(2:6)=cvdt(1:5);

% str=input('Look at numbers');

end

.
.
.

function [b] = CalculateDriftbCoeffs08AB(wMu0dt,dwMu0dtdw,a,SeriesOrder)



b(1)=wMu0dt;

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

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

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

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

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


if(SeriesOrder>=8)
  b(8)=1/(720*7)*(5040* a(8)*dwMu0dtdw(1)+5040* a(4)* a(5)*dwMu0dtdw(2)+5040* a(3)* a(6) *dwMu0dtdw(2)+5040 *a(2)* a(7) *dwMu0dtdw(2)+ ...
      2520* a(3)^2* a(4) *dwMu0dtdw(3)+2520* a(2)* a(4)^2 *dwMu0dtdw(3)+5040 *a(2)* a(3)* a(5) *dwMu0dtdw(3)+2520* a(2)^2 *a(6) *dwMu0dtdw(3)+ ...
      840* a(2)* a(3)^3 *dwMu0dtdw(4)+2520* a(2)^2* a(3)* a(4) *dwMu0dtdw(4)+840* a(2)^3* a(5) *dwMu0dtdw(4)+ ...
      420* a(2)^3* a(3)^2 *dwMu0dtdw(5)+210* a(2)^4* a(4) *dwMu0dtdw(5)+ ...
      42* a(2)^5* a(3) *dwMu0dtdw(6)+a(2)^7 *dwMu0dtdw(7));
  
end
  
if(SeriesOrder>=9)
  b(9)=1/(720*7*8)*(40320* a(9) *dwMu0dtdw(1)+20160* a(5)^2 *dwMu0dtdw(2)+40320* a(4)* a(6) *dwMu0dtdw(2)+40320* a(3) *a(7) *dwMu0dtdw(2)+ ...
      40320* a(2) *a(8) *dwMu0dtdw(2)+ ...
      20160* a(3)* a(4)^2 *dwMu0dtdw(3)+20160* a(3)^2* a(5) *dwMu0dtdw(3)+40320* a(2)* a(4)* a(5) *dwMu0dtdw(3)+40320* a(2)* a(3)* a(6) *dwMu0dtdw(3)+ ...
      20160* a(2)^2* a(7) *dwMu0dtdw(3)+ ...
      1680* a(3)^4 *dwMu0dtdw(4)+20160* a(2)* a(3)^2* a(4) *dwMu0dtdw(4)+10080* a(2)^2* a(4)^2 *dwMu0dtdw(4)+20160* a(2)^2* a(3) *a(5) *dwMu0dtdw(4)+ ...
      6720* a(2)^3* a(6) *dwMu0dtdw(4)+ ...
      3360* a(2)^2* a(3)^3 *dwMu0dtdw(5)+6720* a(2)^3 *a(3)* a(4) *dwMu0dtdw(5)+1680* a(2)^4* a(5) *dwMu0dtdw(5)+ ...
      840* a(2)^4* a(3)^2 *dwMu0dtdw(6)+336* a(2)^5* a(4) *dwMu0dtdw(6)+ ...
      56* a(2)^6* a(3) *dwMu0dtdw(7)+a(2)^8 *dwMu0dtdw(8));
end
% 
% b(10)=1/(720*7*8*9)*(362880* a(10) *dwMu0dtdw(1)+ ...
%     362880* a(5)* a(6) *dwMu0dtdw(2)+ 362880* a(4)* a(7) *dwMu0dtdw(2) +362880* a(3)* a(8) *dwMu0dtdw(2) +362880* a(2) *a(9) *dwMu0dtdw(2) + ...
%     60480* a(4)^3 *dwMu0dtdw(3)+362880* a(3)* a(4)* a(5) *dwMu0dtdw(3)+181440* a(2)* a(5)^2 *dwMu0dtdw(3)+181440* a(3)^2* a(6) *dwMu0dtdw(3) + ...
%     362880* a(2)* a(4)* a(6) *dwMu0dtdw(3)+362880 *a(2)* a(3)* a(7) *dwMu0dtdw(3) +181440 *a(2)^2* a(8) *dwMu0dtdw(3)+ ...
%     60480* a(3)^3* a(4) *dwMu0dtdw(4)+181440* a(2)* a(3)* a(4)^2 *dwMu0dtdw(4)+181440* a(2)* a(3)^2 *a(5) *dwMu0dtdw(4)+ ...
%     181440* a(2)^2* a(4)* a(5) *dwMu0dtdw(4)+181440* a(2)^2* a(3)* a(6) *dwMu0dtdw(4)+60480* a(2)^3* a(7) *dwMu0dtdw(4)+ ...
%     15120* a(2)* a(3)^4 *dwMu0dtdw(5)+90720* a(2)^2* a(3)^2* a(4) *dwMu0dtdw(5)+30240* a(2)^3* a(4)^2 *dwMu0dtdw(5)+ ...
%     60480 *a(2)^3* a(3)* a(5) *dwMu0dtdw(5)+15120* a(2)^4* a(6) *dwMu0dtdw(5)+ ...
%     10080* a(2)^3* a(3)^3 *dwMu0dtdw(6)+15120* a(2)^4 *a(3)* a(4) *dwMu0dtdw(6)+3024* a(2)^5* a(5) *dwMu0dtdw(6)+ ...
%     1512* a(2)^5* a(3)^2 *dwMu0dtdw(7)+504* a(2)^6* a(4) *dwMu0dtdw(7)+ ...
%     72* a(2)^7 *a(3) *dwMu0dtdw(8)+a(2)^9 *dwMu0dtdw(9));
% 
% 
% 
% 
% b(11)=1/(720*7*8*9*10)*(3628800* a(11)* dwMu0dtdw(1)+ ...
%     1814400 * a(6)^2 *dwMu0dtdw(2)+3628800 *a(5)* a(7)* dwMu0dtdw(2)+3628800* a(4)* a(8)* dwMu0dtdw(2)+3628800 *a(3)* a(9)* dwMu0dtdw(2)+ ...
%     3628800* a(2)* a(10)* dwMu0dtdw(2)+ ...
%     1814400* a(4)^2* a(5)* dwMu0dtdw(3)+1814400* a(3)* a(5)^2* dwMu0dtdw(3)+3628800* a(3)* a(4)* a(6)* dwMu0dtdw(3)+3628800* a(2)* a(5)* a(6) *dwMu0dtdw(3)+ ...
%     1814400* a(3)^2* a(7)* dwMu0dtdw(3)+3628800* a(2)* a(4)* a(7)* dwMu0dtdw(3)+3628800 *a(2)* a(3)* a(8)* dwMu0dtdw(3)+1814400* a(2)^2* a(9)* dwMu0dtdw(3)+ ...
%     907200* a(3)^2 *a(4)^2* dwMu0dtdw(4)+604800* a(2)* a(4)^3* dwMu0dtdw(4)+604800* a(3)^3* a(5)* dwMu0dtdw(4)+3628800* a(2)* a(3)* a(4)* a(5)* dwMu0dtdw(4)+ ...
%     907200* a(2)^2* a(5)^2* dwMu0dtdw(4)+1814400* a(2)* a(3)^2* a(6)* dwMu0dtdw(4)+1814400* a(2)^2* a(4)* a(6)* dwMu0dtdw(4)+ ...
%     1814400* a(2)^2* a(3)* a(7)* dwMu0dtdw(4)+ 604800* a(2)^3* a(8) *dwMu0dtdw(4)+ ...
%     30240 *a(3)^5* dwMu0dtdw(5)+604800* a(2)* a(3)^3* a(4)* dwMu0dtdw(5)+907200 *a(2)^2* a(3)* a(4)^2* dwMu0dtdw(5)+907200* a(2)^2* a(3)^2* a(5) *dwMu0dtdw(5)+ ...
%     604800* a(2)^3* a(4)* a(5)* dwMu0dtdw(5)+604800* a(2)^3* a(3)* a(6)* dwMu0dtdw(5)+151200* a(2)^4* a(7)* dwMu0dtdw(5)+ ...
%     75600* a(2)^2* a(3)^4* dwMu0dtdw(6)+302400 *a(2)^3* a(3)^2* a(4)* dwMu0dtdw(6)+75600* a(2)^4 *a(4)^2* dwMu0dtdw(6)+151200* a(2)^4* a(3)* a(5) *dwMu0dtdw(6)+ ...
%     30240* a(2)^5* a(6)* dwMu0dtdw(6)+ ...
%     25200* a(2)^4* a(3)^3* dwMu0dtdw(7)+30240* a(2)^5* a(3)* a(4)* dwMu0dtdw(7)+5040* a(2)^6* a(5)* dwMu0dtdw(7)+ ...
%     2520* a(2)^6* a(3)^2* dwMu0dtdw(8)+720* a(2)^7* a(4)* dwMu0dtdw(8)+ ...
%     90 *a(2)^8* a(3)* dwMu0dtdw(9)+ ...
%     a(2)^10* dwMu0dtdw(10));
% 


end

.
.
.

function [vw0,dvw0] = CalculatevgammaV2AndDerivatives(c0,gamma,gammaV2,SeriesOrder)


    
    vw0=((1-gamma).*c0).^(gammaV2/(1-gamma));
    dvw0(1)=gammaV2*((1-gamma).*c0).^(gammaV2/(1-gamma)-1);
    dvw0(2)=gammaV2*(gammaV2/(1-gamma)-1).*((1-gamma).*c0).^(gammaV2/(1-gamma)-2)*(1-gamma);
    dvw0(3)=gammaV2*(gammaV2/(1-gamma)-1).*(gammaV2/(1-gamma)-2).*((1-gamma).*c0).^(gammaV2/(1-gamma)-3)*(1-gamma)^2;
    dvw0(4)=gammaV2*(gammaV2/(1-gamma)-1).*(gammaV2/(1-gamma)-2)*(gammaV2/(1-gamma)-3)*((1-gamma).*c0).^(gammaV2/(1-gamma)-4)*(1-gamma)^3;
    dvw0(5)=gammaV2*(gammaV2/(1-gamma)-1).*(gammaV2/(1-gamma)-2)*(gammaV2/(1-gamma)-3)*(gammaV2/(1-gamma)-4)*((1-gamma).*c0).^(gammaV2/(1-gamma)-5)*(1-gamma)^4;
    dvw0(6)=gammaV2*(gammaV2/(1-gamma)-1).*(gammaV2/(1-gamma)-2)*(gammaV2/(1-gamma)-3)*(gammaV2/(1-gamma)-4)*(gammaV2/(1-gamma)-5)*((1-gamma).*c0).^(gammaV2/(1-gamma)-6)*(1-gamma)^5;
    

end

.
.
.

function [c] =SeriesProductAB(a,b,SeriesOrder)


c(1)=a(1)*b(1);
c(1:SeriesOrder)=0.0;
for nn=1:SeriesOrder
    c(nn)=c(nn)+a(1)*b(nn);
    c(nn)=c(nn)+b(1)*a(nn);
    for kk=2:nn-1
        c(nn)=c(nn)+a(kk)*b(nn-kk+1);
    end
end


% c0=a0*b0;
% c(1:SeriesOrder)=0.0;
% for nn=1:SeriesOrder
%     c(nn)=c(nn)+a0*b(nn);
%     c(nn)=c(nn)+b0*a(nn);
%     for kk=1:nn-1
%         c(nn)=c(nn)+a(kk)*b(nn-kk);
%     end
% end


end

.
.
.

function [c] = SeriesReciprocalAB(a,SeriesOrder)


Mul=1/(a(1));
a(1:SeriesOrder)=a(1:SeriesOrder)/a(1);
a(1)=1;



c(1)=1;
c(2:SeriesOrder)=0;
for nn=2:SeriesOrder
    c(nn)=c(nn)-a(nn);
    for kk=2:nn-1
        c(nn)=c(nn)-c(kk)*a(nn-kk+1);
    end
end
    

c(1:SeriesOrder)=c(1:SeriesOrder)*Mul;


end

.
.
.

function [aH] = ConvertZCoeffsToHCoeffsAB(a,SeriesOrder)

if(SeriesOrder==4)
aH(1)=a(1)+a(3);
aH(2)=a(2)+3*a(4);
aH(3)=a(3);
aH(4)=a(4);
end


if(SeriesOrder==5)
aH(1)=a(1)+a(3)+3*a(5);
aH(2)=a(2)+3*a(4);
aH(3)=a(3)+6*a(5);
aH(4)=a(4);
aH(5)=a(5);
end


if(SeriesOrder==6)
aH(1)=a(1)+a(3)+3*a(5);
aH(2)=a(2)+3*a(4)+15*a(6);
aH(3)=a(3)+6*a(5);
aH(4)=a(4)+10*a(6);
aH(5)=a(5);
aH(6)=a(6);
end



if(SeriesOrder==7)
aH(1)=a(1)+a(3)+3*a(5)+15*a(7);
aH(2)=a(2)+3*a(4)+15*a(6);
aH(3)=a(3)+6*a(5)+45*a(7);
aH(4)=a(4)+10*a(6);
aH(5)=a(5)+15*a(7);
aH(6)=a(6);
aH(7)=a(7);
end


if(SeriesOrder==8)
aH(1)=a(1)+a(3)+3*a(5)+15*a(7);
aH(2)=a(2)+3*a(4)+15*a(6)+105*a(8);
aH(3)=a(3)+6*a(5)+45*a(7);
aH(4)=a(4)+10*a(6)+105*a(8);
aH(5)=a(5)+15*a(7);
aH(6)=a(6)+21*a(8);
aH(7)=a(7);
aH(8)=a(8);
end



end

.
.
.

function [wMu0dt,dwMu0dtdw,wMu1dt,dwMu1dtdw] = BesselDriftAndDerivatives08(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,ExpnOrder)



NoOfTerms=3;

YqCoeff(1:NoOfTerms)=0.0;
Fp(1:NoOfTerms)=0.0;


YqCoeff(1)=mu1;
YqCoeff(2)=mu2;
YqCoeff(3)=-.5*sigma0^2*gamma;
Fp(1)=(beta1-gamma)/(1-gamma);
Fp(2)=(beta2-gamma)/(1-gamma);
Fp(3)=-1;



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

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


for mm=1:NoOfTerms

    wMu0dt0=YqCoeff(mm).*((1-gamma)*w0).^Fp(mm);
    for nn=1:ExpnOrder+1
        if(nn==1)
            %dwMu0dt(nn)=wMu0dt0*Fp(mm)*(1-gamma)/((1-gamma)*w0);
            dwMu0dt(nn)=wMu0dt0*Fp(mm)/(w0);
        else
        dwMu0dt(nn)=dwMu0dt(nn-1)*(Fp(mm)-(nn-1))/(w0);
        end
    end
    wMu0dt=wMu0dt+wMu0dt0;
    for nn=1:ExpnOrder+1
        dwMu0dtdw(nn)=dwMu0dtdw(nn)+dwMu0dt(nn);
    end
wMu1dt=dwMu0dtdw(1);
dwMu1dtdw(1:ExpnOrder)=dwMu0dtdw(2:ExpnOrder+1);
    
end

%wMu0dt
%dwMu0dtdw

%str=input('Look at drift derivatives');

end

.
.
.

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


NoOfTerms=19;
%NoOfTerms=9;

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

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

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

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

for mm=1:NoOfTerms

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



end

.
.
.

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

.
.
.

Amin · June 21st, 2022, 8:26 am

Before posting the graphs from my work, I would like to take this opportunity to thank those friends who with their help made my research possible. Several years ago, I had to spend a huge amount of time everyday by going into remote parts of Lahore city and try to find good food and water. Mind control agencies would drug food, water and beverages in the entire markets and it was extremely difficult to get good food and water. But as I tried to do good research, several American friends tried to help me by asking mind control agencies to become better and stop employing such tactics. Though food and water is drugged at a lot of places even today, it is easier these days than ever before to go out and be careful and easily be able to get good water and food (though I would still get hit occasionally sometimes). All other forms of torture also decreased with time especially when I would mention the torture here and friends would try to help me by asking secret services to become more civilized. I can never imagine that I could do any of the research I posted on this forum if circumstances were as bad as they were several years ago. I simply cannot have enough gratitude for American friends and others who have indirectly helped me in the past and thus made my research possible.

.
Friends, when I wrote the above lines, I was very careful to try to make it a positive moment that everybody would like and not stress about negative issues and rather downplay them. I was thinking that mind control agents might also become better when people like my work and decrease mind control. On the contrary mind control agencies sharply increased mind control. In past few weeks, if I like some drink or food somewhere, it is made sure to keep that drink or food item permanently drugged for several weeks so if if I try it again, I would be hit with mind control drugs. All the bottled water within several kilometer radius of where I live (and all of the nearby markets) is actively kept drugged with mind control drugs. Petrol is very expensive and it is very difficult to go to remote parts of the city to get food and water once or twice everyday. Once I am hit with drugged food, it is difficult to do anything positive or constructive for next few hours. After having drugged food, I would become very lethargic or have bouts of anxiety and would mostly like to sleep after that to recover and be able to start working again. Meanwhile mind control agents continue to take all sort of liberties with my body (other than trying to control the brain) like targeting my body with lasers, causing very severe itch in my back and sometimes even forcing urine or other such humiliating actions. I am quite used to all of these antics and usually never make a big deal out of them and continue my work like business as usual. But sometimes I may become upset and it was one of those moments when I wrote the following post. I was working on my research and taking notes when mind control agent gave me a very strong itch in my back out of nowhere. I requested him to not do that but he continued to increase the intensity. Unusually I was very upset and wrote the following post right there.

Friends, I have been begging the scum in American army to end their racist practices but it seems that racist scumbags in American army consider it my weakness and want to somehow continue my persecution forcefully through the rest of my life. Though many pseudo-patriotic(Patriotism is true and intelligent love for your country and people of your country and please do not equate it for love for scum armed forces and their weaponry). Americans may not like my speaking the truth about scum American army, matter of the fact is that American army is a large group of racist scumbags. Basically most large armies if they are not kept properly in check become scum. Is it not interesting that many Americans love their army but think of Pakistani army as full of scum bags. Pakistanis love their army and call American army full of scumbags. No, all large armies, yours and mine are equally full of scumbags. Basically scumbags in American army suck the dick of rich and racist jews so as to gain jobs paying millions of dollars in their companies after their retirement. Friends believe it or not it is really all about Benjamins. When tens of millions of Benjamins are given by corrupt racist crooks in American army to scumbag generals in Pakistani army, they start retarding hundreds of talented people in my country and they love to do it because they can never make such a large number of Benjamins any other way. So this chain of scumbags starting from racist scumbags in American army to scumbags in Pakistani army(other countries army) all wanting big millions of Benjamins. destroys talent all over the world. Please thwart the scum in American army and ask them to respect human beings.

.
I was also very heartbroken and all my enthusiasm was gone and I did only some lukewarm research on SV SDEs (even though I had several other very good ideas) after that and started to work on my market trading research. Also because I was hit with drugged food several times. For example, Food I had this morning was also drugged (they had drugged food in the bakery before I reached there this morning) and I slept for a few hours after that despite that I had a more than usual sleep last night.
Though I have started to like research on market trading, research on SDEs and stochastic processes is really close to my heart and I love working on it. In contrast to market trading research which can possibly have huge utility in the future, I do not get paid for my research on SDEs and I do not have any clients and it is mostly a thankless business other than that I love working on research with SDEs and I really want to continue it . I want to request friends to please force American mind control agency to please treat me with my human dignity intact, stop taking liberties with my body and stop drugging food and particularly water in all of the surrounding areas of the city. And end my mind control.

Amin · June 21st, 2022, 10:23 am

As addition to previous post, they always have mind-control gas in the air-conditioners in the rooms I use whenever I turn them on.

Amin · June 22nd, 2022, 5:09 am

Friends, I will be posting a new version of analytic solution to pair of SV SDEs. It is simpler since it depends only on addition of squared variances across hermite polynomials and yields somewhat better results than solution of Fokker-planck equation I posted earlier. I hope it would be interesting for friends to play around and explore it. It is late now but I would be posting it tomorrow.

Friends, I had another version of 2D SV SDEs solution based on addition of squared coefficients of hermite polynomials in the form of a 2D series. When I worked on that solution, it was not working properly and was giving similar or rather slightly better results than our general 2D SV solution. I believe the reason again was that dz-integral had to be derived from corresponding proper stochastic dt-integral. In next few days, I will be working with this solution. I hope that it would work perfectly for CEV and lognormal type asset SDEs but for more general SDEs, we could use this method of addition of squared coefficients at the start(for a few initial steps) without the need for a miniscule step for the asset SDE at the start. This is very interesting and I hope most friends would really like it. Since I have still not completed it perfectly there is a slight chance that something may go wrong but I am very sure that it will work perfectly. If it does work, I think it would be faster and more accurate than a full-fledged 2D SV SDEs solution implementation that we already have at least for first year or first few years. I will resume work on it today and will come back with results and the new matlab program in a few days.

Amin · June 24th, 2022, 4:24 pm

Friends, I am taking a break from my research on SDEs for one or two weeks. I continued to work on SDEs till yesterday but I was not getting any real work done and I was feeling burnt out. In past whenever I switched from my research on market trading to work on SDEs and vice versa, my productivity sharply increased. Especially for the program for SDEs that I mentioned in my previous post, I have to first work out the stochastic dt-integral of functions of the mean-reverting SDE variable. We have previously learnt how to do it for SDE variable but I have to extend it to functions of SDE variable before I I could complete the work in previous post.
I have several other ideas for research in stochastics that I really want to try to find out if some success could be found. For example, I have some ideas about finding the Z-series representation of a stochastic variable once we are given its first few moments. I have some vague ideas that such a solution could possibly be found out of the zeros of an appropriate polynomial(zeros would be solved for numerically) but of course, I could not be sure if the idea actually works out. If it does work out, it would be very interesting since we could provide an estimate of moments or alternatively cumulants and our function would simply find the Z-series representation of the stochastic variable which is very handy for a lot of analytics.
I will be spending one or two weeks on market trading research and then come back to research on SDEs and related stochastic processes.
Another thing I want to tell friends that I lost my phone several days ago. I rarely use my phone and do most of the work on my laptop. I do not have many friends who would call me and mostly 99% of the time, I receive calls from members of my family only. I was robbed of my earlier phone when I was walking in the morning and three guys on a motorcycle came close, showed me a pistol and took my phone with them. I bought a new phone but would leave it in table drawers in my room while going for morning walk. Since I do not use the phone much, many times, it would remain in the drawers for the whole day after I come back. Today when I realized that I do not have used my phone for past few days, I looked for it in the drawers and everywhere else but could not find it. I have no idea where it is gone. I am telling friends because there is a possibility that my phone might be with mind control crooks and as they are desperate to concoct some type of wrongdoing or evil intention from my side so that they could continue my mind control so there is a possibility they might have used my phone in some way that could help them find a valid reason for my persecution. I want to tell friends that I do not use my phone on internet and have not talked to anyone other than my family for several months. If my phone has been used in some of these ways, it has been plotted by mind control agents and I have not done it especially I have not used my phone at all in past three or four days. Since mind control agents are desperate to continue my persecution, they can go to any limit to concoct evidence to falsely establish that I am a bad person or I have bad intentions.

Amin · June 26th, 2022, 9:44 am

Friends, I have found a surprisingly simple formula for a probability distribution based on gaussian distribution and hermite polynomials where you can fit first few (six or eight or more) moments very easily using this distribution and it will be ensured that our probability distribution would have first few moments matching exactly the specified moments. I will write a detailed post tonight or early tomorrow.

Amin · June 27th, 2022, 2:16 pm

Some notations and preliminaries first.

Let us suppose we denote a standard gaussian density as [$]p(X) \, = \, \frac{1}{2 \pi} \exp(- \, \frac{ X^2 }{2})[$]

We also have hermite polynomials and first few hermite polynomials are denoted as
[$]H_0(X) \, = \, 1[$]
[$]H_1(X) \, = \, X[$]
[$]H_2(X) \, = \, X^2 \, - 1[$]
[$]H_3(X) \, = \, X^3 \, - 3 X[$]
[$]H_4(X) \, = \, X^4 \, - 6 X^2 \, +3[$]
and so on.
And I want to share an interesting property related to following integral
[$]\int_{-\infty}^{+\infty} X^m \, H_n(X) \, p(X) dx [$] Equation(A)
We know when m is odd and n is even or m is even and n is odd, the above integral always goes to zero.
But when m is odd and n is also odd, the above integral would always go to zero when n > m.
Similarly when m is even and n is also even, the above integral would again always go to zero, when n > m.
To explain it in words with simple examples, when fourth power of X is integrated in the above integral (eq A) with zeroth hermite, 2nd hermite or fourth hermite, there will be a finite value but when fourth power of X is integrated with sixth hermite, eighth hermite or tenth hermite, it will always go to zero. Similarly for odd powers, when fifth power of X is integrated there is a finite value when we integrate it with first hermite, third hermite or fifth hermite in the equation (A), but when we integrate fifth power of X with seventh, ninth or 11th hermite or higher odd hermites, the value would always go to zero.
The above property is interesting since in our forthcoming framework when we have fixed some moment in our equation, adding higher hermites(for higher moments) would not change the lower moments at all.

Now we want to fix first few moments to our proposed probability distribution so that all the moments are perfectly matched (a good property as I mentioned earlier would be that adding further higher moments would not change our lower moments that have already been matched).

The proposed probability distribution has a very simple form (upto sixth order) as

[$] q(X) \, = \, p(X) \, \big[ \, 1 \, + \, h_1 \, H_1(X)+ \, h_2/2 \, H_2(X)+ \, h_3/6 \, H_3(X)+ \, h_4/24 \, H_4(X)+ \, h_5/120 \, H_5(X)+ \, h_6/720 \, H_6(X) \, \big] \, [$]

here [$]h_1[$] , [$]h_2[$] etc are coefficients of hermites that will be adjusted to match the moments of above probability density form to our desired moments.
1st thing we know by inspection that above probability density equation integrates to one since all hermites (other than the first term) have an expected value of zero when integrated over normal density.
Now we write the first few raw moments of the above equation that can be verified by full-fledged integration(as I have done in my notes)

[$]E[X]\,= \int_{-\infty}^{+\infty} X \, q(X) dX \, = \, h_1[$]
[$]E[X^2]\,= \int_{-\infty}^{+\infty} X^2 \, q(X) dX \, = \,1 \, + \, h_2[$]
[$]E[X^3]\,= \int_{-\infty}^{+\infty} X^3 \, q(X) dX \, = \,3 \, h_1 \, + \, h_3[$]
[$]E[X^4]\,= \int_{-\infty}^{+\infty} X^4 \, q(X) dX \, = \, 3 \, +\,6 \, h_2 \, + \, h_4[$]
[$]E[X^5]\,= \int_{-\infty}^{+\infty} X^5 \, q(X) dX \, = \, 15 \, +\,10 \, h_3 \, + \, h_5[$]
[$]E[X^6]\,= \int_{-\infty}^{+\infty} X^6 \, q(X) dX \, = \, 15 \, +\,45 \, h_2 \, + \,15 \, h_4 \, + h_6[$]

So we see that we can keep on adding more hermite polynomials to fix higher moments and it will not alter our lower moments. However higher odd/even moments take a contribution from hermite polynomial terms associated with lower odd/even moments(that have already been fixed) up to the order of moment being calculated.
It seems a better idea to try matching moments on standard gaussian as opposed to non-standard gaussian as we do when we match moments with gram-charlier or edgeworth.
Please pardon any errors.

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