Serving the Quantitative Finance Community

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

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

August 30th, 2021, 4:03 pm

Friends, we are just getting familiar with the problem of SV option pricing with this new method and getting our orientation right with this new technique. The matlab program I posted is merely an initial framework on which we will build to solve the problem in a very good manner and it requires many many improvements in a lot of aspects. So we will continue to study and try to better solve the problem of SV option pricing for next few weeks and I will continue to post re-worked versions of the program after every few days. The present program has a lot of pitfalls but we will continue to improve it for at least several weeks until everything is properly finessed.
 
User avatar
Amin
Topic Author
Posts: 2853
Joined: July 14th, 2002, 3:00 am

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

August 30th, 2021, 9:37 pm

I have gained some more insight towards the solution to the problem. In previous programs, the shape of density is systematically different from monte carlo when viewed very closely. I think I know the reason why the density is slightly but systematically off from true density and am trying to find a robust solution. I tried a decent fix and made the graph below with changes to previous method. Now density fit is a lot better. I am attaching the graph of asset density. Please notice how density changes (asymmetrically) close to peak to stay close to the true density. 

I will be posting new matlab programs in another day or two.

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

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

August 31st, 2021, 6:33 pm

I am going to share the analytics for solution of 2D stochastic volatility equations. Before I give the solution, I want to mention that in typical SV system of SDEs, there are two SDEs out of which one (usually Stochastic volatility) is independent of the other SDE and its one dimensional PDE can be solved independently. If we had a system of two SDEs which are both dependent on each other, in that case the solution though still possible would be more difficult.

I have taken both asset and SV Sdes in bessel coordinates. SV SDE can be solved independently and its PDE is given as

[$]\frac{\partial P(w,t)}{\partial t} \, = \, -\frac{\partial [\mu_w \, P(w,t)]}{\partial w} \, +\, .5 \, \frac{\partial^2 [{\sigma_w}^2 \, P(w,t)]}{\partial w^2} \,[$]

When we solve this equation with conservation of mass as I showed in a previous post, we get


[$] -\mu_w \, \frac{\partial Z_1}{\partial w} - .5 {\sigma_w}^2 \big[ -Z_1  \, {(\frac{\partial Z_1}{\partial w})}^2- P(Z_1)\, {(\frac{\partial Z_1}{\partial w})}^3  \frac{\partial^2 w}{\partial {Z_1}^2} \, \big]  \,[$]
[$] +  \frac{\partial w_b}{\partial t} \frac{\partial Z_1}{\partial w} =0 [$]

In the above equation b in [$]w_b[$] denotes that it is the boundary point that is moving so as to keep the probability conserved up till that point. We will need the above equation because when we move to two dimensional stochastic volatility PDE, the whole of above equation would be integrated in another dimension but since it equates to zero, it will also remain zero even after an extra integration and we will eliminate it from 2D equation.

our two dimensional SV and asset partial differential equation in bessel coordinates that we want to solve is 
[$]\frac{\partial P(w,y,t)}{\partial t} \, = \, -\frac{\partial \big[\mu_y \, P(w,y,t) \big]}{\partial y} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(w,y,t) \big]}{\partial y^2}  \, -\frac{\partial \big[\mu_w \, P(w,y,t) \big]}{\partial w} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_w}^2 \, P(w,y,t) \big]}{\partial w^2} \,+ \,\rho \, \frac{\partial^2 \big[ \sigma_y \, \sigma_w  \, {V(w)}^{ \gamma_V} \, P(w,y,t) \big]}{\partial w \, \partial y} \, [$]

though correlation term is conceptually simple and when we do a two dimensional integration to preserve probability mass, this term becomes even more easy but correlation has some other issues when taking derivatives with respect of density of y with respect to Z_1 in diffusion terms. When there is no correlation our jacobian for two dimensional density when taken with respect to two independent standard nromals becomes

[$]\frac{\partial (Z_1,Z_2)}{\partial (w,y)}=\frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}[$]

since [$]Z_1[$] and [$]Z_2[$] are independent, we can write

[$]P(w,y)=P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} [$]

Now we set ourselves in the probability conservation set up, and we have to write a two dimensional integral from infinity(or rather zero ) to moving boundary (boundary might be a misnomer since this is not grid boundary but rather boundary of mass conserved till a certain grid point and this so called boundary is moving or alternatively the grid point is moving) points [$]w_b[$] and [$]y_b[$]. Please note that on a two dimensional grid there are a series of these points [$]w_b[$]'s and [$]y_b[$]'s. The grid in inependent SV dimension is given by index n and the grid on dependent asset dimension is given by (n,m) and respective moving boundaries in conservation set up are indicated as  [$]w_b(n)[$](that conserves mass till nth point on SV grid) and [$]y_b(n,m)[$] (that conserves mass till (n,m) points on the grid) . Please note that  [$]y_b(n,m)[$] is not a straight line and is curved in both indices.

We apply the conservation of probability mass till arbitrary nth and mth points on the grid as

[$] \frac{\partial \Big[\int_0^{w_b(n)} \, \int_0^{y_b(n,m)} \, P(w,y,t) \,\, dy \, dw  \Big]}{\partial t} \, =0 [$]
[$]=\int_0^{w_b(n)} \, \int_0^{y_b(n,m)} \, \frac{P(w,y,t)}{\partial t} \, \, dy \, dw   \, [$]
[$]+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(w_b(n),y,t) dy [$] 
[$]+ \int_0^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t} \, P(w,y_b,t) dw =0 [$]

Substituting the 2D uncorrelated PDE in the first RHS integral, we get the expression

[$]=\int_0^{w_b(n)} \, \int_0^{y_b(n,m)} \, \Big[\, -\frac{\partial \big[\mu_y \, P(w,y,t) \big]}{\partial y} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(w,y,t) \big]}{\partial y^2}  \, -\frac{\partial \big[\mu_w \, P(w,y,t) \big]}{\partial w} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_w}^2 \, P(w,y,t) \big]}{\partial w^2} \, \Big] \, \, dy \, dw   \, [$]
[$]+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(w_b(n),y,t) dy [$] 
[$]+ \int_0^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t} \, P(w,y_b,t) dw =0 [$]

We do the integrations in the first term to get
[$]=\int_0^{w_b(n)}  \, \Big[\, -\mu_y \, P(w,y,t)  \, +\, .5 \, \frac{\partial \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(w,y,t) \big]}{\partial y}  \, \Big] dw [$]
[$]\int_0^{y_b(n,m)} \Big[ -\mu_w \, P(w,y,t)  \, +\, .5 \, \frac{\partial \big[{\sigma_w}^2 \, P(w,y,t) \big]}{\partial w} \, \Big] \, \, dy   \, [$]
[$]+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(w_b(n),y,t) dy [$] 
[$]+ \int_0^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t} \, P(w,y_b,t) dw =0 [$]

Now we convert our coordinates to [$]Z_1[$] and [$]Z_2[$] instead of w and y as

[$]=\int_0^{w_b(n)}  \, \Big[\, -\mu_y \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}  \, +\, .5 \, \frac{\partial \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial y}  \, \Big] dw [$]
[$]\int_0^{y_b(n,m)} \Big[ -\mu_w \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}  \, +\, .5 \, \frac{\partial \big[{\sigma_w}^2 \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial w} \, \Big] \, \, dy   \, [$]
[$]+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(Z_1(n)) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} dy [$] 
[$]+ \int_0^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1(n)}{\partial w} \, \frac{\partial Z_2(m)}{\partial y} dw =0 [$]

which can be further simplified as

[$]=\int_0^{w_b(n)}  \, \Big[\, -\mu_y \,  \frac{\partial Z_2}{\partial y}  \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V}   \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big]  \, \Big] dZ_1 [$]
[$]\int_0^{y_b(n,m)} \Big[ -\mu_w  \, \frac{\partial Z_1}{\partial w}   \, +\, .5 \, {\sigma_w}^2   \, \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2 - {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \big]  \Big] \, \, dZ_2   \, [$]
[$]+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \,  \frac{\partial Z_1}{\partial w} \,  dZ_2 [$] 
[$]+ \int_0^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t}  \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0 [$]

But we notice that second and third term together are merely the solution of 1D independent PDE integrated in an extra dimension and since the original solution of both terms for independent 1D PDE was zero, the integrated term should also be zero. (Even if it were not zero, we could have found this solution from solving the 1D independent PDE), therefore
[$]\int_0^{y_b(n,m)} \Big[ -\mu_w  \, \frac{\partial Z_1}{\partial w}   \, +\, .5 \, {\sigma_w}^2   \, \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2 - {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \big]  \Big] \, \, dZ_2   \,[$]
[$]+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \,  \frac{\partial Z_1}{\partial w} \,  dZ_2=0 [$] 

Now we are left with the terms to solve as

[$]=\int_0^{w_b(n)}  \, \Big[\, -\mu_y \,  \frac{\partial Z_2}{\partial y}  \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V}   \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big]  \, \Big] dZ_1 [$]
[$]+ \int_0^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t}  \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0 [$]

For a given value of m, we can recursively solve the above PDE for each value of n. To make matters simple, I write the above equation as

[$]=\int_0^{w_b(n)}  \, \Big[\, -\mu_y \,  \frac{\partial Z_2}{\partial y}  \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V}   \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big]  \, \Big] dZ_1 [$]
[$]+ \int_0^{w_b(n-1)}  \frac{\partial y_b(n,m)}{\partial t}  \, \frac{\partial Z_2}{\partial y} dZ_1 [$]
[$]+ \int_{w_b(n-1)}^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t}  \, \frac{\partial Z_2}{\partial y} dZ_1 =0[$]

where in the above equation, I have expanded the last terms as sum of all previous terms and one current term to solve and now the equation can be solved as

[$] \frac{\partial y_b(n,m)}{\partial t} =\Bigg[ \int_0^{w_b(n)}  \, \Big[\, \mu_y \,  \frac{\partial Z_2}{\partial y}  \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V}   \, \big[ Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 + {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big]  \, \Big] dZ_1 [$]
[$]- \int_0^{w_b(n-1)}  \frac{\partial y_b(n,m)}{\partial t}  \, \frac{\partial Z_2}{\partial y} dZ_1 \Bigg] / \Big[ \int_{w_b(n-1)}^{w_b(n)}  \, \frac{\partial Z_2}{\partial y} dZ_1 \Big][$]

please note that [$]\frac{\partial y_b(n,m)}{\partial t} \Delta t[$] is the change in coordinates y(n,m) so as to conserve the mass in a small time step according to given PDE dynamics.
The last equation solves the 2D PDE in the asset dimension. Adding correlation term should not be difficult but then we have to make a few changes as correlation will appear in jacobian. Also derivatives of asset with respect to [$]Z_1[$] would have to be included something that was making my equation unstable when I tried it last time but I will try it again to see if I could make it work.
I will be posting the program for uncorrelated PDE based on these analytics tomorrow.
Friends, in the term
Now we convert our coordinates to [$]Z_1[$] and [$]Z_2[$] instead of w and y as

[$]=\int_0^{w_b(n)}  \, \Big[\, -\mu_y \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}  \, +\, .5 \, \frac{\partial \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial y}  \, \Big] dw [$]
[$]\int_0^{y_b(n,m)} \Big[ -\mu_w \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}  \, +\, .5 \, \frac{\partial \big[{\sigma_w}^2 \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial w} \, \Big] \, \, dy   \, [$]
[$]+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(Z_1(n)) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} dy [$] 
[$]+ \int_0^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1(n)}{\partial w} \, \frac{\partial Z_2(m)}{\partial y} dw =0 [$]

I have not properly handled the derivative  given by expression  [$]\frac{\partial \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial y} [$]
y is changing in both [$]Z_1[$] and [$]Z_2[$] directions. Also [$]{V(w)}^{2 \gamma_V}[$] also takes a derivative with respect to y through inverse derivative of y to w.  I wanted to fix these things tonight but I have been given mind control chemicals in food and I find it impossible to work tonight. I am in a very poor state mentally and totally unable to work. I also received an antipsychotic injection today in my arm. I will try fixing the errors in next 2-3 days and then post the program here. Please protest against my mind control. It is very painful to be unable to work when I need to work hard to solve the problem. 
 
User avatar
Amin
Topic Author
Posts: 2853
Joined: July 14th, 2002, 3:00 am

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

September 4th, 2021, 8:15 am

I am going to share the analytics for solution of 2D stochastic volatility equations. Before I give the solution, I want to mention that in typical SV system of SDEs, there are two SDEs out of which one (usually Stochastic volatility) is independent of the other SDE and its one dimensional PDE can be solved independently. If we had a system of two SDEs which are both dependent on each other, in that case the solution though still possible would be more difficult.

I have taken both asset and SV Sdes in bessel coordinates. SV SDE can be solved independently and its PDE is given as

[$]\frac{\partial P(w,t)}{\partial t} \, = \, -\frac{\partial [\mu_w \, P(w,t)]}{\partial w} \, +\, .5 \, \frac{\partial^2 [{\sigma_w}^2 \, P(w,t)]}{\partial w^2} \,[$]

When we solve this equation with conservation of mass as I showed in a previous post, we get


[$] -\mu_w \, \frac{\partial Z_1}{\partial w} - .5 {\sigma_w}^2 \big[ -Z_1  \, {(\frac{\partial Z_1}{\partial w})}^2- P(Z_1)\, {(\frac{\partial Z_1}{\partial w})}^3  \frac{\partial^2 w}{\partial {Z_1}^2} \, \big]  \,[$]
[$] +  \frac{\partial w_b}{\partial t} \frac{\partial Z_1}{\partial w} =0 [$]

In the above equation b in [$]w_b[$] denotes that it is the boundary point that is moving so as to keep the probability conserved up till that point. We will need the above equation because when we move to two dimensional stochastic volatility PDE, the whole of above equation would be integrated in another dimension but since it equates to zero, it will also remain zero even after an extra integration and we will eliminate it from 2D equation.

our two dimensional SV and asset partial differential equation in bessel coordinates that we want to solve is 
[$]\frac{\partial P(w,y,t)}{\partial t} \, = \, -\frac{\partial \big[\mu_y \, P(w,y,t) \big]}{\partial y} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(w,y,t) \big]}{\partial y^2}  \, -\frac{\partial \big[\mu_w \, P(w,y,t) \big]}{\partial w} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_w}^2 \, P(w,y,t) \big]}{\partial w^2} \,+ \,\rho \, \frac{\partial^2 \big[ \sigma_y \, \sigma_w  \, {V(w)}^{ \gamma_V} \, P(w,y,t) \big]}{\partial w \, \partial y} \, [$]

though correlation term is conceptually simple and when we do a two dimensional integration to preserve probability mass, this term becomes even more easy but correlation has some other issues when taking derivatives with respect of density of y with respect to Z_1 in diffusion terms. When there is no correlation our jacobian for two dimensional density when taken with respect to two independent standard nromals becomes

[$]\frac{\partial (Z_1,Z_2)}{\partial (w,y)}=\frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}[$]

since [$]Z_1[$] and [$]Z_2[$] are independent, we can write

[$]P(w,y)=P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} [$]

Now we set ourselves in the probability conservation set up, and we have to write a two dimensional integral from infinity(or rather zero ) to moving boundary (boundary might be a misnomer since this is not grid boundary but rather boundary of mass conserved till a certain grid point and this so called boundary is moving or alternatively the grid point is moving) points [$]w_b[$] and [$]y_b[$]. Please note that on a two dimensional grid there are a series of these points [$]w_b[$]'s and [$]y_b[$]'s. The grid in inependent SV dimension is given by index n and the grid on dependent asset dimension is given by (n,m) and respective moving boundaries in conservation set up are indicated as  [$]w_b(n)[$](that conserves mass till nth point on SV grid) and [$]y_b(n,m)[$] (that conserves mass till (n,m) points on the grid) . Please note that  [$]y_b(n,m)[$] is not a straight line and is curved in both indices.

We apply the conservation of probability mass till arbitrary nth and mth points on the grid as

[$] \frac{\partial \Big[\int_0^{w_b(n)} \, \int_0^{y_b(n,m)} \, P(w,y,t) \,\, dy \, dw  \Big]}{\partial t} \, =0 [$]
[$]=\int_0^{w_b(n)} \, \int_0^{y_b(n,m)} \, \frac{P(w,y,t)}{\partial t} \, \, dy \, dw   \, [$]
[$]+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(w_b(n),y,t) dy [$] 
[$]+ \int_0^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t} \, P(w,y_b,t) dw =0 [$]

Substituting the 2D uncorrelated PDE in the first RHS integral, we get the expression

[$]=\int_0^{w_b(n)} \, \int_0^{y_b(n,m)} \, \Big[\, -\frac{\partial \big[\mu_y \, P(w,y,t) \big]}{\partial y} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(w,y,t) \big]}{\partial y^2}  \, -\frac{\partial \big[\mu_w \, P(w,y,t) \big]}{\partial w} \, +\, .5 \, \frac{\partial^2 \big[{\sigma_w}^2 \, P(w,y,t) \big]}{\partial w^2} \, \Big] \, \, dy \, dw   \, [$]
[$]+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(w_b(n),y,t) dy [$] 
[$]+ \int_0^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t} \, P(w,y_b,t) dw =0 [$]

We do the integrations in the first term to get
[$]=\int_0^{w_b(n)}  \, \Big[\, -\mu_y \, P(w,y,t)  \, +\, .5 \, \frac{\partial \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(w,y,t) \big]}{\partial y}  \, \Big] dw [$]
[$]\int_0^{y_b(n,m)} \Big[ -\mu_w \, P(w,y,t)  \, +\, .5 \, \frac{\partial \big[{\sigma_w}^2 \, P(w,y,t) \big]}{\partial w} \, \Big] \, \, dy   \, [$]
[$]+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(w_b(n),y,t) dy [$] 
[$]+ \int_0^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t} \, P(w,y_b,t) dw =0 [$]

Now we convert our coordinates to [$]Z_1[$] and [$]Z_2[$] instead of w and y as

[$]=\int_0^{w_b(n)}  \, \Big[\, -\mu_y \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}  \, +\, .5 \, \frac{\partial \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial y}  \, \Big] dw [$]
[$]\int_0^{y_b(n,m)} \Big[ -\mu_w \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}  \, +\, .5 \, \frac{\partial \big[{\sigma_w}^2 \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial w} \, \Big] \, \, dy   \, [$]
[$]+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(Z_1(n)) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} dy [$] 
[$]+ \int_0^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1(n)}{\partial w} \, \frac{\partial Z_2(m)}{\partial y} dw =0 [$]

which can be further simplified as

[$]=\int_0^{w_b(n)}  \, \Big[\, -\mu_y \,  \frac{\partial Z_2}{\partial y}  \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V}   \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big]  \, \Big] dZ_1 [$]
[$]\int_0^{y_b(n,m)} \Big[ -\mu_w  \, \frac{\partial Z_1}{\partial w}   \, +\, .5 \, {\sigma_w}^2   \, \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2 - {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \big]  \Big] \, \, dZ_2   \, [$]
[$]+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \,  \frac{\partial Z_1}{\partial w} \,  dZ_2 [$] 
[$]+ \int_0^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t}  \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0 [$]

But we notice that second and third term together are merely the solution of 1D independent PDE integrated in an extra dimension and since the original solution of both terms for independent 1D PDE was zero, the integrated term should also be zero. (Even if it were not zero, we could have found this solution from solving the 1D independent PDE), therefore
[$]\int_0^{y_b(n,m)} \Big[ -\mu_w  \, \frac{\partial Z_1}{\partial w}   \, +\, .5 \, {\sigma_w}^2   \, \big[ -Z_1 \, {(\frac{\partial Z_1}{\partial w})}^2 - {(\frac{\partial Z_1}{\partial w})}^3 \frac{\partial^2 w}{\partial {Z_1}^2} \big]  \Big] \, \, dZ_2   \,[$]
[$]+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \,  \frac{\partial Z_1}{\partial w} \,  dZ_2=0 [$] 

Now we are left with the terms to solve as

[$]=\int_0^{w_b(n)}  \, \Big[\, -\mu_y \,  \frac{\partial Z_2}{\partial y}  \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V}   \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big]  \, \Big] dZ_1 [$]
[$]+ \int_0^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t}  \, \frac{\partial Z_2(m)}{\partial y} dZ_1 =0 [$]

For a given value of m, we can recursively solve the above PDE for each value of n. To make matters simple, I write the above equation as

[$]=\int_0^{w_b(n)}  \, \Big[\, -\mu_y \,  \frac{\partial Z_2}{\partial y}  \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V}   \, \big[ -Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 - {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big]  \, \Big] dZ_1 [$]
[$]+ \int_0^{w_b(n-1)}  \frac{\partial y_b(n,m)}{\partial t}  \, \frac{\partial Z_2}{\partial y} dZ_1 [$]
[$]+ \int_{w_b(n-1)}^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t}  \, \frac{\partial Z_2}{\partial y} dZ_1 =0[$]

where in the above equation, I have expanded the last terms as sum of all previous terms and one current term to solve and now the equation can be solved as

[$] \frac{\partial y_b(n,m)}{\partial t} =\Bigg[ \int_0^{w_b(n)}  \, \Big[\, \mu_y \,  \frac{\partial Z_2}{\partial y}  \, +\, .5 \, {\sigma_y}^2 \, {V(w)}^{2 \gamma_V}   \, \big[ Z_2 \, {(\frac{\partial Z_2}{\partial y})}^2 + {(\frac{\partial Z_2}{\partial y})}^3 \frac{\partial^2 y}{\partial {Z_2}^2} \big]  \, \Big] dZ_1 [$]
[$]- \int_0^{w_b(n-1)}  \frac{\partial y_b(n,m)}{\partial t}  \, \frac{\partial Z_2}{\partial y} dZ_1 \Bigg] / \Big[ \int_{w_b(n-1)}^{w_b(n)}  \, \frac{\partial Z_2}{\partial y} dZ_1 \Big][$]

please note that [$]\frac{\partial y_b(n,m)}{\partial t} \Delta t[$] is the change in coordinates y(n,m) so as to conserve the mass in a small time step according to given PDE dynamics.
The last equation solves the 2D PDE in the asset dimension. Adding correlation term should not be difficult but then we have to make a few changes as correlation will appear in jacobian. Also derivatives of asset with respect to [$]Z_1[$] would have to be included something that was making my equation unstable when I tried it last time but I will try it again to see if I could make it work.
I will be posting the program for uncorrelated PDE based on these analytics tomorrow.
.
.
Friends, in the term
Now we convert our coordinates to [$]Z_1[$] and [$]Z_2[$] instead of w and y as

[$]=\int_0^{w_b(n)}  \, \Big[\, -\mu_y \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}  \, +\, .5 \, \frac{\partial \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial y}  \, \Big] dw [$]
[$]\int_0^{y_b(n,m)} \Big[ -\mu_w \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y}  \, +\, .5 \, \frac{\partial \big[{\sigma_w}^2 \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial w} \, \Big] \, \, dy   \, [$]
[$]+\frac{\partial w_b}{\partial t} \int_0^{y_b(n,m)} \, P(Z_1(n)) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} dy [$] 
[$]+ \int_0^{w_b(n)}  \frac{\partial y_b(n,m)}{\partial t} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1(n)}{\partial w} \, \frac{\partial Z_2(m)}{\partial y} dw =0 [$]

I have not properly handled the derivative  given by expression  [$]\frac{\partial \big[{\sigma_y}^2 \, {V(w)}^{2 \gamma_V} \, P(Z_1) \, P(Z_2) \, \frac{\partial Z_1}{\partial w} \, \frac{\partial Z_2}{\partial y} \big]}{\partial y} [$]
y is changing in both [$]Z_1[$] and [$]Z_2[$] directions. Also [$]{V(w)}^{2 \gamma_V}[$] also takes a derivative with respect to y through inverse derivative of y to w.  I wanted to fix these things tonight but I have been given mind control chemicals in food and I find it impossible to work tonight. I am in a very poor state mentally and totally unable to work. I also received an antipsychotic injection today in my arm. I will try fixing the errors in next 2-3 days and then post the program here. Please protest against my mind control. It is very painful to be unable to work when I need to work hard to solve the problem. 
.
.
Friends, I have been working on the model. The above derivatives in [$]Z_1[$] directions are not relevant for uncorrelated problem but they are very relevant when specify a  positive or negative correlation between SV and asset. I have worked with all these correlated terms and related derivatives and now working to complete the model with correlations. When we are dealing with correlation coordinates, the nature of problem is very similar to as if we were in original non-bessel coordinates
and we have to add (apart from other terms) an extra drift like term that appeared in solution with original coordinates whose sigh now depends on sign of correlation. I hope to post the ( worked out but still preliminary)  program with correlations in another two days possibly on monday.
 
User avatar
Amin
Topic Author
Posts: 2853
Joined: July 14th, 2002, 3:00 am

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

September 6th, 2021, 6:34 pm

Though I have settled most of the issues with correlated stochastic volatility, there are still points that need a lot of work and remain to be understood so it would be at least 3-4 days before I could post a properly working correlated SV model. Since I have made a lot of changes to non-correlated version, I will be posting a properly commented version of non-correlated SV model tomorrow.
 
User avatar
Amin
Topic Author
Posts: 2853
Joined: July 14th, 2002, 3:00 am

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

September 7th, 2021, 6:26 pm

Friends, I meant to write comments on the program but rather continued to work on some stability issues. Earlier versions were working with moderate volatilities only and if we increased volatility scaling factor of the asset sigma1 greater than .65, the program would not work. Today I worked on these issues and now you can run SV pricing scenarios with high volatility and scaling factor sigma1 equal to one. It is a far more stable program and can also take mildly negative values in the asset. Because I continued to work on stability issues, I did not get time to write comments. I will post the commented program tomorrow. 
 
User avatar
Amin
Topic Author
Posts: 2853
Joined: July 14th, 2002, 3:00 am

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

September 8th, 2021, 3:30 pm

Friends, I continued to work to improve the behavior of program at zero and made some more improvements. When I was working on the program in my room, mind control agents continued to repeatedly release very very sickening gas in my room and now I have very little energy left to do any more work. The gas continued to sap all energy out of me. During past few weeks, mind control agents have got far bolder in releasing the gas in my room ,out of the ACs in every room in the house. Several rooms in our house had better air-conditioning but after I would work there for a few days, the would start releasing large amounts of gas from the ACs. Now I do not turn on air-conditioning despite very high humidity but they continue to release large amount of gas in my room and other places where I work. 
I continued to write for a week to ten days about mind control torture on me but the torture was not decreasing so I just simply gave up on writing about it everyday.
As a succession of what I revealed a few weeks ago, I will like to tell friends that a jewish billionaire is Godfather of mind control. Though there are far richer billionaires in US who are staunchly opposed to mind control, they are not ready to shell out hundreds of millions of dollars to bribe anyone to end mind control while mind control Godfather  continues to give large bribes to people to continue mind control. Many of the people related to defense in mind control are given extremely high paying jobs in Godfather's companies after they retire from mind control and Godfather is extremely generous towards them and it is well established and well known to people in mind control after retarding muslims, blacks and others in their career, they will move to paradise after they retire. No wonder there is no let up in mind control whatever good people try. So much so for goodness in American society where at the end of the day it is money that controls everything.
 
User avatar
Amin
Topic Author
Posts: 2853
Joined: July 14th, 2002, 3:00 am

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

September 9th, 2021, 4:44 pm

Friends here is the new improved program with comments. It needs many further improvements but still lot better than the earlier version.
.
function [] = SVSDEAnalyticAndMonteCarlo2DWilmott40()

%I have improved the program and if some of the asset goes to zero in a
%small fashion, it works but if volatility goes to zero, it does not work.
%There are still many improvements needed. It works fine with high
%stochastic volatility if there is very small mass at zero.
%Sometimes for high volatilities graph do not remain smooth and that has to
%be solved by an additional interpolation in next version.

%If parameters are extreme and the analytic distribution is smooth but has
%mismatch with monte carlo, try decreasing the step size of analytic
%simulation and monte carlo simulation.

% dX(t)=(a X(t)^alpha1 + b X(t)^alpha2) dt + sqrt(1-rho^2) sigma1 V(t)^gammaV X(t)^gammaX dZ1(t) + rho sigma1 V(t)^gammaV X(t)^gammaX dZ2(t)
% dV(t)=(mu1 V(t)^beta1 + mu2 V(t)^beta2)*dt + sigma0 V(t)^gamma dZ2(t)


%%%%%%%%%%%%%
X0=1.0;
alpha1=0;
alpha2=1;
kappaX=0.0;
thetaX=01.00;
a=kappaX*thetaX;
b=-kappaX;
rho=-.7*0;
sigma1=.85;
gammaV=.5;
gammaX=.95;
%%%%%%%%%%%%%
V0=.25;
thetaV=.075;
kappaV=2.0;
mu1=kappaV*thetaV;
mu2=-kappaV;
beta1=0;
beta2=1;
gamma=.95;
sigma0=.950;
%%%%%%%%%%%%

dtM=.03125;
TtM=32;



dt=dtM/32;
Tt=TtM*32;
Tt=256*4;
dtM=dt*16;
TtM=Tt/16;
T=TtM*dtM;
paths=100000;



dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=40;  % No of normal density subdivisions
NnMidl=20;%One half density Subdivision left from mid of normal density(low)
NnMidh=21;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;


dMm=.25/1;   % Normal density subdivisions width. would change with number of subdivisions
Mm=32;  % No of normal density subdivisions
MmMidl=16;%One half density Subdivision left from mid of normal density(low)
MmMidh=17;%One half density subdivision right from the mid of normal density(high)
MmMid=4.0;

  
w(1:Nn)=V0^(1-gamma)/(1-gamma); %Stochastic Volatility in bessel coordinates.
y(1:Nn,1:Mm)=X0^(1-gammaX)/(1-gammaX);%Asset in bessel coordinates

Z(1:Nn)=(((1:Nn)-.5)*dNn-NnMid);   %Z is standard gaussian associated with Stochastic volatility.


Z0(1:Mm)=(((1:Mm)-.5)*dMm-MmMid); %Z0 is standard gaussian associated with asset given a certain level of stochastic volatility.
Z
Z0
str=input('Look at Z');

%Belwo ZProb is integrated probability within each Z grid cell for
%stochastic volatility.
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);

%ZCProb(1:Nn)=normcdf(Z(1:Nn));
%ZProb0(1)=normcdf(Z(1),0,1);
%ZProb0(2:Nn)=normcdf(Z(2:Nn),0,1)-normcdf(Z(1:Nn-1),0,1);

%Belwo Z0Prob is integrated probability within each Z0 grid cell for
%asset given a certain level of stochastic volatility.
Z0Prob(1)=normcdf(.5*Z0(1)+.5*Z0(2),0,1)-normcdf(.5*Z0(1)+.5*Z0(2)-dMm,0,1);
Z0Prob(Mm)=normcdf(.5*Z0(Mm)+.5*Z0(Mm-1)+dMm,0,1)-normcdf(.5*Z0(Mm)+.5*Z0(Mm-1),0,1);
Z0Prob(2:Mm-1)=normcdf(.5*Z0(2:Mm-1)+.5*Z0(3:Mm),0,1)-normcdf(.5*Z0(2:Mm-1)+.5*Z0(1:Mm-2),0,1);

%Z0Prob0(1)=normcdf(Z0(1),0,1);
%Z0Prob0(2:Mm)=normcdf(Z0(2:Mm),0,1)-normcdf(Z0(1:Mm-1),0,1);


%Z0CProb(1:Mm)=normcdf(Z0(1:Mm));


%%%%%%%%%%%%%%%%%%%%%%%%5
%Below wnStart is starting grid cell for volatility. This remains one
%throughout the implementation something that will change in next version
%where we will tackle stochastic volatility going to zero.

wnStart=1;

%Below ynStart is starting grid cell for asset. This remains one
%throughout the implementation and is used where it is difficult to apply
%volatility dependent variable initial grid point.
ynStart0=1;

%Below ynStart is starting grid cell for asset and this starting grid cell
%varies with volatility level especially when at high volatility levels,
%the initial asset diffusion grid points have reached zero so they have to
%be treated differently while there are low Stochastic volatility grids
%where diffusion has not reached zero and it remains one at those places.
ynStart(1:Nn)=1;
ynStart(wnStart:Nn)=1;
dwdZ(wnStart:Nn)=0.0;
tic
for tt=1:Tt
%Below we calculate volatility and asset in original coordinates and then
%calculate wMu0dt which is drift in SV in bessel coordinates and yMu0dt
%which is drift in asset in bessel coordinates.
    V(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
    X(wnStart:Nn,ynStart0:Mm)=((1-gammaX)*y(wnStart:Nn,ynStart0:Mm)).^(1/(1-gammaX));

for nn=wnStart:Nn
    %yMu0dt is asset drift in bessel coordinates.
yMu0dt(nn,ynStart0:Mm)= (a*X(nn,ynStart0:Mm).^(alpha1-gammaX)+b*X(nn,ynStart0:Mm).^(alpha2-gammaX) - ...
     .5*gammaX*sigma1^2*V(nn).^(2*gammaV).*X(nn,ynStart0:Mm).^(gammaX-1))*dt + ...
     (a*(alpha1-gammaX).*X(nn,ynStart0:Mm).^(alpha1-gammaX-1)+b*(alpha2-gammaX)*X(nn,ynStart0:Mm).^(alpha2-gammaX-1) - ...
     .5*gammaX*sigma1^2*V(nn).^(2*gammaV).*(gammaX-1)*X(nn,ynStart0:Mm).^(gammaX-2)).* ...
     (a*X(nn,ynStart0:Mm).^(alpha1)+b*X(nn,ynStart0:Mm).^(alpha2))*dt^2/2+ ...
     -(.5*gammaX*sigma1^2*(2*gammaV).*V(nn).^(2*gammaV-1).*X(nn,ynStart0:Mm).^(gammaX-1)).* ...
     (mu1*V(nn).^(beta1)+mu2*V(nn).^(beta2))*dt^2/2 + ...
    (a*(alpha1-gammaX).*(alpha1-gammaX-1).*X(nn,ynStart0:Mm).^(alpha1-gammaX-2)+b*(alpha2-gammaX-1).*(alpha2-gammaX)*X(nn,ynStart0:Mm).^(alpha2-gammaX-2) - ...
    .5*gammaX*sigma1^2*(gammaX-1)*(gammaX-2)*V(nn).^(2*gammaV).*X(nn,ynStart0:Mm).^(gammaX-3)).* ...
    .5*sigma1^2*V(nn).^(2*gammaV).*X(nn,ynStart0:Mm).^(2*gammaX)*dt^2/2+ ...
    -(.5*gammaX*sigma1^2*(2*gammaV).*(2*gammaV-1).*V(nn).^(2*gammaV-2).*X(nn,ynStart0:Mm).^(gammaX-1)).* ...
     (sigma0^2*V(nn).^(2*gamma))*dt^2/2;
 %Below is Z0 dependent expanded volatility whose square needs to be used instead of sigma^2*dt;
    c1y(nn,ynStart0:Mm)=sigma1*V(nn).^(gammaV).*sqrt(dt)+ ...
        (a*(alpha1-gammaX).*X(nn,ynStart0:Mm).^(alpha1-gammaX-1)+b*(alpha2-gammaX)*X(nn,ynStart0:Mm).^(alpha2-gammaX-1) - ...
        .5*gammaX*sigma1^2*V(nn).^(2*gammaV).*(gammaX-1)*X(nn,ynStart0:Mm).^(gammaX-2)).* ...
        sigma1*V(nn).^(gammaV).*X(nn,ynStart0:Mm).^(gammaX).*(1-1/sqrt(3)).*dt^1.5;
 %Below is Z (and not Z0) dependent volatility. Its square is not O(dt) but higher order but I believe it has to be added 
 %but more on it in next version of the program. It makes a very small difference in shape of density.    
    c2y(nn,ynStart0:Mm)=-(.5*gammaX*sigma1^2*(2*gammaV).*V(nn).^(2*gammaV-1).*X(nn,ynStart0:Mm).^(gammaX-1)).* ...
        sigma0*V(nn).^gamma .*(1-1/sqrt(3)).*dt^1.5;
 
end

  
  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Below we calculate drift of stochastic volatility in bessel coordinates.
    wMu0dt(wnStart:Nn)= (mu1*V(wnStart:Nn).^(beta1-gamma)+mu2*V(wnStart:Nn).^(beta2-gamma) - ...
    .5*gamma*sigma0^2*V(wnStart:Nn).^(gamma-1))*dt + ...
    (mu1*(beta1-gamma).*V(wnStart:Nn).^(beta1-gamma-1)+mu2*(beta2-gamma)*V(wnStart:Nn).^(beta2-gamma-1) - ...
    .5*gamma*sigma0^2*(gamma-1)*V(wnStart:Nn).^(gamma-2)).* ...
    (mu1*V(wnStart:Nn).^(beta1-gamma)+mu2*V(wnStart:Nn).^(beta2-gamma))*dt^2/2+ ...
    (mu1*(beta1-gamma).*(beta1-gamma-1).*V(wnStart:Nn).^(beta1-gamma-2)+mu2*(beta2-gamma-1).*(beta2-gamma)*V(wnStart:Nn).^(beta2-gamma-2) - ...
    .5*gamma*sigma0^2*(gamma-1)*(gamma-2).*V(wnStart:Nn).^(gamma-3)).* ...
    .5.*sigma0^2.*V(wnStart:Nn).^(2*gamma)*dt^2/2;
%Below we calculate volatility of stochastic volatility in bessel
%coordinates. Its square will be used in place of sigm0^2 dt.
    c1(wnStart:Nn)=sigma0.*sqrt(dt) + ...
        (mu1.*(beta1-gamma).*V(wnStart:Nn).^(beta1-gamma-1) + ...    
        mu2.*(beta2-gamma).*V(wnStart:Nn).^(beta2-gamma-1) - ...
        .5*gamma* sigma0^2 .*(gamma-1).*V(wnStart:Nn).^(gamma-2)) .* sigma0.* V(wnStart:Nn).^(gamma).*(1-1/sqrt(3)) .*dt^1.5;
    
    %Below is the proper way to add volatility innovation in existing
    %density.
    %It has two parts a:squared addition  and b:diffusion smoothing. 
    %Below is squared addition part.
    %wMid is point on w density that is associated with Z=0 or median point
    %of the density.
    [wMid] = InterpolateOrderN8(8,0,Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidl),Z(NnMidh),Z(NnMidh+1),Z(NnMidh+2),Z(NnMidh+3),w(NnMidl-3),w(NnMidl-2),w(NnMidl-1),w(NnMidl),w(NnMidh),w(NnMidh+1),w(NnMidh+2),w(NnMidh+3));
        
    Zt1(wnStart:Nn)=w(wnStart:Nn)-wMid;
        
    [dwdZ,d2wdZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt1,Z);
        
    C0(wnStart:Nn)=Zt1(wnStart:Nn)-dwdZ(wnStart:Nn).*Z(wnStart:Nn);    
    
    Zt2(wnStart:Nn)=C0(wnStart:Nn)+sign((dwdZ(wnStart:Nn))+c1(wnStart:Nn)).* ...
        sqrt(abs(sign((dwdZ(wnStart:Nn))).*(dwdZ(wnStart:Nn)).^2+ ...
        sign(c1(wnStart:Nn)).*c1(wnStart:Nn).^2)).*Z(wnStart:Nn);
    
    [dw2dZ,d2w2dZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt2,Z);
    %below the line ".5*c1(wnStart:Nn).^2.*(dw2dZ(wnStart:Nn).^(-2)).*d2w2dZ2A(wnStart:Nn);
    %is equivalent to diffusion smoothing part.
    %SO three terms below in dwb are contributions from drift, squared addition and
    %diffusion smoothing. dwb is net change in displacement of nth grid
    %point.
    if(tt>20)
        
        dwb(wnStart:Nn)= wMu0dt(wnStart:Nn)+(Zt2(wnStart:Nn)-Zt1(wnStart:Nn))+ ...
            .5*c1(wnStart:Nn).^2.*(dw2dZ(wnStart:Nn).^(-2)).*d2w2dZ2A(wnStart:Nn);
    else

        dwb(wnStart:Nn)=wMu0dt(wnStart:Nn)+Zt2(wnStart:Nn)-Zt1(wnStart:Nn);%+ ... 
    end
    dwb(isnan(dwb)==1)=0.0;
    dwb(isinf(dwb)==1)=0.0;
    %dwb(dwb>.5)=.50;
    w(wnStart:Nn)=w(wnStart:Nn)+dwb(wnStart:Nn);
    
    
  
  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%555
    %dwdZPrev(wnStart:Nn)=dwdZ(wnStart:Nn); 
    %%d2wdZ2Prev(wnStart:Nn)=d2wdZ2(wnStart:Nn);
   % 
   % [dwdZNew,d2wdZ2New] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,w,Z); 
  
tt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%Below we repeat(squared addition and diffusion smoothing) for every volatility conditional value of asset
%what we did for volatility itself in one dimension.
for nn=wnStart:Nn
    
    
    [yMid(nn)] = InterpolateOrderN8(8,0,Z0(MmMidl-3),Z0(MmMidl-2),Z0(MmMidl-1),Z0(MmMidl),Z0(MmMidh),Z0(MmMidh+1),Z0(MmMidh+2),Z0(MmMidh+3),y(nn,MmMidl-3),y(nn,MmMidl-2),y(nn,MmMidl-1),y(nn,MmMidl),y(nn,MmMidh),y(nn,MmMidh+1),y(nn,MmMidh+2),y(nn,MmMidh+3));
    Zy1(nn,ynStart0:Mm)=y(nn,ynStart0:Mm)-yMid(nn);
    Zy10(ynStart0:Mm)= Zy1(nn,ynStart0:Mm);   
    [dydZ0A,d2ydZ02A] = First2Derivatives2ndOrderEqSpacedA(ynStart(nn),Mm,dMm,Zy10,Z0);
    if(tt>1)
    dydZ0Prev(nn,ynStart(nn):Mm)=dydZ0(nn,ynStart(nn):Mm);
    end
    dydZ0(nn,ynStart(nn):Mm)=dydZ0A(ynStart(nn):Mm);    
    Cy0(nn,ynStart(nn):Mm)=Zy1(nn,ynStart(nn):Mm)-dydZ0A(ynStart(nn):Mm).*Z0(ynStart(nn):Mm);
        
    %Zt2(wnStart:Nn)=C0(wnStart:Nn)+abs(sqrt((dwdZ(wnStart:Nn)).^2+sigma0^2*dt)).*Z(wnStart:Nn);
%    Zy2(nn,ynStart:Mm)=Cy0(nn,ynStart:Mm)+abs(sqrt((dydZ0A(ynStart:Mm)).^2+c1y(nn,ynStart:Mm).^2)).*Z0(ynStart:Mm);
    
    
    Zy2(nn,ynStart(nn):Mm)=Cy0(nn,ynStart(nn):Mm)+sign((dydZ0A(ynStart(nn):Mm))+c1y(nn,ynStart(nn):Mm)).* ...
        sqrt(abs(sign((dydZ0A(ynStart(nn):Mm))).*(dydZ0A(ynStart(nn):Mm)).^2+ ...
        sign(c1y(nn,ynStart(nn):Mm)).*c1y(nn,ynStart(nn):Mm).^2)).*Z0(ynStart(nn):Mm);

    
    
    
    Zy20(ynStart(nn):Mm)= Zy2(nn,ynStart(nn):Mm);
    [dy2dZ0B,d2y2dZ02B] = First2Derivatives2ndOrderEqSpacedA(ynStart(nn),Mm,dMm,Zy20,Z0);
    dy2dZ0(nn,ynStart(nn):Mm)=dy2dZ0B(ynStart(nn):Mm);
    dy2dZ0Inv(nn,ynStart(nn):Mm)=1.0 ./dy2dZ0(nn,ynStart(nn):Mm);
    dy2dZ0Inv(isnan(dy2dZ0Inv)==1)=0;
    
    d2y2dZ02(nn,ynStart(nn):Mm)=d2y2dZ02B(ynStart(nn):Mm);
    d2y2dZ02(isnan(d2y2dZ02)==1)=0.0;
    dydZ0(nn,ynStart(nn):Mm)=dy2dZ0(ynStart(nn):Mm);
    if(tt>10)
    %yDiff0(nn,ynStart:Mm)=Zy2(nn,ynStart:Mm)-Zy1(nn,ynStart:Mm) + ...
    %    .5*sigma1^2* V(nn).^(2*gammaV).*dt.*(dy2dZ0(nn,ynStart:Mm).^(-2)).*d2y2dZ02(nn,ynStart:Mm);
    yDiff0(nn,ynStart(nn):Mm)= .5*c1y(nn,ynStart(nn):Mm).^2.*(dy2dZ0Inv(nn,ynStart(nn):Mm).^(2)).*d2y2dZ02(nn,ynStart(nn):Mm);
    else
        yDiff0(nn,ynStart(nn):Mm)=0.0;%Zy2(nn,ynStart:Mm)-Zy1(nn,ynStart:Mm);
    end
    yDiff0(isnan(yDiff0)==1)=0.0;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%  for mm=ynStart0:Mm
% %     
% %     c1y2A(wnStart:Nn)=c1y(wnStart:Nn,mm).^2;
% %     [dc1y2AdZ,d2c1y2AdZ2] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,c1y2A,Z);
% %     dc1y2dZ(wnStart:Nn,mm)=dc1y2AdZ(wnStart:Nn);
% % 
%      [ywMid(mm)] = InterpolateOrderN8(8,0,Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidl),Z(NnMidh),Z(NnMidh+1),Z(NnMidh+2),Z(NnMidh+3),y(NnMidl-3,mm),y(NnMidl-2,mm),y(NnMidl-1,mm),y(NnMidl,mm),y(NnMidh,mm),y(NnMidh+1,mm),y(NnMidh+2,mm),y(NnMidh+3,mm));
% %         
%      Zyw1(wnStart:Nn,mm)=y(wnStart:Nn,mm)-ywMid(mm);
%      Zyw10(wnStart:Nn)= Zyw1(wnStart:Nn,mm);       
%      if(tt>1)
%         dywdZPrev(wnStart:Nn,mm)=dywdZ(wnStart:Nn,mm);
%      end
%      [dywdZA,d2ywdZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zyw10,Z);
% %     
%      dywdZ(wnStart:Nn,mm)=dywdZA(wnStart:Nn);
% %     d2ywdZ2(wnStart:Nn,mm)=d2ywdZ2A(wnStart:Nn);
%  end




%In order to make interpolation easier, we sum together the drift, squared
%addition and diffusion smoothing in one term instead of three separate
%terms. Multiplication by .*dydZ0(nn,ynStart(nn):Mm).^-1 is in the
%algorithm.

for nn=wnStart:Nn
    for mm=wnStart:Mm
        if(y(nn,mm)>0)
            dy0(nn,ynStart(nn):Mm)=(yMu0dt(nn,ynStart(nn):Mm)+ Zy2(nn,ynStart(nn):Mm)-Zy1(nn,ynStart(nn):Mm)+yDiff0(nn,ynStart(nn):Mm) ...
                ).*dydZ0(nn,ynStart(nn):Mm).^-1;
        else
            dy0(nn,ynStart(nn):Mm)=real((yMu0dt(nn,ynStart(nn):Mm)+ Zy2(nn,ynStart(nn):Mm)-Zy1(nn,ynStart(nn):Mm))).*real(dydZ0(nn,ynStart(nn):Mm)).^-1; 
        end     
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

%All the values below with capital I in them are interpolated values. Two additional grids
%are interpolated between any two existing grids so that n=1:Nn is 3*n-2:3*Nn-2 on interpolated grid. 
%values with a like dy0Ia are starting values at start of grid cell.
%values withhout a and b like dy0I are middle values at center of grid cell.
%values with a like dy0Ib are values at end of grid cell.
%ZI are interpolated Z values at middle of grid cells.
%NnI = 3*Nn +2;
%dNnI = dNn/3.
%ZIProb is integrated probability in each interpolated grid cell.
[dy0I,dy0Ia,dy0Ib,ZI,NnI,dNnI,ZIProb] = InterpolateAssetAlongVolatilityWithEnds(dy0,Z,wnStart,Nn,dNn,ynStart0,Mm);
[dydZ0I,dydZ0Ia,dydZ0Ib,ZI,NnI,dNnI,ZIProb] = InterpolateAssetAlongVolatilityWithEnds(dydZ0,Z,wnStart,Nn,dNn,ynStart0,Mm);


if(tt>10)
 
%copying from above.
%%All the values with capital I in them are interpolated values. Two additional grids
%are interpolated between any two existing grids so that n=1:Nn is 3*n-2:3*Nn-2 on interpolated grid.
%Below is the probbility mass conservation algorithm on interpolated stochastic volatility grid.
    cc1=0;
    cc2=.5;
    cc3=.5;
    Tdy0(1:NnI,1:Mm)=0.0;
    for nn=wnStart:Nn
        ynStartI(3*nn-2)=ynStart(nn);
        ynStartI(3*nn-1)=ynStart(nn);
        ynStartI(3*nn)=ynStart(nn);
    end



    for nn=wnStart:NnI

        if(nn==wnStart)
            Tdy0(nn,ynStartI(nn):Mm)=(cc1*dy0Ia(nn,ynStartI(nn):Mm)+cc2*dy0I(nn,ynStartI(nn):Mm)+cc3*dy0Ib(nn,ynStartI(nn):Mm)).*ZIProb(nn);

            dyb(nn,ynStartI(nn):Mm)= (Tdy0(nn,ynStartI(nn):Mm))./(cc1*dydZ0Ia(nn,ynStartI(nn):Mm).^(-1).*ZIProb(nn)+ ...
                cc2*dydZ0I(nn,ynStartI(nn):Mm).^(-1).*ZIProb(nn)+cc3*dydZ0Ib(nn,ynStartI(nn):Mm).^(-1).*ZIProb(nn));
        
            Tyb(nn,ynStartI(nn):Mm)=dyb(nn,ynStartI(nn):Mm).*(cc1*dydZ0Ia(nn,ynStartI(nn):Mm).^(-1) + ...
                cc2*dydZ0I(nn,ynStartI(nn):Mm).^(-1)+cc3*dydZ0Ib(nn,ynStartI(nn):Mm).^(-1)).*ZIProb(nn); 


        else
   
            Tdy0(nn,ynStartI(nn):Mm)=Tdy0(nn-1,ynStartI(nn):Mm)+ ...
                cc1*dy0Ia(nn,ynStartI(nn):Mm).*ZIProb(nn)+ ...
                cc2*dy0I(nn,ynStartI(nn):Mm).*ZIProb(nn)+...
                cc3*dy0Ib(nn,ynStartI(nn):Mm).*ZIProb(nn);
            dyb(nn,ynStartI(nn):Mm)= (Tdy0(nn,ynStartI(nn):Mm)-Tyb(nn-1,ynStartI(nn):Mm))./ ...
                (cc1*dydZ0Ia(nn,ynStartI(nn):Mm).^(-1).*ZIProb(nn)+cc2*dydZ0I(nn,ynStartI(nn):Mm).^(-1).*ZIProb(nn)+ ...
                +cc3*dydZ0Ib(nn,ynStartI(nn):Mm).^(-1).*ZIProb(nn));
            Tyb(nn,ynStartI(nn):Mm)=Tyb(nn-1,ynStartI(nn):Mm)+ ...
                cc1*dyb(nn,ynStartI(nn):Mm).*dydZ0Ia(nn,ynStartI(nn):Mm).^(-1).*ZIProb(nn) + ... 
                cc2*dyb(nn,ynStartI(nn):Mm).*dydZ0I(nn,ynStartI(nn):Mm).^(-1).*ZIProb(nn)+ ...
                cc3*dyb(nn,ynStartI(nn):Mm).*dydZ0Ib(nn,ynStartI(nn):Mm).^(-1).*ZIProb(nn);
        end
    
    end

    %dyb
    %size(dyb)
    %str=input('Look at numbers_2');
    dyb(isnan(dyb)==1)=0.0;

    for nn=wnStart:Nn
    
        for mm=ynStart(nn):Mm
            
            ytemp(nn,mm)=y(nn,mm)+dyb(3*nn-2,mm);
                if(real(ytemp(nn,mm))<0)
                    ytemp(nn,mm)=0;
                end
                if((ytemp(nn,mm)<=0) && (mm>=ynStart(nn)) && (mm <=4)    )
                    ynStart(nn)=mm;
                end
        end
   
        y(nn,1:Mm)=0.0; 
        y(nn,ynStart(nn):Mm)=(ytemp(nn,ynStart(nn):Mm));
    end
            
    
    %y(real(y)<=0)

else

    
             y(wnStart:Nn,ynStart0:Mm)=y(wnStart:Nn,ynStart0:Mm)+ ...
             yMu0dt(wnStart:Nn,ynStart0:Mm)+yDiff0(wnStart:Nn,ynStart0:Mm)+ ...
            Zy2(wnStart:Nn,ynStart0:Mm)-Zy1(wnStart:Nn,ynStart0:Mm);% + ...
%            1*yDiff4(wnStart:Nn,ynStart:Mm);

end
   
 
end


Vw(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
DfVw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
    
    DfVw(nn) = (Vw(nn + 1) - Vw(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    %Change of variable derivative for densities
end

pVw(1:Nn)=0;
for nn = wnStart:Nn-1
    
   pVw(nn) = (normpdf(Z(nn),0, 1))/abs(DfVw(nn));
end


minX=10;
maxX=.001;
for mm=ynStart0:Mm
%Xy(mm)=sum(((1-gammaX).*y(wnStart:Nn,mm)).^(1/(1-gammaX)).*ZProb(wnStart:Nn));
    Xy_limits(mm)=sum(((1-gammaX).*y(37,mm)).^(1/(1-gammaX)));

end
y(real(y)<0)=0;

w
y
Xy_limits
if(Xy_limits(2)<minX)
    minX=real(Xy_limits(2));
end
if(Xy_limits(28)>maxX)
    maxX=real(Xy_limits(28));
end



minX=minX-.05;
maxX=maxX+1.250;
minX=max(.001,minX);





Mm00=300;
dX=(maxX-minX)/(Mm00-1)

minX
maxX
X5(1:Mm00)=minX+((1:Mm00)-1).*dX
%dX0(1:100-1)=dX
%dX0(100)=dX
%dX0(101:200)=dX
%X5(101:300)=X5(100)+(1:200).*dX
Mm0=Mm00;
%dX=(maxX-minX)/(Mm0-1);
X5
dX
str=input('Look at X');
%[X,XCdf] = ConstructCDFOnNewGrid02a(X5,y,wnStart,Nn,ynStart,Mm,Mm0,ZProb,Z0,gammaX)

[X,XCdf] = ConstructCDFOnNewGridWithInterpolation(X5,y,wnStart,Nn,dNn,Z,ynStart0,Mm,Mm0,Z0,gammaX);


Zx(2:Mm0-1)=norminv(XCdf(2:Mm0-1));%is calculated at boundaries of cells
%ZX(2:Mm0)=(Zb(2:Mm0)+Zb(1:Mm0-1))/2.0;%is calculated at center of grid cells.
%Zx(2:Mm0-1)=Zb(2:Mm0-1);%is calculated at center of grid cells.
[Zx(1)] = InterpolateOrderN6(4,X(1),X(2),X(3),X(4),X(5),0,0,Zx(2),Zx(3),Zx(4),Zx(5),0,0);
dXdZx(3:Mm0-2)=(X(4:Mm0-1)-X(2:Mm0-3))./(Zx(4:Mm0-1)-Zx(2:Mm0-3))
[dXdZx(2)] = InterpolateOrderN6(4,Zx(2),Zx(3),Zx(4),Zx(5),Zx(6),0,0,dXdZx(3),dXdZx(4),dXdZx(5),dXdZx(6),0,0);
[dXdZx(1)] = InterpolateOrderN6(4,Zx(1),Zx(3),Zx(4),Zx(5),Zx(6),0,0,dXdZx(3),dXdZx(4),dXdZx(5),dXdZx(6),0,0);

Zb(1:Mm0-3)=(Zx(1:Mm0-3)+Zx(2:Mm0-2))/2;

ZbProb(1)=normcdf(Zb(1));
ZbProb(2:Mm0-3)=normcdf(Zb(2:Mm0-3))-normcdf(Zb(1:Mm0-4));

for mm=1:Mm0-10    
    pX(mm) = (normpdf(Zx(mm),0, 1))/abs(dXdZx(mm));
end

toc 
%Zx
 
%pX
    pX(isnan(pX)==1)=0.0;
    pX(isinf(pX)==1)=0.0;
%pX
%XProb
XX=X;
plot(XX(1:Mm0-10),XCdf(1:Mm0-10),'r');
str=input('Look at graph');
plot(XX(2:Mm0-10),pX(2:Mm0-10),'r');
str=input('Look at graph');


ItoHermiteSVMean=sum(Vw(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1))
TrueSVMean=thetaV+(V0-thetaV)*exp(-kappaV*dt*Tt)


ItoHermiteAssetMeanX2=sum(X(2:Mm0-3).*ZbProb(2:Mm0-3))
TrueAssetMean=thetaX+(X0-thetaX)*exp(-kappaX*dt*Tt)



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%rng(29079137, 'twister')
rng(15898837, 'twister')
paths=150000;
V(1:paths)=V0;  %Original process monte carlo.
X=0.0;
X(1:paths)=X0;
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;
    

    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;


end

SVolMeanAnalytic=thetaV+(V0-thetaV)*exp(-kappaV*dt*Tt)
SVolMeanMC=sum(V(1:paths))/paths
AssetMeanAnalytic=thetaX+exp(-kappaX*T)*(X0-thetaX)
AssetMeanMC=sum(X(1:paths))/paths



MaxCutOff=30;
%NoOfBins=round(600*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too 
NoOfBins=750;
[XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti_NEW(X,paths,NoOfBins,MaxCutOff );
plot(IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g',XX(1:Mm0-10),pX(1:Mm0-10),'r');
%plot(IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'g',XX(3:80),pX(3:80),'r');
 %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,thetaX=%.3f,kappaX=%.2f,gammaX=%.3f,gammaV=%.3f,rho=%.2f,sigma1=%.2f,T=%.2f,dt=%.5f,AMean=%.4f,McMean=%.4f,TMean=%.4f', X0,thetaX,kappaX,gammaX,gammaV,rho,sigma1,T,dt,ItoHermiteAssetMeanX2,AssetMeanMC,AssetMeanAnalytic));%,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('Green line is Monte carlo density of Asset SDE while red line is analytic density.');

MaxCutOff=30;
%NoOfBins=round(600*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too 
NoOfBins=1000;
[VDensity,IndexOutV,IndexMaxV] = MakeDensityFromSimulation_Infiniti_NEW(V,paths,NoOfBins,MaxCutOff );
plot(IndexOutV(1:IndexMaxV),VDensity(1:IndexMaxV),'g',Vw(wnStart+1:Nn-1),pVw(wnStart+1:Nn-1),'r');
 %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');
%plot(v(wnStart+1:Nn-1),pv(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
 
title(sprintf('V0 = %.4f,thetaV=%.3f,kappaV=%.2f,gamma=%.3f,sigma0=%.2f,T=%.2f,dt=%.5f,SVolMean=%.4f,SVolMeanAnalytic=%.4f', V0,thetaV,kappaV,gamma,sigma0,T,dt,SVolMeanMC,SVolMeanAnalytic));%,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('Green line is Monte carlo density of volatility SDE while red line is analytic density.');

end

.
.
.
function [yI,yaI,ybI,ZI,NnI,dNnI,ZIProb] = InterpolateAssetAlongVolatilityWithEnds(y,Z,wnStart,Nn,dNn,ynStart,Mm)
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here

    
    x01=Z(wnStart)+1/3*dNn;
    x02=Z(wnStart)+2/3*dNn;
    
    x0a=Z(wnStart)-1/6*dNn;
    x0b=Z(wnStart)+1/6*dNn;
    x0c=Z(wnStart)+1/2*dNn;
    
    
    
    x1=Z(wnStart);
    x2=Z(wnStart+1);
    x3=Z(wnStart+2);
    x4=Z(wnStart+3);
    x5=0;
    x6=0;
 

    y1=y(wnStart,ynStart:Mm);
    y2=y(wnStart+1,ynStart:Mm);
    y3=y(wnStart+2,ynStart:Mm);
    y4=y(wnStart+3,ynStart:Mm);
    y5=0;
    y6=0;
    
    N=4;
    [y01] = InterpolateOrderN6(N,x01,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y02] = InterpolateOrderN6(N,x02,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);

    [y0a] = InterpolateOrderN6(N,x0a,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y0b] = InterpolateOrderN6(N,x0b,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y0c] = InterpolateOrderN6(N,x0c,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    

    
    yI(wnStart,ynStart:Mm)=y(wnStart,ynStart:Mm);
    yI(wnStart+1,ynStart:Mm)=y01(ynStart:Mm);
    yI(wnStart+2,ynStart:Mm)=y02(ynStart:Mm);
    
    yaI(wnStart,ynStart:Mm)=y0a(ynStart:Mm);
    yaI(wnStart+1,ynStart:Mm)=y0b(ynStart:Mm);
    yaI(wnStart+2,ynStart:Mm)=y0c(ynStart:Mm);
    
    
    
    x0a=Z(wnStart+1)-1/6*dNn;
    x0b=Z(wnStart+1)+1/6*dNn;
    x0c=Z(wnStart+1)+1/2*dNn;
    
    x01=Z(wnStart+1)+1/3*dNn;
    x02=Z(wnStart+1)+2/3*dNn;
    x1=Z(wnStart);
    x2=Z(wnStart+1);
    x3=Z(wnStart+2);
    x4=Z(wnStart+3);
    x5=Z(wnStart+4);
    x6=0;
    

    y1=y(wnStart,ynStart:Mm);
    y2=y(wnStart+1,ynStart:Mm);
    y3=y(wnStart+2,ynStart:Mm);
    y4=y(wnStart+3,ynStart:Mm);
    y5=y(wnStart+4,ynStart:Mm);
    y6=0;
    N=5;
    [y01] = InterpolateOrderN6(N,x01,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y02] = InterpolateOrderN6(N,x02,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    
    [y0a] = InterpolateOrderN6(N,x0a,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y0b] = InterpolateOrderN6(N,x0b,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y0c] = InterpolateOrderN6(N,x0c,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);

    yI(wnStart+3,ynStart:Mm)=y(wnStart+1,ynStart:Mm);
    yI(wnStart+4,ynStart:Mm)=y01(ynStart:Mm);
    yI(wnStart+5,ynStart:Mm)=y02(ynStart:Mm);
    
    yaI(wnStart+3,ynStart:Mm)=y0a(ynStart:Mm);
    yaI(wnStart+4,ynStart:Mm)=y0b(ynStart:Mm);
    yaI(wnStart+5,ynStart:Mm)=y0c(ynStart:Mm);
    

for nn=wnStart+2:Nn-3
    x01=Z(nn)+1/3*dNn;
    x02=Z(nn)+2/3*dNn;
    
    x0a=Z(nn)-1/6*dNn;
    x0b=Z(nn)+1/6*dNn;
    x0c=Z(nn)+1/2*dNn;
    
    
    x1=Z(nn-2);
    x2=Z(nn-1);
    x3=Z(nn);
    x4=Z(nn+1);
    x5=Z(nn+2);
    x6=Z(nn+3);

    y1=y(nn-2,ynStart:Mm);
    y2=y(nn-1,ynStart:Mm);
    y3=y(nn,ynStart:Mm);
    y4=y(nn+1,ynStart:Mm);
    y5=y(nn+2,ynStart:Mm);
    y6=y(nn+3,ynStart:Mm);
    
    N=6;
    
    [y01] = InterpolateOrderN6(N,x01,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y02] = InterpolateOrderN6(N,x02,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    
    [y0a] = InterpolateOrderN6(N,x0a,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y0b] = InterpolateOrderN6(N,x0b,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y0c] = InterpolateOrderN6(N,x0c,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    
    
    yI(3*nn-2,ynStart:Mm)=y(nn,ynStart:Mm);
    yI(3*nn-1,ynStart:Mm)=y01(ynStart:Mm);
    yI(3*nn,ynStart:Mm)=y02(ynStart:Mm);
    
    yaI(3*nn-2,ynStart:Mm)=y0a(ynStart:Mm);
    yaI(3*nn-1,ynStart:Mm)=y0b(ynStart:Mm);
    yaI(3*nn,ynStart:Mm)=y0c(ynStart:Mm);
    
    
end

    x01=Z(Nn-2)+1/3*dNn;
    x02=Z(Nn-2)+2/3*dNn;
    
    x0a=Z(Nn-2)-1/6*dNn;
    x0b=Z(Nn-2)+1/6*dNn;
    x0c=Z(Nn-2)+1/2*dNn;
    
    
    
    x1=Z(Nn-4);
    x2=Z(Nn-3);
    x3=Z(Nn-2);
    x4=Z(Nn-1);
    x5=Z(Nn);
    x6=0;
    

    y1=y(Nn-4,ynStart:Mm);
    y2=y(Nn-3,ynStart:Mm);
    y3=y(Nn-2,ynStart:Mm);
    y4=y(Nn-1,ynStart:Mm);
    y5=y(Nn,ynStart:Mm);
    y6=0;
    N=5;
    [y01] = InterpolateOrderN6(N,x01,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y02] = InterpolateOrderN6(N,x02,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    
    [y0a] = InterpolateOrderN6(N,x0a,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y0b] = InterpolateOrderN6(N,x0b,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y0c] = InterpolateOrderN6(N,x0c,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    

    yI(3*(Nn-2)-2,ynStart:Mm)=y(Nn-2,ynStart:Mm);
    yI(3*(Nn-2)-1,ynStart:Mm)=y01(ynStart:Mm);
    yI(3*(Nn-2),ynStart:Mm)=y02(ynStart:Mm);
    
    yaI(3*(Nn-2)-2,ynStart:Mm)=y0a(ynStart:Mm);
    yaI(3*(Nn-2)-1,ynStart:Mm)=y0b(ynStart:Mm);
    yaI(3*(Nn-2),ynStart:Mm)=y0c(ynStart:Mm);
    
    x01=Z(Nn-1)+1/3*dNn;
    x02=Z(Nn-1)+2/3*dNn;
    
    x0a=Z(Nn-1)-1/6*dNn;
    x0b=Z(Nn-1)+1/6*dNn;
    x0c=Z(Nn-1)+1/2*dNn;
    x0d=Z(Nn-1)+5/6*dNn;
    x0e=Z(Nn-1)+7/6*dNn;
    
    x1=Z(Nn-3);
    x2=Z(Nn-2);
    x3=Z(Nn-1);
    x4=Z(Nn);
    x5=0;
    x6=0;
    

    y1=y(Nn-3,ynStart:Mm);
    y2=y(Nn-2,ynStart:Mm);
    y3=y(Nn-1,ynStart:Mm);
    y4=y(Nn,ynStart:Mm);
    y5=0;
    y6=0;
    N=4;
    [y01] = InterpolateOrderN6(N,x01,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y02] = InterpolateOrderN6(N,x02,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    
    [y0a] = InterpolateOrderN6(N,x0a,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y0b] = InterpolateOrderN6(N,x0b,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y0c] = InterpolateOrderN6(N,x0c,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y0d] = InterpolateOrderN6(N,x0d,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y0e] = InterpolateOrderN6(N,x0e,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    

    yI(3*(Nn-1)-2,ynStart:Mm)=y(Nn-1,ynStart:Mm);
    yI(3*(Nn-1)-1,ynStart:Mm)=y01(ynStart:Mm);
    yI(3*(Nn-1),ynStart:Mm)=y02(ynStart:Mm);
    
    yaI(3*(Nn-1)-2,ynStart:Mm)=y0a(ynStart:Mm);
    yaI(3*(Nn-1)-1,ynStart:Mm)=y0b(ynStart:Mm);
    yaI(3*(Nn-1),ynStart:Mm)=y0c(ynStart:Mm);
    
    
    yaI(3*Nn-2,ynStart:Mm)=y0d(ynStart:Mm);
    ybI(3*Nn-2,ynStart:Mm)=y0e(ynStart:Mm);
    
    ybI(1:3*(Nn-1),ynStart:Mm)=yaI(2:3*Nn-2,ynStart:Mm);
    
    yI(3*Nn-2,ynStart:Mm)=y(Nn,ynStart:Mm);
    
    
    for nn=1:Nn-1
       ZI(nn*3-2)=Z(nn);
       ZI(nn*3-1)=Z(nn)+1/3*dNn;
       ZI(nn*3)=Z(nn)+2/3*dNn;
    end
    ZI(3*Nn-2)=Z(Nn);
    
    NnI=3*Nn-2;
    dNnI=1/3*dNn;
    
    ZIProb(wnStart)=normcdf((.5*ZI(wnStart)+.5*ZI(wnStart+1)));
    for nn=wnStart+1:NnI-1
        ZIProb(nn)=normcdf((.5*ZI(nn)+.5*ZI(nn+1)))-normcdf((.5*ZI(nn)+.5*ZI(nn-1)));
    end
    
    ZIProb(NnI)=1-normcdf((.5*ZI(NnI)+.5*ZI(NnI-1)));
    
    
    end
.
.
.
function [X,XCdf] = ConstructCDFOnNewGridWithInterpolation(X,y,wnStart,Nn,dNn,Z,ynStart,Mm,Mm0,Z1,gammaX)


Xyy0(wnStart:Nn,ynStart:Mm)=((1-gammaX).*y(wnStart:Nn,ynStart:Mm)).^(1/(1-gammaX));

[yI,ZI,NnI,dNnI,ZIProb] = InterpolateAssetAlongVolatility(Xyy0,Z,wnStart,Nn,dNn,ynStart,Mm);


Xyy(wnStart:NnI,ynStart:Mm)=yI(wnStart:NnI,ynStart:Mm);

XCdf(1:Mm0)=0.0;


for nn=wnStart:NnI
    Xy(ynStart:Mm)=real(Xyy(nn,ynStart:Mm));
    %str=input('Look at Xy');
    mm=2;
    mm0=1;
    while(mm0<Mm0-2)
        if((mm==2))
           while((X(mm0)<=Xy(mm)) &&(mm0<Mm0))
                Zx(mm0) =InterpolateOrderN6(4,X(mm0),Xy(mm),Xy(mm+1),Xy(mm+2),Xy(mm+3),0,0,Z1(mm),Z1(mm+1),Z1(mm+2),Z1(mm+3),0,0);
                XCdf(mm0)=XCdf(mm0)+normcdf(real(Zx(mm0))).*ZIProb(nn);
                mm0=mm0+1;
           end
           mm=mm+1;
        elseif( (mm==3))
            while((X(mm0)<=Xy(mm))&&(mm0<Mm0))
                Zx(mm0) =InterpolateOrderN6(5,X(mm0),Xy(mm-1),Xy(mm),Xy(mm+1),Xy(mm+2),Xy(mm+3),0,Z1(mm-1),Z1(mm),Z1(mm+1),Z1(mm+2),Z1(mm+3),0);
                if((mm0>1)&&(Zx(mm0)<Zx(mm0-1)))
                    Zx(mm0)=Zx(mm0-1);
                end
                XCdf(mm0)=XCdf(mm0)+normcdf(real(Zx(mm0))).*ZIProb(nn);
                mm0=mm0+1;
            end
            mm=mm+1;
        elseif( (mm==4))
            while((X(mm0)<=Xy(mm))&&(mm0<Mm0))
                Zx(mm0) =InterpolateOrderN6(5,X(mm0),Xy(mm-2),Xy(mm-1),Xy(mm),Xy(mm+1),Xy(mm+2),0,Z1(mm-2),Z1(mm-1),Z1(mm),Z1(mm+1),Z1(mm+2),0);
                if((mm0>1)&&(Zx(mm0)<Zx(mm0-1)))
                    Zx(mm0)=Zx(mm0-1);
                end
                XCdf(mm0)=XCdf(mm0)+normcdf(real(Zx(mm0))).*ZIProb(nn);
                mm0=mm0+1;
            end
            mm=mm+1;    
        elseif( (mm>4) && (mm<=Mm-3))
            while((X(mm0)<=Xy(mm))&&(mm0<Mm0))
                Zx(mm0) =InterpolateOrderN6(6,X(mm0),Xy(mm-3),Xy(mm-2),Xy(mm-1),Xy(mm),Xy(mm+1),Xy(mm+2),Z1(mm-3),Z1(mm-2),Z1(mm-1),Z1(mm),Z1(mm+1),Z1(mm+2));
                if((mm0>1)&&(Zx(mm0)<Zx(mm0-1)))
                    Zx(mm0)=Zx(mm0-1);
                end
                XCdf(mm0)=XCdf(mm0)+normcdf(real(Zx(mm0))).*ZIProb(nn);
                mm0=mm0+1;
            end
            mm=mm+1;
        elseif(mm==Mm-2)
            while((X(mm0)<=Xy(mm))&&(mm0<Mm0))
                Zx(mm0) =InterpolateOrderN6(5,X(mm0),Xy(mm-3),Xy(mm-2),Xy(mm-1),Xy(mm),Xy(mm+1),0,Z1(mm-3),Z1(mm-2),Z1(mm-1),Z1(mm),Z1(mm+1),0);
                XCdf(mm0)=XCdf(mm0)+normcdf(real(Zx(mm0))).*ZIProb(nn);
                if((mm0>1)&&(Zx(mm0)<Zx(mm0-1)))
                    Zx(mm0)=Zx(mm0-1);
                end
                mm0=mm0+1;
            end
            mm=mm+1;
        elseif(mm==Mm-1)
            while((X(mm0)<=Xy(mm))&&(mm0<Mm0))
                Zx(mm0) =InterpolateOrderN6(4,X(mm0),Xy(mm-3),Xy(mm-2),Xy(mm-1),Xy(mm),0,0,Z1(mm-3),Z1(mm-2),Z1(mm-1),Z1(mm),0,0);
                if((mm0>1)&&(Zx(mm0)<Zx(mm0-1)))
                    Zx(mm0)=Zx(mm0-1);
                end
                XCdf(mm0)=XCdf(mm0)+normcdf(real(Zx(mm0))).*ZIProb(nn);
                mm0=mm0+1;
            end
            while((X(mm0)>Xy(mm))&&(mm0<Mm0))
                Zx(mm0) =InterpolateOrderN6(4,X(mm0),Xy(mm-3),Xy(mm-2),Xy(mm-1),Xy(mm),0,0,Z1(mm-3),Z1(mm-2),Z1(mm-1),Z1(mm),0,0);
                if((mm0>1)&&(Zx(mm0)<Zx(mm0-1)))
                    Zx(mm0)=Zx(mm0-1);
                end
                XCdf(mm0)=XCdf(mm0)+normcdf(real(Zx(mm0))).*ZIProb(nn);
                mm0=mm0+1;
            end
        end   
    end
    
end

end
.
.
.
function [yI,ZI,NnI,dNnI,ZIProb] = InterpolateAssetAlongVolatility(y,Z,wnStart,Nn,dNn,ynStart,Mm)
%UNTITLED4 Summary of this function goes here
%   Detailed explanation goes here

    x01=Z(wnStart)+1/3*dNn;
    x02=Z(wnStart)+2/3*dNn;
    x1=Z(wnStart);
    x2=Z(wnStart+1);
    x3=Z(wnStart+2);
    x4=Z(wnStart+3);
    x5=0;
    x6=0;
 

    y1=y(wnStart,ynStart:Mm);
    y2=y(wnStart+1,ynStart:Mm);
    y3=y(wnStart+2,ynStart:Mm);
    y4=y(wnStart+3,ynStart:Mm);
    y5=0;
    y6=0;
    
    N=4;
    [y01] = InterpolateOrderN6(N,x01,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y02] = InterpolateOrderN6(N,x02,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);

    yI(wnStart,ynStart:Mm)=y(wnStart,ynStart:Mm);
    yI(wnStart+1,ynStart:Mm)=y01(ynStart:Mm);
    yI(wnStart+2,ynStart:Mm)=y02(ynStart:Mm);
    
    
    x01=Z(wnStart+1)+1/3*dNn;
    x02=Z(wnStart+1)+2/3*dNn;
    x1=Z(wnStart);
    x2=Z(wnStart+1);
    x3=Z(wnStart+2);
    x4=Z(wnStart+3);
    x5=Z(wnStart+4);
    x6=0;
    

    y1=y(wnStart,ynStart:Mm);
    y2=y(wnStart+1,ynStart:Mm);
    y3=y(wnStart+2,ynStart:Mm);
    y4=y(wnStart+3,ynStart:Mm);
    y5=y(wnStart+4,ynStart:Mm);
    y6=0;
    N=5;
    [y01] = InterpolateOrderN6(N,x01,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y02] = InterpolateOrderN6(N,x02,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);

    yI(wnStart+3,ynStart:Mm)=y(wnStart+1,ynStart:Mm);
    yI(wnStart+4,ynStart:Mm)=y01(ynStart:Mm);
    yI(wnStart+5,ynStart:Mm)=y02(ynStart:Mm);
    

for nn=wnStart+2:Nn-3
    x01=Z(nn)+1/3*dNn;
    x02=Z(nn)+2/3*dNn;
    x1=Z(nn-2);
    x2=Z(nn-1);
    x3=Z(nn);
    x4=Z(nn+1);
    x5=Z(nn+2);
    x6=Z(nn+3);

    y1=y(nn-2,ynStart:Mm);
    y2=y(nn-1,ynStart:Mm);
    y3=y(nn,ynStart:Mm);
    y4=y(nn+1,ynStart:Mm);
    y5=y(nn+2,ynStart:Mm);
    y6=y(nn+3,ynStart:Mm);
    
    N=6;
    
    [y01] = InterpolateOrderN6(N,x01,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y02] = InterpolateOrderN6(N,x02,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    
    yI(3*nn-2,ynStart:Mm)=y(nn,ynStart:Mm);
    yI(3*nn-1,ynStart:Mm)=y01(ynStart:Mm);
    yI(3*nn,ynStart:Mm)=y02(ynStart:Mm);
    
    
end

    x01=Z(Nn-2)+1/3*dNn;
    x02=Z(Nn-2)+2/3*dNn;
    x1=Z(Nn-4);
    x2=Z(Nn-3);
    x3=Z(Nn-2);
    x4=Z(Nn-1);
    x5=Z(Nn);
    x6=0;
    

    y1=y(Nn-4,ynStart:Mm);
    y2=y(Nn-3,ynStart:Mm);
    y3=y(Nn-2,ynStart:Mm);
    y4=y(Nn-1,ynStart:Mm);
    y5=y(Nn,ynStart:Mm);
    y6=0;
    N=5;
    [y01] = InterpolateOrderN6(N,x01,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y02] = InterpolateOrderN6(N,x02,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);

    yI(3*(Nn-2)-2,ynStart:Mm)=y(Nn-2,ynStart:Mm);
    yI(3*(Nn-2)-1,ynStart:Mm)=y01(ynStart:Mm);
    yI(3*(Nn-2),ynStart:Mm)=y02(ynStart:Mm);
    
    x01=Z(Nn-1)+1/3*dNn;
    x02=Z(Nn-1)+2/3*dNn;
    x1=Z(Nn-3);
    x2=Z(Nn-2);
    x3=Z(Nn-1);
    x4=Z(Nn);
    x5=0;
    x6=0;
    

    y1=y(Nn-3,ynStart:Mm);
    y2=y(Nn-2,ynStart:Mm);
    y3=y(Nn-1,ynStart:Mm);
    y4=y(Nn,ynStart:Mm);
    y5=0;
    y6=0;
    N=4;
    [y01] = InterpolateOrderN6(N,x01,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);
    [y02] = InterpolateOrderN6(N,x02,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6);

    yI(3*(Nn-1)-2,ynStart:Mm)=y(Nn-1,ynStart:Mm);
    yI(3*(Nn-1)-1,ynStart:Mm)=y01(ynStart:Mm);
    yI(3*(Nn-1),ynStart:Mm)=y02(ynStart:Mm);
    
    yI(3*Nn-2,ynStart:Mm)=y(Nn,ynStart:Mm);
    
    
    for nn=1:Nn-1
       ZI(nn*3-2)=Z(nn);
       ZI(nn*3-1)=Z(nn)+1/3*dNn;
       ZI(nn*3)=Z(nn)+2/3*dNn;
    end
    ZI(3*Nn-2)=Z(Nn);
    
    NnI=3*Nn-2;
    dNnI=1/3*dNn;
    
    ZIProb(wnStart)=normcdf((.5*ZI(wnStart)+.5*ZI(wnStart+1)));
    for nn=wnStart+1:NnI-1
        ZIProb(nn)=normcdf((.5*ZI(nn)+.5*ZI(nn+1)))-normcdf((.5*ZI(nn)+.5*ZI(nn-1)));
    end
    
    ZIProb(NnI)=1-normcdf((.5*ZI(NnI)+.5*ZI(NnI-1)));
    
    
    end
 
User avatar
Amin
Topic Author
Posts: 2853
Joined: July 14th, 2002, 3:00 am

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

September 9th, 2021, 5:02 pm

Here is the last output you would have when you run the posted program.

ItoHermiteSVMean =
   0.098525095594031
TrueSVMean =
   0.098683674566407
ItoHermiteAssetMeanX2 =
   0.999916714578683
TrueAssetMean =
     1
SVolMeanAnalytic =
   0.098683674566407
SVolMeanMC =
   0.098605764361773
AssetMeanAnalytic =
     1
AssetMeanMC =
   1.000114617811523
IndexMax =
   751


Image

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

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

September 9th, 2021, 6:05 pm

Please in the main program, replace this 
for nn=wnStart:Nn
    for mm=wnStart:Mm
        if(y(nn,mm)>0)
            dy0(nn,ynStart(nn):Mm)=(yMu0dt(nn,ynStart(nn):Mm)+ Zy2(nn,ynStart(nn):Mm)-Zy1(nn,ynStart(nn):Mm)+yDiff0(nn,ynStart(nn):Mm) ...
                ).*dydZ0(nn,ynStart(nn):Mm).^-1;
        else
            dy0(nn,ynStart(nn):Mm)=real((yMu0dt(nn,ynStart(nn):Mm)+ Zy2(nn,ynStart(nn):Mm)-Zy1(nn,ynStart(nn):Mm))).*real(dydZ0(nn,ynStart(nn):Mm)).^-1; 
        end     
    end
end
with this
,
for nn=wnStart:Nn
    
            dy0(nn,ynStart(nn):Mm)=(yMu0dt(nn,ynStart(nn):Mm)+ Zy2(nn,ynStart(nn):Mm)-Zy1(nn,ynStart(nn):Mm)+yDiff0(nn,ynStart(nn):Mm) ...
                ).*dydZ0(nn,ynStart(nn):Mm).^-1;
           
end
 
User avatar
Amin
Topic Author
Posts: 2853
Joined: July 14th, 2002, 3:00 am

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

September 10th, 2021, 4:56 am

Friends, I updated my resume and I am posting a copy of it here:

Ahsan Amin
Email: anan2999@yahoo.com
Phone: +92-336-2602125
Website: https://ahsanamin2999.wordpress.com/
Research Website: https://papers.ssrn.com/sol3/cf_dev/Abs ... _id=435366
Education
New York University, New York, NY, USA.
Masters of Arts in Economics with a concentration in Mathematical Finance, September 2002.
Wrote my thesis on pricing of Bermudan Swaptions and other Bermudan structured Derivatives within the framework of Libor Market Model. Libor Market Model was the newly introduced and most modern interest rate model at that time. My thesis was titled, “Pricing Bermudan fixed income derivatives in multi-factor extended LIBOR market model.” My thesis was cited in several other research papers/theses.
Northwestern University, Evanston, IL, USA.

Bachelor of Science in Electrical Engineering, March 97

Work
CEO, Infiniti Derivatives Technologies, Lahore, Pakistan
Feb 2012 - Present
1.       My company  web site is: http://www.infinitiderivatives.com
 
2.      Worked on research on the analytic solution of ordinary differential equations. I presented methods for analytic solution of general nth-order ordinary differential equations and also for the analytic solution of systems of nth-order ordinary differential equations. My research is discussed in the research paper, "On the General Solution of Initial Value Problems of Ordinary Differential Equations Using the Method of Iterated Integrals" downloadable at https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2872598
 
3.      I found solution to the problem of high order monte carlo simulations of Stochastic differential equations and now you can easily do accurate monte carlo simulations to arbitrarily higher order based on my work. This is probably my best and most useful work. I have been able to simulate with great precision a correlated stochastic volatility system of two SDEs with arbitrary drift and volatility parameters in each of the SDE with this method.
 
4.      I worked for a very accurate method to find the numerical solution of Fokker-Planck equation for SDEs in original and Bessel coordinates. The method is extremely fast, explicit and works where many other methods fail. The new method is based on principle of conservation of probability mass and results from an equation that I call sister equation to continuity equation of physics.
 
5.      I have extended the method of solution of Fokker-planck equations by probability mass conservation to two dimensions to find out the solution of a general 2D PDE where underlying SDEs have arbitrary parameters. This method nests all models like Heston, SABR, ZABR, Cheyette and many other two dimensional stochastic volatility models but some more work remains to be done to nest all these models.

6.      I suggested a new method to evolve and discretize the solution of SDEs on a transition probabilities grid. My method represents the next generation version of quantizers algorithm. My next research goal is to do simulation of smooth and non-smooth functionals of SDEs on this transition probabilities grid.
 
7.       I suggested a new method to construct the densities that can match moments of any order to the density with precision. The method ensures that density is perfectly positive, completely normalized, and simultaneously matches all given moments. Now you would not need Edgeworth expansions, or Gram-Charlier because my suggested method works better.
 
8.      My research paper titled, “Calibration, Simulation and Hedging in a Heston Libor Market Model with Stochastic Basis” was an invited publication in a Risk book “Interest Rate Modelling after the Financial Crisis.” Here is the link to Risk book web site: http://riskbooks.com/interest-rate-mode ... ial-crisis  
 
9.      I also successfully completed another research project which resulted in the research paper titled “Solution of Stochastic Volatility Models Using Variance Transition Probabilities and Path Integrals.” This paper can be downloaded at https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2149231
 
10.  I worked with a method to predict an arbitrary volatile time series of the kind that appears in stocks and other financial instruments trading using slowly moving parameters method applied through Kalman filtering. This method can also be used to model the degree of mean-reversion in financial instruments and helps in mean-reversion and convergence trades.
 
 
Quantitative Research Analyst, UP-FRONT Inc., Tokyo, Japan
Feb 2003 – Oct 2010


It was a research and development Role and I was based in Lahore, Pakistan. These were the major research projects I completed with Upfront.
 
1.      I developed a stochastic volatility displaced diffusion LIBOR Market Model.  Worked on global Calibration of the model to more than 2000 swaptions with forward curve out to sixty years. Priced callable structured derivatives and more exotic deals as well as their deltas and vegas in this stochastic volatility model.
 
2.      Worked on advanced multi-factor skew extended calibration of LIBOR Market model. Very robust smoothing constraints were implemented which were required for stability of hedges. Worked with calibration of four different versions of skew extended LIBOR Market Models. Successfully dealt with the problem of simultaneous calibration of the whole swaptions matrix.

3.      Worked extensively with pricing of most complex callable LIBOR exotics. The instruments include Bermudans swaptions, Bermudan callable reverse floaters, Bermudan captions, Bermudan CMS floaters, Bermudan CMS reverse floaters, Bermudan CMS spread options, CMS TARN structures, Bermudan Snowballs, snowblades, and Bermudan callable PRDC.
 
4.      Extended Longstaff and Schwartz method to deal with complex callable fixed income structured derivatives in my research paper, “Multi-Factor Cross Currency Libor Market Models: Implementation, Calibration and Examples.” The paper can be downloaded here: https://papers.ssrn.com/sol3/papers.cfm ... id=1214042
 
5.      Worked on the pricing of Bermudan Callable and Trigger Power Reverse Dual Currency Note and other FX/IR hybrid derivatives. Developed a multi-factor cross-currency LIBOR Market Model for the pricing of these instruments. Implemented the advanced analytics of the model on Excel using VBA and C++.
 
6.      Worked on a three factor Hull White Cross Currency model in Cheyette framework for fast pricing of Bermudan trigger PRDCs. This involved analytic calibration of the model to FX volatilities and to swaption volatilities of local and foreign economy.
 
User avatar
Amin
Topic Author
Posts: 2853
Joined: July 14th, 2002, 3:00 am

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

September 10th, 2021, 9:14 am

Friends, I continued to work to improve the behavior of program at zero and made some more improvements. When I was working on the program in my room, mind control agents continued to repeatedly release very very sickening gas in my room and now I have very little energy left to do any more work. The gas continued to sap all energy out of me. During past few weeks, mind control agents have got far bolder in releasing the gas in my room ,out of the ACs in every room in the house. Several rooms in our house had better air-conditioning but after I would work there for a few days, the would start releasing large amounts of gas from the ACs. Now I do not turn on air-conditioning despite very high humidity but they continue to release large amount of gas in my room and other places where I work. 
I continued to write for a week to ten days about mind control torture on me but the torture was not decreasing so I just simply gave up on writing about it everyday.
As a succession of what I revealed a few weeks ago, I will like to tell friends that a jewish billionaire is Godfather of mind control. Though there are far richer billionaires in US who are staunchly opposed to mind control, they are not ready to shell out hundreds of millions of dollars to bribe anyone to end mind control while mind control Godfather  continues to give large bribes to people to continue mind control. Many of the people related to defense in mind control are given extremely high paying jobs in Godfather's companies after they retire from mind control and Godfather is extremely generous towards them and it is well established and well known to people in mind control after retarding muslims, blacks and others in their career, they will move to paradise after they retire. No wonder there is no let up in mind control whatever good people try. So much so for goodness in American society where at the end of the day it is money that controls everything.
We know in almost all classic folk lore stories evil is too strong and too sure of itself while good is mostly weak and just trying to protect itself from evil who is bent on damaging and eradicating the weak and the poor good but it is the weak and poor good that finally prevails. I want to tell good American friends that evil at this stage has become too bold and thousands of talented people are retarded every year. And it is just not muslims and blacks that are on the line, several hundreds of brilliant whites are retarded across United States by mind control agency every year that has become too emboldened due to lack of any accountability. Thousands of foreigners are also kept retarded in dozens of countries including many European countries. 
I want to tell friends something I have said before here that I am sure American society would eventually end the menace of mind control from their and other countries. I want to request all good Americans and Europeans to come forward and invite any investigative journalist you know to read my threads and ask them to help end the menace of mind control in our societies through their investigative journalism. Without any accountability evil will continue to grow stronger and then mind control would start to threaten even those people who consider themselves safe due to affiliations, status or wealth.
I want to request CNN, New York Times, Washington Post and large reputed European media outlets to run a comprehensive story on mind control torture, retarding and victimization of intelligent people of US and other nations by crooks in US army on behest of some powerful people in United States and due to rightwing extremist biases of these crooks. If you would like to do investigative research about animal practices of US army, one great resource would be mainland European embassies in Muslim countries who keep a detailed account of mind control persecution of intelligent muslims in these countries by crooks in US army. Many of the staff in mainland European embassies are very good human beings who abhor such practices and would love to cooperate with good journalists in exposing the evil animal practices of US army crooks. Only the accounts of people at mainland European embassies in muslim countries thorough animal practices of American army's mind control wing would be enough to drop a great bombshell in the media and general public all across the world. I want to warn all the good people who try to have a civilized dialogue with crooks in American army to end their evil practices that their being nice with crooks is a very misguided approach. Since good people do not have the power to forcefully end the evil practices of US army crooks, only way to end the evil practices of US army crooks is by exposing them openly in US public and all across the world. 
It seems that crooks in American army and mind control  are unable to pay a heed to civilized calls by good American people and others to end their animal cruel practices and remain bent on mischief. If crooks in American army and mind control remain adamant on continuing evil practices, the only solution seems that some brave investigative journalists expose them by running a comprehensive story on mind control practices by American army. It seems that doing a civilized dialogue with hardened crooks in mind control gives them a strong sense of weakness of good people and makes the crooks even more adamant to continue their evil practices.  Reminds me of  Abu-Gharib. Very similarly, Crooks in defense had no conception that they were doing any wrong thing when there was a culture of openly urinating on human captives. It was an open thing in army and most crooks thought it was indeed a very right thing to do( and I am sure some of those crooks still believe that it was a right thing they did even after being reprimanded by the broader American society). It was not until brave people at CNN did a daring story against animal practices by crooks that evil practices stopped. Though American army crooks retard intelligent people of all color and creed, these ultra right wing army crooks love to retard blacks and muslims with great relish. Blacks are only 8-10% of US population but make a very large proportion of victims in united states since many ultra right wing crooks in US army would rather die than let those intelligent blacks succeed in American society in a big way.
 
User avatar
Amin
Topic Author
Posts: 2853
Joined: July 14th, 2002, 3:00 am

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

September 10th, 2021, 6:00 pm

Friends, I am posting a transition probabilities calculation program for a 1 dimensional SDE. Friends would recall, at the start of this year, when I was trying to do the same thing, we had to develop a lot of infrastructure for that but now with precise calculations of  1D SDEs with probability conservation method, it has become almost trivial to construct transition probabilities on this grid. I quickly wrote this program in past two hours and the purpose is to do a full scale monte carlo like simulations with path integrals of smooth or non-smooth functionals of SDEs like forward densities of swaps or other structured payoffs and later we could do a backward dynamic programming on the same grid to find the values of American/bermudan/callable instruments and their hedges. I want to do this on 1D grid first and then take it to 2D grid after we have successfully tackled 1D SDEs.
I have calculated transition probabilities in bessel coordinates but please note that they remain the same across bessel transformation and can directly be used for original coordinates.
Here is the intermediate program I quickly wrote with transition probabilities calculations. Please note how simple these calculations are as compared to what we were trying to do at the start of this year.
function [] = FPERevisitedTransProb08ABwNew04DownloadedAndPI()

%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/2/2/2;   % Simulation time interval.%Fodiffusions close to zero
             %decrease dt for accuracy.
Tt=128*2*2*2*2;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year; 
T=Tt*dt;
OrderA=4;  %
OrderM=4;  %
dtM=.0625/2;
TtM=T/dtM;


dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=50;  % No of normal density subdivisions
NnMidl=25;%One half density Subdivision left from mid of normal density(low)
NnMidh=26;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;


x0=1;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.95;%50;   % volatility power.                                                                                                                                                                                                                                                                     
kappa=4.0;%.950;   %mean reversion parameter.
theta=.25;%mean reversion target
sigma0=1.50;%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;                
w(1:Nn)=x0^(1-gamma)/(1-gamma);
ws2(1:Nn)=w(1:Nn);

%Z(1:Nn)=(((1:Nn)-5.5)*dNn-NnMid);
Z(1:Nn)=(((1:Nn)-5.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);
    %Above calculate probability mass in each probability subdivision.

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

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

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


wnStart=1;%
ss=0;
tic

Zt1(wnStart:Nn)=0.0;
Zt2(wnStart:Nn)=0.0;
for tt=1:Tt  
    
    [wMu0dt,dwMu0dtdw,c1] = CalculateDriftAndVolA4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
        
    [wMid] = InterpolateOrderN8(8,0,Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidl),Z(NnMidh),Z(NnMidh+1),Z(NnMidh+2),Z(NnMidh+3),w(NnMidl-3),w(NnMidl-2),w(NnMidl-1),w(NnMidl),w(NnMidh),w(NnMidh+1),w(NnMidh+2),w(NnMidh+3));
        
    Zt1(wnStart:Nn)=w(wnStart:Nn)-wMid;
        
    [dwdZ,d2wdZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt1,Z);
        
    C0(wnStart:Nn)=Zt1(wnStart:Nn)-dwdZ(wnStart:Nn).*Z(wnStart:Nn);
        
    %Zt2(wnStart:Nn)=C0(wnStart:Nn)+abs(sqrt((dwdZ(wnStart:Nn)).^2+sigma0^2*dt)).*Z(wnStart:Nn);
    Zt2(wnStart:Nn)=C0(wnStart:Nn)+abs(sqrt((dwdZ(wnStart:Nn)).^2+c1(wnStart:Nn).^2)).*Z(wnStart:Nn);
    
    [dw2dZ,d2w2dZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt2,Z);
    if(tt>20)
        w(wnStart:Nn)=w(wnStart:Nn)+ ...
            wMu0dt(wnStart:Nn)+Zt2(wnStart:Nn)-Zt1(wnStart:Nn)+ ...
            0*.5*sigma0^2*(dwdZ(wnStart:Nn).^(-2)).*d2wdZ2A(wnStart:Nn).*dt+ ...
            .5*sigma0^2*(dw2dZ(wnStart:Nn).^(-2)).*d2w2dZ2A(wnStart:Nn).*dt;
    else
        w(wnStart:Nn)=w(wnStart:Nn)+ ...
            wMu0dt(wnStart:Nn)+Zt2(wnStart:Nn)-Zt1(wnStart:Nn);%+ ... 
    end
    %yy(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
        
    [wE] = InterpolateOrderN6(6,Z(Nn)+dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));
     
     w1(1:Nn-1)=w(1:Nn-1);
     w1(Nn)=w(Nn);
     w2(1:Nn-1)=w(2:Nn);
     w2(Nn)=wE;          
     w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
      %Change 3:7/25/2020: I have improved zero correction in above.
     w(w<0)=0.0;
     for nn=1:Nn
         if(w(nn)<=0)
            wnStart=nn+1;
         end
     end
     
     
     if(rem(tt,10)==0) %Since original grid is too fine in time, we want to find transition probabilities between grids every ten time steps apart.
        ss=ss+1;
        ds=dt*10;
        ws1(wnStart:Nn)=ws2(wnStart:Nn);
        ws2(wnStart:Nn)=w(wnStart:Nn);
     
        [ws2a,ws2b] = CalculateGridStartsAndEndsBoundedEnds(ws2,Z,wnStart,Nn,dNn);
        
        
        [Muw,Sigmaw,Sigma2w,Sigma3w] = CalculateDriftAndVolA404B4Drift02(ws1,wnStart,Nn,YqCoeff0,Fp1,gamma,ds);
        Pn(1:Nn)=0.0;
        for mm=wnStart:Nn
            Zt(wnStart:Nn)=(2*(-ws2b(wnStart:Nn)+ws1(mm)+Muw(mm)-Sigma2w(mm)))./(-Sigmaw(mm) - sqrt(Sigmaw(mm).^2-4*Sigma2w(mm).*(-ws2b(wnStart:Nn)+ws1(mm)+Muw(mm)-Sigma2w(mm))));
                       
            Zt(wnStart:Nn)=real(Zt(wnStart:Nn));
            Pm0n(wnStart:Nn)=normcdf(Zt(wnStart:Nn),0,1);
            %Pmn below is integrated transition probability between grid
            %cell mm at time (ss-1)*ds to cells wnStart:Nn on next grid at
            %time ss*ds.
            Pmn(mm,wnStart)=Pm0n(wnStart).*ZProb(mm);
            Pmn(mm,wnStart+1:Nn-1)=(Pm0n(wnStart+1:Nn-1)-Pm0n(wnStart:Nn-2)).*ZProb(mm);
            Pmn(mm,Nn)=(1-Pm0n(Nn-1)).*ZProb(mm);
        %    Pn(wnStart:Nn)=Pn(wnStart:Nn)+Pmn(mm,wnStart:Nn);
        end
        
        %plot((wnStart:Nn),Pn(wnStart:Nn),'r',(wnStart:Nn),ZProb(wnStart:Nn),'g');
        %str=input('Look at graph comparison')
     end
     
end
yy(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(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

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

toc

ItoHermiteMean=sum(yy(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %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*dt*Tt)%Mean reverting SDE original variable true average

theta1=1;
rng(29079137, 'twister')
paths=200000;
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+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %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*dt*Tt)%Mean reverting SDE original variable true average


MaxCutOff=30;
NoOfBins=round(300*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
.
.
.
Here is the output of this program.
ItoHermiteMean =
   0.250230651876615
true Mean only applicable to standard SV mean reverting type models otherwise disregard
TrueMean =
   0.250251596970927
Original process average from monte carlo
MCMean =
   0.249542499753951
Original process average from our simulation
ItoHermiteMean =
   0.250230651876615
true Mean only applicble to standard SV mean reverting type models otherwise disregard
TrueMean =
   0.250251596970927
IndexMax =
   651
red line is density of SDE from Ito-Hermite method, green is monte carlo.
 
User avatar
Amin
Topic Author
Posts: 2853
Joined: July 14th, 2002, 3:00 am

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

September 12th, 2021, 8:37 pm

Friends, I was able to make great progress on simulation of densities of smooth and non-smooth functionals on transition probabilities grid. I was able to write very simple matrix equations that perfectly retrieve the expectations of functionals of SDE forward into time on the transition probabilities grid. I had continued to try this earlier without any luck but now it turned out to be  very easy and quite simple. There are still some very minor issues that I hope to resolve in next 2-3 days. I have made the major breakthrough and I will be sharing the analytics through latex equations very soon and of course post the program as well.
 
User avatar
Amin
Topic Author
Posts: 2853
Joined: July 14th, 2002, 3:00 am

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

September 13th, 2021, 7:24 pm

Friends, I have been able to match analytic densities of path integrals of SDEs from transition probability grid with path integral densities from monte carlo. Now it works quite well for mean-reverting SDEs. I will be finessing things further and match analytic densities of path integrals of functions of SDEs with their monte carlo densities tomorrow. And also work on weights on path integrals that are not smooth in time and change arbitrarily across time. I hope to be posting initial programs in two to three days.