SERVING THE QUANTITATIVE FINANCE COMMUNITY

Amin
Topic Author
Posts: 2421
Joined: July 14th, 2002, 3:00 am

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

Dear friends, I was trying to think how to derive the proper equations for the evolution of density in Ito-hermite method. I was able to find some other interesting equations that I decided to share here with friends.
Before I write the main equations, I want to mention that this idea is based on equivalence between two densities. I take one density that is standard normal called Z (could also be taken as brownian motion) and then I try to find the evolution of a particular point on the density of SDE variable X(t,Z) and this particular point is related to a specific fixed value of Z which is already known. Z is a static density that does not change with time at all.
I will write a few equations first that we will use later in the main derivation. Here are the equations

$P(Z)=P(X,t) |{\frac{\partial X}{\partial Z}}|$   Follows from change of variable formula for densities.
$\frac{\partial P(X,t)}{\partial X}= \frac{\partial p(Z)}{\partial Z} (\frac{\partial Z}{\partial X})^2$

$\frac{d}{dt} \big[ {\frac{\partial X}{\partial Z}} \big]=\frac{\partial}{\partial X}\big[ {\frac{\partial X}{\partial Z}} \big] \frac{\partial X}{\partial t}$
$=-{(\frac{\partial X}{\partial Z})}^2 \frac{\partial^2 Z}{\partial X^2} \frac{\partial X}{\partial t}$

$\frac{\partial P(X,t)}{\partial t} = f(X)$ where f(X) is the fokker planck equation.

Now we come towards the main equation derivation(which might as well possibly be wrong)

Since standard normal is a static density and it does not change with time, we can write

$\frac{d}{dt}[p(Z)]=0$
since we can write densities in terms of each other using change of variable formula, we can write
$\frac{d}{dt}[p(X,t)\frac{\partial X}{\partial Z}]=0$
taking the derivatives, we can expand the above equation as
$\big [ \frac{\partial p(X,t)}{\partial t}+\frac{\partial p(X,t)}{\partial X}\frac{\partial X}{\partial t} ]\frac{\partial X}{\partial Z}$
$-p(X,t){(\frac{\partial X}{\partial Z})}^2 \frac{\partial^2 Z}{\partial X^2} \frac{\partial X}{\partial t}=0$

using
$P(Z)=P(X,t) |{\frac{\partial X}{\partial Z}}|$
$\frac{\partial P(X,t)}{\partial X}= \frac{\partial p(Z)}{\partial Z} (\frac{\partial Z}{\partial X})^2$
and $\frac{\partial P(X,t)}{\partial t} = f(X)$
in the above equation ,we get
$\frac{\partial X(Z,t)}{\partial t} =\frac{ \frac{\partial p(X,t)}{\partial t}}{ \big[ (p(Z) \frac{\partial^2 Z}{\partial X^2}+ Z p(Z) ({\frac{\partial Z}{\partial X}})^2\big]}$

using the above equation, we can probably solve every SDE so that we can solve for evolution of X as a function of Z and t. More detailed post and notes tomorrow.
Sorry It was too late yesterday when I wrote the above post and it took me more than one and a half hour to write it and I slept and woke up several times on my desk in between. So I resume yesterday's post.
copying from the above post we have
$\frac{\partial X(Z,t)}{\partial t} =\frac{ \frac{\partial p(X,t)}{\partial t}}{ \big[ (p(Z) \frac{\partial^2 Z}{\partial X^2}+ Z p(Z) ({\frac{\partial Z}{\partial X}})^2\big]}$

where $\frac{\partial p(X,t)}{\partial t}$ is the fokker planck equation for the related SDE.
We can use the change of variable formula in fokker-planck again to get rid of the densities and get an evolution equation solely in terms of Z and SDE variable X(Z,t). Suppose the structure of fokker planck is given as

$\frac{\partial p(X,t)}{\partial t}=a(X,t) + b(X,t) p(X,t) + c(X,t) \frac{\partial p(X,t)}{\partial X} + d(X,t) \frac{\partial^2 p(X,t)}{\partial X^2}$

when we have $a(X,t)=0$ we can eliminate the densities from the equation by converting p(X,t) into p(Z) and cancelling the p(Z) in the equation. using the following formulas, we get

$p(X)=p(Z) |\frac{\partial Z}{\partial X}|$  This follows from the standard change of variables in a density.
$\frac{\partial p(X)}{\partial X}=\frac{\partial p(Z)}{\partial Z} {(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2} =-Z p(Z){(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2}$
$\frac{\partial^2 p(X)}{\partial X^2}=\frac{\partial^2 p(Z)}{\partial Z^2} {(\frac{\partial Z}{\partial X})}^3+3 \frac{\partial p(Z)}{\partial Z} \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2}+p(Z) \frac{\partial^3 Z}{\partial X^3}$
$=(Z^2-1) p(Z) {(\frac{\partial Z}{\partial X})}^3 -3 Z p(Z) \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} +p(Z) \frac{\partial^3 Z}{\partial X^3}$

$\frac{\partial p(X,t)}{\partial t}= b(X,t) p(Z) |\frac{\partial Z}{\partial X}| + c(X,t) \big[-Z p(Z){(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2}\big]$
$+ d(X,t) \big[ (Z^2-1) p(Z) {(\frac{\partial Z}{\partial X})}^3 -3 Z p(Z) \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} +p(Z) \frac{\partial^3 Z}{\partial X^3} \big]$

putting the above value of fokker planck in the first equation and cancelling p(Z) we get
$\frac{\partial X(Z,t)}{\partial t} =\frac{ b(X,t) |\frac{\partial Z}{\partial X}| + c(X,t) \big[-Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2}\big] + d(X,t) \big[ (Z^2-1) {(\frac{\partial Z}{\partial X})}^3 -3 Z \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} + \frac{\partial^3 Z}{\partial X^3} \big]}{ \big[ \frac{\partial^2 Z}{\partial X^2}+ Z ({\frac{\partial Z}{\partial X}})^2\big]}$
Please be careful with signs as I might have missed some absolute in change of probability derivatives.

In order to work out a proper evolution equation, We can take the value of Z and its derivatives $\frac{\partial Z}{\partial X}$ and $\frac{\partial^2 Z}{\partial X^2}$ and $\frac{\partial Z^3}{\partial X^3}$ as constant at start of the time step and expand the integrals of the kind like $\int_0^t b(X,s) H_n(Z) ds$,  $\int_0^t c(X,s)H_n(Z) ds$ and $\int_0^t d(X,s) H_n(Z) ds$ for a particular value of Z related to each grid point using stochastics just like we did for monte carlo simulations and then form a valid evolution equation. In the above integrals $H_n(Z)$ is a standard normal hermite polynomial of nth degree.
These are very basic thoughts that I just worked out quickly and there is also some possibility that we might have to use the equivalence relation with brownian motion density instead of standard normal density but if required doing the same exercise on brownian motion should also be straightforward.
I will be doing my experiments and coming out with a more detailed post in a day.
Again since I recently did this and never numerically tried it there is a chance that we might have to use brownian motion density instead of standard normal density but that should also be straightforward.

Amin
Topic Author
Posts: 2421
Joined: July 14th, 2002, 3:00 am

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

Dear friends, the method proposed in previous post to solve for the densities of SDEs by doing a change of variable with static density of standard normal will work for other SDEs only if it works for the simplest SDE which is that of brownian motion. And doing it for brownian motion SDE will give us insights about how to solve more complicated SDEs since mathematics there can become very tedious and knowing basic principles towards solution of brownian motion SDE will be very helpful. Therefore I tried to solve the brownian motion SDE with this method and it does work. Here are the details.
In what follows, we use the notation $B(t)$ for standard brownian motion and $Z$ for standard normal.

The fokker planck equation for Brownian motion is given as
$\frac{\partial P(B,t)}{\partial t}=.5 \frac{\partial ^2 P(B,t)}{\partial B^2}$  Equation(1)

and we apply time derivative on static standard normal density as
$\frac{d}{dt}[P(Z)]=0$     Equation(2)
Changing to standard brownian motion density inside the brackets by application of change of variable formula
$\frac{d}{dt}[P(B,t)|\frac{\partial B(t)}{\partial Z}|]=0$    Equation(3)
applying the time derivative to equation (3), we get the equation
$[\frac{\partial P(B,t)}{\partial t} +\frac{\partial P(B,t)}{\partial B}\frac{\partial B(t)}{\partial t}]\frac{\partial B(t)}{\partial Z}$
$- P(B,t) (\frac{\partial B(t)}{\partial Z})^2 \frac{\partial^2 Z}{\partial B^2} \frac{\partial B(t)}{\partial t} =0$   Equation(4)
I have removed absolute signs since the derivative between brownian motion and standard normal will be positive.
We put the expression for standard brownian fokker-planck equation from equation(1) to above equation (4)
$[.5 \frac{\partial^2 P(B,t)}{\partial B^2} +\frac{\partial P(B,t)}{\partial B}\frac{\partial B(t)}{\partial t}]\frac{\partial B(t)}{\partial Z}$
$- P(B,t) (\frac{\partial B(t)}{\partial Z})^2 \frac{\partial^2 Z}{\partial B^2} \frac{\partial B(t)}{\partial t} =0$  Equation(5)
now we want to convert the derivatives $\frac{\partial^2 P(B,t)}{\partial B^2}$, $\frac{\partial P(B,t)}{\partial B}$ and $P(B,t)$ again to standard normal density so that we can cancel the densities from our equation.
We use the following equations
$P(B,t)=P(Z)\frac{\partial Z}{\partial B}$
$\frac{\partial P(B,t)}{\partial B}=\frac{\partial P(Z)}{\partial Z}{(\frac{\partial Z}{\partial B})}^2+p(Z)\frac{\partial^2 Z}{\partial B^2}$
$=P(Z)[-Z{(\frac{\partial Z}{\partial B})}^2+\frac{\partial^2 Z}{\partial B^2}]$
$\frac{\partial^2 P(B,t)}{\partial B^2}=\frac{\partial^2 P(Z)}{\partial Z^2}{(\frac{\partial Z}{\partial B})}^3+3 \frac{\partial P(Z)}{\partial Z}\frac{\partial Z}{\partial B} \frac{\partial^2 Z}{\partial B^2}+p(Z)\frac{\partial^3 Z}{\partial B^3}$
$=p(Z)[(Z^2-1) {(\frac{\partial Z}{\partial B})}^3-3 Z \frac{\partial Z}{\partial B} \frac{\partial^2 Z}{\partial B^2}+\frac{\partial^3 Z}{\partial B^3}]$
where I have used
$\frac{\partial^2 P(Z)}{\partial Z^2}=(Z^2-1)p(Z)$ and
$\frac{\partial P(Z)}{\partial Z}=-Z P(Z)$ which are derivatives of standard normal density.
After applying the above formulas in previous equation(5) and cancelling P(Z) that is common throughout, we get
$[.5 [(Z^2-1) {(\frac{\partial Z}{\partial B})}^3-3 Z \frac{\partial Z}{\partial B} \frac{\partial^2 Z}{\partial B^2}+\frac{\partial^3 Z}{\partial B^3}]\frac{\partial B(t)}{\partial Z}$
$+[-Z{(\frac{\partial Z}{\partial B})}^2+\frac{\partial^2 Z}{\partial B^2}]\frac{\partial B(t)}{\partial t}\frac{\partial B(t)}{\partial Z}$
$- (\frac{\partial B(t)}{\partial Z})^2 \frac{\partial^2 Z}{\partial B^2} \frac{\partial B(t)}{\partial t} =0$  Equation(6)
For general solution of SDEs, we will have to find out the derivatives like $\frac{\partial B(t)}{\partial Z}$ and $\frac{\partial^2 B(t)}{\partial Z^2}$ and the third derivative using finite differences or some other analytic procedure. But in case of brownian motion, we know these derivatives and apply them in the above equation
$\frac{\partial B(t)}{\partial Z}=\sqrt{t}$
$\frac{\partial^2 B(t)}{\partial Z^2}=0$
$\frac{\partial^3 B(t)}{\partial Z^3}=0$
and now our above equation (6)  becomes
$.5 (Z^2-1) \frac{\partial Z}{\partial B}-Z \frac{\partial B(t)}{\partial t}=0$ Equation(7)
But $\frac{\partial Z}{\partial B}=\frac{1}{\sqrt{t}}$   and $Z \sqrt{t}=B(t)$
So we get from Equation(7)
$.5 (Z^2-1) -Z \sqrt{t} \frac{\partial B(t)}{\partial t}=0$
$=.5 (Z^2-1) -B(t) \frac{\partial B(t)}{\partial t}=0$
taking dt to the other side of the equation(i.e multiply by dt.), we get
$B(t)dB(t)= .5 (Z^2-1) dt$
or
$B(t)dB(t)+.5 dt = .5 Z^2 dt$
The expression on the left can be solved(by integration on both sides) to get
${B(t)}^2=Z^2 t$ or in general
${B(t_2)}^2={B(t_1)}^2+Z^2 \int_{t_1}^{t_2} ds$   Equation(8)

Of course we will have to take care of signs on Z when taking root on the above formula.
So we can solve for Brownian motion with this method and the good thing is that we learnt that we will have to do some mathematical manipulation of the equations to solve the general SDE and turn the solution expression into squared form and this seems to be very helpful that we know in advance that complicated SDEs will also have to be solved for squared expressions.

Amin
Topic Author
Posts: 2421
Joined: July 14th, 2002, 3:00 am

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

Dear friends, I was trying to think how to derive the proper equations for the evolution of density in Ito-hermite method. I was able to find some other interesting equations that I decided to share here with friends.
Before I write the main equations, I want to mention that this idea is based on equivalence between two densities. I take one density that is standard normal called Z (could also be taken as brownian motion) and then I try to find the evolution of a particular point on the density of SDE variable X(t,Z) and this particular point is related to a specific fixed value of Z which is already known. Z is a static density that does not change with time at all.
I will write a few equations first that we will use later in the main derivation. Here are the equations

$P(Z)=P(X,t) |{\frac{\partial X}{\partial Z}}|$   Follows from change of variable formula for densities.
$\frac{\partial P(X,t)}{\partial X}= \frac{\partial p(Z)}{\partial Z} (\frac{\partial Z}{\partial X})^2$

$\frac{d}{dt} \big[ {\frac{\partial X}{\partial Z}} \big]=\frac{\partial}{\partial X}\big[ {\frac{\partial X}{\partial Z}} \big] \frac{\partial X}{\partial t}$
$=-{(\frac{\partial X}{\partial Z})}^2 \frac{\partial^2 Z}{\partial X^2} \frac{\partial X}{\partial t}$

$\frac{\partial P(X,t)}{\partial t} = f(X)$ where f(X) is the fokker planck equation.

Now we come towards the main equation derivation(which might as well possibly be wrong)

Since standard normal is a static density and it does not change with time, we can write

$\frac{d}{dt}[p(Z)]=0$
since we can write densities in terms of each other using change of variable formula, we can write
$\frac{d}{dt}[p(X,t)\frac{\partial X}{\partial Z}]=0$
taking the derivatives, we can expand the above equation as
$\big [ \frac{\partial p(X,t)}{\partial t}+\frac{\partial p(X,t)}{\partial X}\frac{\partial X}{\partial t} ]\frac{\partial X}{\partial Z}$
$-p(X,t){(\frac{\partial X}{\partial Z})}^2 \frac{\partial^2 Z}{\partial X^2} \frac{\partial X}{\partial t}=0$

using
$P(Z)=P(X,t) |{\frac{\partial X}{\partial Z}}|$
$\frac{\partial P(X,t)}{\partial X}= \frac{\partial p(Z)}{\partial Z} (\frac{\partial Z}{\partial X})^2$
and $\frac{\partial P(X,t)}{\partial t} = f(X)$
in the above equation ,we get
$\frac{\partial X(Z,t)}{\partial t} =\frac{ \frac{\partial p(X,t)}{\partial t}}{ \big[ (p(Z) \frac{\partial^2 Z}{\partial X^2}+ Z p(Z) ({\frac{\partial Z}{\partial X}})^2\big]}$

using the above equation, we can probably solve every SDE so that we can solve for evolution of X as a function of Z and t. More detailed post and notes tomorrow.
Sorry It was too late yesterday when I wrote the above post and it took me more than one and a half hour to write it and I slept and woke up several times on my desk in between. So I resume yesterday's post.
copying from the above post we have
$\frac{\partial X(Z,t)}{\partial t} =\frac{ \frac{\partial p(X,t)}{\partial t}}{ \big[ (p(Z) \frac{\partial^2 Z}{\partial X^2}+ Z p(Z) ({\frac{\partial Z}{\partial X}})^2\big]}$

where $\frac{\partial p(X,t)}{\partial t}$ is the fokker planck equation for the related SDE.
We can use the change of variable formula in fokker-planck again to get rid of the densities and get an evolution equation solely in terms of Z and SDE variable X(Z,t). Suppose the structure of fokker planck is given as

$\frac{\partial p(X,t)}{\partial t}=a(X,t) + b(X,t) p(X,t) + c(X,t) \frac{\partial p(X,t)}{\partial X} + d(X,t) \frac{\partial^2 p(X,t)}{\partial X^2}$

when we have $a(X,t)=0$ we can eliminate the densities from the equation by converting p(X,t) into p(Z) and cancelling the p(Z) in the equation. using the following formulas, we get

$p(X)=p(Z) |\frac{\partial Z}{\partial X}|$  This follows from the standard change of variables in a density.
$\frac{\partial p(X)}{\partial X}=\frac{\partial p(Z)}{\partial Z} {(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2} =-Z p(Z){(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2}$
$\frac{\partial^2 p(X)}{\partial X^2}=\frac{\partial^2 p(Z)}{\partial Z^2} {(\frac{\partial Z}{\partial X})}^3+3 \frac{\partial p(Z)}{\partial Z} \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2}+p(Z) \frac{\partial^3 Z}{\partial X^3}$
$=(Z^2-1) p(Z) {(\frac{\partial Z}{\partial X})}^3 -3 Z p(Z) \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} +p(Z) \frac{\partial^3 Z}{\partial X^3}$

$\frac{\partial p(X,t)}{\partial t}= b(X,t) p(Z) |\frac{\partial Z}{\partial X}| + c(X,t) \big[-Z p(Z){(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2}\big]$
$+ d(X,t) \big[ (Z^2-1) p(Z) {(\frac{\partial Z}{\partial X})}^3 -3 Z p(Z) \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} +p(Z) \frac{\partial^3 Z}{\partial X^3} \big]$

putting the above value of fokker planck in the first equation and cancelling p(Z) we get
$\frac{\partial X(Z,t)}{\partial t} =\frac{ b(X,t) |\frac{\partial Z}{\partial X}| + c(X,t) \big[-Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2}\big] + d(X,t) \big[ (Z^2-1) {(\frac{\partial Z}{\partial X})}^3 -3 Z \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} + \frac{\partial^3 Z}{\partial X^3} \big]}{ \big[ \frac{\partial^2 Z}{\partial X^2}+ Z ({\frac{\partial Z}{\partial X}})^2\big]}$
Please be careful with signs as I might have missed some absolute in change of probability derivatives.

In order to work out a proper evolution equation, We can take the value of Z and its derivatives $\frac{\partial Z}{\partial X}$ and $\frac{\partial^2 Z}{\partial X^2}$ and $\frac{\partial Z^3}{\partial X^3}$ as constant at start of the time step and expand the integrals of the kind like $\int_0^t b(X,s) H_n(Z) ds$,  $\int_0^t c(X,s)H_n(Z) ds$ and $\int_0^t d(X,s) H_n(Z) ds$ for a particular value of Z related to each grid point using stochastics just like we did for monte carlo simulations and then form a valid evolution equation. In the above integrals $H_n(Z)$ is a standard normal hermite polynomial of nth degree.
These are very basic thoughts that I just worked out quickly and there is also some possibility that we might have to use the equivalence relation with brownian motion density instead of standard normal density but if required doing the same exercise on brownian motion should also be straightforward.
I will be doing my experiments and coming out with a more detailed post in a day.
Again since I recently did this and never numerically tried it there is a chance that we might have to use brownian motion density instead of standard normal density but that should also be straightforward.
Though the solution of brownian motion density from fokker-planck turned out to be very simple, the solution of general SDEs should not be very difficult.
As I gave the general expression for the "pre-solution" form which is
$\frac{d X(Z,t)}{dt} =\frac{ b(X,t) |\frac{\partial Z}{\partial X}| + c(X,t) \big[-Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2}\big] + d(X,t) \big[ (Z^2-1) {(\frac{\partial Z}{\partial X})}^3 -3 Z \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} + \frac{\partial^3 Z}{\partial X^3} \big]}{ \big[ \frac{\partial^2 Z}{\partial X^2}+ Z ({\frac{\partial Z}{\partial X}})^2\big]}$

I suppose we have calculated all required derivatives numerically or analytically. Then the above equation can be written as

$\frac{dX(Z,t)}{dt}{ \big[ \frac{\partial^2 Z}{\partial X^2}+ Z ({\frac{\partial Z}{\partial X}})^2\big]}$
$= b(X,t) |\frac{\partial Z}{\partial X}| + c(X,t) \big[-Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2}\big]$
$+ d(X,t) \big[ (Z^2-1) {(\frac{\partial Z}{\partial X})}^3 -3 Z \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} + \frac{\partial^3 Z}{\partial X^3} \big]$

Please note that in our Ito-hermite expanding contracting grid framework, each point on the grid takes a constant value of Z so in X(t,Z) the value of X changes as a function of time and Z remains constant throughout the time evolution of the grid.
We want to change the right hand side of the equation so that after multiplication by dt and taking it to right hand side of previous equation, we have on left hand side of the equation
$dX(Z,t)\big[ \frac{\partial^2 Z}{\partial X^2}+ Z ({\frac{\partial Z}{\partial X}})^2 \big]$
$= d\big [C_0 {X(Z,t)}^2 + C_1 X(Z,t) \big]$

please note that $C_0$ and $C_1$ though different for every point on the X-grid and therefore depends upon X(Z) in a cross-sectional way, it has no role to play in time evolution dynamics as it is not an explicit function of X. Both are like constants that are different for each point on the grid.
We expand
$= d\big [C_0 {X(Z,t)}^2 + C_1 X(Z,t) \big]$
$=C_0 [2 X dX + <dX^2>] + C_1 dX$ here zero in subscript $g_0$ means that it has no role to play in time evolution and depends upon initial cross-sectional value of X.
Now equating first terms of  relevant equations
$dX(Z,t) Z ({\frac{\partial Z}{\partial X}})^2=C_0 2 X dX(Z,t)$
meaning
$C_0=\frac{Z ({\frac{\partial Z}{\partial X}})^2}{2 X(Z) }$
and then we equate
$[C_0 <dX^2> + C_1 dX(Z,t)]=dX(Z,t) \frac{\partial^2 Z}{\partial X^2}$
which would after a little algebra give us the value of $C_1$
so after integration our solution of the equation becomes
$\int d\big [C_0 {X(Z,t)}^2 + C_1 X(Z,t) \big]$
$=\int_0^t \Big[ b(X,t) |\frac{\partial Z}{\partial X}| + c(X,t) \big[-Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2}\big]$
$+ d(X,t) \big[ (Z^2-1) {(\frac{\partial Z}{\partial X})}^3 -3 Z \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} + \frac{\partial^3 Z}{\partial X^3} \big] \Big ] ds$

which can probably be solved after a bit of effort as left side is a complete integral though we will have to do further manipulations on the line of solution of ODEs to get a valid numerical answer.
I will come up with more numerical experiments in a day or two.

Amin
Topic Author
Posts: 2421
Joined: July 14th, 2002, 3:00 am

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

I will be writing a formal paper on the solution of densities of SDEs using my proposed way to solve the fokker-planck equation. Since so many work has already been done, work on the paper should not take more than 10-14 days and I will be updating this thread with numerical material almost everyday. I hope to publish the paper in Wilmott magazine since I have been writing here for so long and Paul was so kind to let me continue this thread despite so many errors I made in the past.
For the sake of reference of friends, I tried to solve the fokker-planck equation several months ago in October 2019 and almost nailed it but I when I tried to test the equations for brownian motion, I failed to complete the total differential of square of brownian motion and therefore thought that there might be some error and it is the "complete differential of square of brownian motion part" that I completed now. You do not have to take the time derivative of normal density as I did in past few posts written very recently. You can try the change of density variable derivatives(with respect to standard normal density) directly in fokker-planck to arrive at exactly the same equations and that is what I had done in early October 2019.   In October last year I was being given antipsychotic injections regularly and therefore I was not able to think very clearly and I really want to thank those friends who pressed on mind control agencies to be better and I was able to get off antipsychotic injections and do better work and became able to suggest the solution to FP equation by making a complete differential as in the simple case of brownian motion and later for general equations.
Here are two posts that I wrote and I copy them here for friends where I had done all the work to solve FP equation in October last year(2019)
Here is my post written on Sun Oct 06, 2019 1:53 pm this is post 812 on this thread. Here is the post link: https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=810#p849953
Amin:
I have made some quick notes on solution of fokker planck equation with hermite polynomials. I have not numerically verified anything so please pardon any errors but the main ideas would stay the same.

Let us suppose we want to solve a general SDE given as
$dX(t)=\mu(X) dt + \sigma(X) dz(t)$
The fokker planck equation of this SDE is given as
$\frac{\partial p(X,t)}{\partial t} = -\frac{\partial [\mu(X) p(X,t)]}{\partial X} + .5 \frac{\partial^2 [{\sigma(X)}^2 p(X,t)]}{\partial X^2}$

First a very brief primer on derivatives of normal and derivatives of the SDE variable.
I denote the standard normal density as p(Z)
We want to convert derivatives of the SDE density of X(t) into derivatives of standard normal density. Here are the main equations
$p(Z)=\frac{1}{\sqrt{2 \pi}} exp[-.5 Z^2]$
derivatives of standard normal density are given in terms of hermite polynomials.
$\frac{\partial p(Z)}{\partial Z}=-Z p(Z)$
$\frac{\partial^2 p(Z)}{\partial Z^2}=(Z^2-1) p(Z)$
now we want to convert the density of X(t) and its derivatives in terms of density of Z. In the density of X, each value of X(t) is associated with a particular value of Z through constant CDF points on corresponding densities. So we calculate the density and derivatives of X(t) in terms of Z.
$p(X)=p(Z) |\frac{\partial Z}{\partial X}|$  This follows from the standard change of variables in a density.
$\frac{\partial p(X)}{\partial X}=\frac{\partial p(Z)}{\partial Z} {(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2} =-Z p(Z){(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2}$
$\frac{\partial^2 p(X)}{\partial X^2}=\frac{\partial^2 p(Z)}{\partial Z^2} {(\frac{\partial Z}{\partial X})}^3+3 \frac{\partial p(Z)}{\partial Z} \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2}+p(Z) \frac{\partial^3 Z}{\partial X^3}$
$=(Z^2-1) p(Z) {(\frac{\partial Z}{\partial X})}^3 -3 Z p(Z) \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} +p(Z) \frac{\partial^3 Z}{\partial X^3}$

Now I write the fokker planck equation after applying the derivatives as
$\frac{\partial p(X,t)}{\partial t} = -\mu(X) \frac{\partial p(X,t)}{\partial X} - \frac{\partial \mu(X)}{\partial X} p(X,t)$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) p(X,t)$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} \frac{\partial p(X,t)}{\partial X}+ .5 {\sigma(X)}^2 \frac{\partial^2 p(X,t)}{\partial X^2}$

Now we write the density of X(t) in terms of density of standard normal Z by substituting from above equations

$\frac{\partial [p(Z) \frac{\partial Z}{\partial X} ]}{\partial t} = -\mu(X) (-Z p(Z){(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2}) - \frac{\partial \mu(X)}{\partial X} p(Z) |\frac{\partial Z}{\partial X}|$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) p(Z) |\frac{\partial Z}{\partial X}|$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} (-Z p(Z){(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2})$
$+ .5 {\sigma(X)}^2 ((Z^2-1) p(Z) {(\frac{\partial Z}{\partial X})}^3 -3 Z p(Z) \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} +p(Z) \frac{\partial^3 Z}{\partial X^3})$

We see that $p(Z)$ which is a static density is common throughout the equations. We eliminate it and write the new equation as

$\frac{\partial [ \frac{\partial Z}{\partial X} ]}{\partial t} = \mu(X) (-Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2}) - \frac{\partial \mu(X)}{\partial X} |\frac{\partial Z}{\partial X}|$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) |\frac{\partial Z}{\partial X}|$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} (-Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2})$
$+ .5 {\sigma(X)}^2 ((Z^2-1) {(\frac{\partial Z}{\partial X})}^3 -3 Z \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} + \frac{\partial^3 Z}{\partial X^3})$

So we get an explicit equation in terms of hermite polynomials describing how  $\frac{\partial [ \frac{\partial Z}{\partial X} ]}{\partial t}$ is locally expanding/contracting in time. We can easily use this to write a numerical program for the evolution of constant CDF lines associated with a particular value of Z ranging from -5 to +5. There can be an infinite amount of analytics that can be done on this form of equation. So we have effectively eliminated the density and we are now directly dealing with how the density is changing/scaling as a function of Z.
Please bookmark this thread as I will be coming with more numerical examples and analytic techniques.
I will also be distributing a numerical program based on these analytics. Hope friends enjoyed reading this.
And here is where I tried to find the solution for brownian motion from the above equations but failed to complete the last step to form a complete differential of the square of brownian motion and therefore did not know what to make of the solution. This post No. 813 was written on  Tue Oct 08, 2019 3:18 pm. Here is the address to the post: https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=810#p850017
Since I was not being able to make acomplete differential of square of brownian motion, I made several silly mistakes (please be kind and pardon the mistakes) to somehow equate the equations and then dropped the idea since it did not seem to work(as I failed to make a complete differential on the left)

I copy the second last equation from the previous post.

$\frac{\partial [p(Z) \frac{\partial Z}{\partial X} ]}{\partial t} = -\mu(X) (-Z p(Z){(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2}) - \frac{\partial \mu(X)}{\partial X} p(Z) |\frac{\partial Z}{\partial X}|$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) p(Z) |\frac{\partial Z}{\partial X}|$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} (-Z p(Z){(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2})$
$+ .5 {\sigma(X)}^2 ((Z^2-1) p(Z) {(\frac{\partial Z}{\partial X})}^3 -3 Z p(Z) \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} +p(Z) \frac{\partial^3 Z}{\partial X^3})$

I had not expnded the left hand side derivative with respect to t. This expansion is
$\frac{\partial p(X)}{\partial t}=\frac{\partial [p(Z) \frac{\partial Z}{\partial X}]}{\partial t}= \frac{\partial p(Z)}{\partial Z} {(\frac{\partial Z}{\partial X})}^2 \frac{\partial X}{\partial t}+ p(Z) \frac{\partial^2 Z}{\partial X^2} \frac{\partial X}{\partial t}$
$=-Zp(Z) {(\frac{\partial Z}{\partial X})}^2 \frac{\partial X}{\partial t}+ p(Z) \frac{\partial^2 Z}{\partial X^2} \frac{\partial X}{\partial t}$

putting this expansion on left side of the first equation, we get
$-Zp(Z) {(\frac{\partial Z}{\partial X})}^2 \frac{\partial X}{\partial t}+ p(Z) \frac{\partial^2 Z}{\partial X^2} \frac{\partial X}{\partial t}$
$= -\mu(X) (-Z p(Z){(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2}) - \frac{\partial \mu(X)}{\partial X} p(Z) |\frac{\partial Z}{\partial X}|$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) p(Z) |\frac{\partial Z}{\partial X}|$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} (-Z p(Z){(\frac{\partial Z}{\partial X})}^2+p(Z) \frac{\partial^2 Z}{\partial X^2})$
$+ .5 {\sigma(X)}^2 ((Z^2-1) p(Z) {(\frac{\partial Z}{\partial X})}^3 -3 Z p(Z) \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} +p(Z) \frac{\partial^3 Z}{\partial X^3})$

eliminating $p(Z)$ from the equation we get an explicit formula for evolution of $X(t)$ as

$\frac{\partial X}{\partial t} (-Z {(\frac{\partial Z}{\partial X})}^2 + \frac{\partial^2 Z}{\partial X^2})$
$= -\mu(X) (-Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2}) - \frac{\partial \mu(X)}{\partial X} |\frac{\partial Z}{\partial X}|$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) |\frac{\partial Z}{\partial X}|$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} (-Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2})$
$+ .5 {\sigma(X)}^2 ((Z^2-1) {(\frac{\partial Z}{\partial X})}^3 -3 Z \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} + \frac{\partial^3 Z}{\partial X^3})$

Let us apply it on toy problem of brownian motion where these partial derivatives with respect to standard normal are known. Please pardon me for some formulas as we are not in a formal brownian motion framework. We are in brownian motion density simulation framework where density is evolving in time along constant CDF lines. Here are some formulas for these constant CDF lines
$B(t)=\sigma \sqrt{t} Z$
$\frac{\partial B}{\partial Z}=\sigma \sqrt{t}$
$\frac{\partial B}{\partial t}=\frac{\sigma}{2 \sqrt{t}} Z$
$\frac{\partial^2 B}{\partial Z^2}=0$
$\frac{\partial^3 B}{\partial Z^3}=0$

Fokker planck for brownian motion is given as
$\frac{\partial p(X)}{\partial t}=.5 \sigma^2 \frac{\partial^2 p(X)}{\partial X^2}$
expanding both sides in our framework into derivatives with respect to density of standard normal and after eliminating the density on both sides we get
$\frac{\partial X}{\partial t} (-Z {(\frac{\partial Z}{\partial X})}^2 + \frac{\partial^2 Z}{\partial X^2})$
$= .5 {\sigma}^2 ((Z^2-1) {(\frac{\partial Z}{\partial X})}^3 -3 Z \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} + \frac{\partial^3 Z}{\partial X^3})$

substituting the known values of partial derivatives, we get

$\frac{\partial B}{\partial t} (-Z {( \frac{1}{\sigma \sqrt{t}} )}^2 )=.5 {\sigma}^2 ((Z^2-1) {( \frac{1}{\sigma \sqrt{t}} )}^3)$
from which follows
$\frac{\partial B}{\partial t} = \frac{(Z^2-1) }{-Z} .5 {\sigma}^2 {( {\sigma \sqrt{t}} )}^2{( \frac{1}{\sigma \sqrt{t}} )}^3)=\frac{(Z^2-1) }{-Z} .5 \frac{\sigma}{\sqrt{t}}$
I believe $\frac{(Z^2-1) }{-Z}$ approaches Z from the properties of hermite polynomials and normal density(I am looking on it if it is an error). Then
$\frac{\partial B}{\partial t} = .5 \frac{\sigma}{\sqrt{t}} Z$ and integrating it we retrieve the original equation of brownian motion along constant CDF lines as
$B(t) = \sigma \sqrt{t} Z$

Amin
Topic Author
Posts: 2421
Joined: July 14th, 2002, 3:00 am

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

Here I re-did the calculations to solve the fokker-planck in Lamperti coordinates. I will post a numerical program in another day possibly tomorrow.
We have the equation from previous posts that we workd out and I copy it here in terms of SDE of w in lamperti coordinates

$\frac{dw}{dt}[Z {(\frac{\partial Z}{\partial w})}^2 + \frac{\partial Z^2}{\partial w^2}]$
$=-\mu(w)[-Z {(\frac{\partial Z}{\partial w})}^2 + \frac{\partial Z^2}{\partial w^2}] -\frac{\partial \mu(w)}{\partial w}{(\frac{\partial Z}{\partial w})}$
$+.5 {\sigma}^2 \big[ (Z^2-1) {(\frac{\partial Z}{\partial w})}^3-3 Z {(\frac{\partial Z}{\partial w})} {(\frac{\partial^2 Z}{\partial w^2})}+{(\frac{\partial Z^3}{\partial w^3})}$

We multiply both sides of the equation by ${(\frac{\partial w}{\partial Z})}^3$ and taking dt to the RHS of the equation, we get

${dw}[Z {(\frac{\partial w}{\partial Z})} +{(\frac{\partial w}{\partial Z})}^3 \frac{\partial Z^2}{\partial w^2}]$
$=\Big[-\mu(w)[-Z {(\frac{\partial w}{\partial Z})} + {(\frac{\partial w}{\partial Z})}^3 \frac{\partial Z^2}{\partial w^2}] -\frac{\partial \mu(w)}{\partial w}{(\frac{\partial w}{\partial Z})}^2$
$+.5 {\sigma}^2 \big[ (Z^2-1)-3 Z{(\frac{\partial w}{\partial Z})}^2 {(\frac{\partial^2 Z}{\partial w^2})}+{(\frac{\partial w}{\partial Z})}^3 {(\frac{\partial Z^3}{\partial w^3})}] \Big] dt$

Manipulating to make the left hand side a complete integral, we first write
$w(Z)=[Z {(\frac{\partial w}{\partial Z})} + {(\frac{\partial w}{\partial Z})}^3 \frac{\partial Z^2}{\partial w^2}]-C_0(Z)$
which means
$[Z {(\frac{\partial w}{\partial Z})} + {(\frac{\partial w}{\partial Z})}^3 \frac{\partial Z^2}{\partial w^2}]=w(Z) + C_0(Z)$
where w(Z) means that w is different on different grid points as a function of Z and for every individual grid point the value of Z remains constant.

putting the above expression in the LHS of fokker-planck equation, we get
$dw[w(Z) + C_0(Z)]$
$=\Big[-\mu(w)[-Z {(\frac{\partial w}{\partial Z})} + {(\frac{\partial w}{\partial Z})}^3 \frac{\partial Z^2}{\partial w^2}] -\frac{\partial \mu(w)}{\partial w}{(\frac{\partial w}{\partial Z})}^2$
$+.5 {\sigma}^2 \big[ (Z^2-1)-3 Z{(\frac{\partial w}{\partial Z})}^2 {(\frac{\partial^2 Z}{\partial w^2})}+{(\frac{\partial w}{\partial Z})}^3 {(\frac{\partial Z^3}{\partial w^3})}\big] \Big]dt$

Now the left hand side is still not a complete integral since we still need to add quadratic variations in Lamperti coordinates $.5 {\sigma}^2 dt$  for square of w. We notice that in the RHS, the $.5 {\sigma}^2 dt$ part in term $.5 {\sigma}^2 (Z^2-1) dt$ has the precisely right quadratic variations for $.5 w^2$ in lamperti coordinates and shift the $.5 {\sigma}^2 dt$ part to the left hand side of the equation and leave  $.5 {\sigma}^2 Z^2 dt$ on the RHS as such and
doing that we get
$w(Z)dw+.5{\sigma}^2 dt + C_0(Z)dw$
$=\Big[-\mu(w)[-Z {(\frac{\partial w}{\partial Z})} + {(\frac{\partial w}{\partial Z})}^3 \frac{\partial Z^2}{\partial w^2}] -\frac{\partial \mu(w)}{\partial w}{(\frac{\partial w}{\partial Z})}^2$
$+.5 {\sigma}^2 \big[ (Z^2)-3 Z{(\frac{\partial w}{\partial Z})}^2 {(\frac{\partial^2 Z}{\partial w^2})}+{(\frac{\partial w}{\partial Z})}^3 {(\frac{\partial Z^3}{\partial w^3})}\big] \Big]dt$

Now the left hand side is a complete integral as
$d[.5 w^2 + C_0 w] =\Big[-\mu(w)[-Z {(\frac{\partial w}{\partial Z})} + {(\frac{\partial w}{\partial Z})}^3 \frac{\partial Z^2}{\partial w^2}] -\frac{\partial \mu(w)}{\partial w}{(\frac{\partial w}{\partial Z})}^2$
$+.5 {\sigma}^2 \big[ (Z^2)-3 Z{(\frac{\partial w}{\partial Z})}^2 {(\frac{\partial^2 Z}{\partial w^2})}+{(\frac{\partial w}{\partial Z})}^3 {(\frac{\partial Z^3}{\partial w^3})}\big] \Big]dt$

And then we get the solution as
$.5 {w(t_1)}^2 + C_0 w(t_1) =.5 {w(t_0)}^2 + C_0 w(t_0)$
$+ \int_{t_0}^{t_1} \Big[ -\mu(w) [-Z {(\frac{\partial w}{\partial Z})} + {(\frac{\partial w}{\partial Z})}^3\frac{\partial^2 Z}{\partial w^2}]$
$-\frac{\partial \mu(w)}{\partial w}{(\frac{\partial w}{\partial Z})}^2$
$+.5 {\sigma}^2 \big[ (Z^2)-3 Z{(\frac{\partial w}{\partial Z})}^2 {(\frac{\partial^2 Z}{\partial w^2})}+{(\frac{\partial w}{\partial Z})}^3 {(\frac{\partial Z^3}{\partial w^3})}\big] \Big]ds$

Amin
Topic Author
Posts: 2421
Joined: July 14th, 2002, 3:00 am

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

I am writing about progress in my work. As I showed in previous posts copied from October 2019, I have developed this approach to solve the fokker-planck equation. I do not keep any research secret and distribute it with friends regularly. Many of the people working with Statistics and probability theory have brilliant minds both in academia and industry. But I will request friends to let me take this research on Fokker-planck towards conclusion and a formal paper. I would certainly love other friends doing research but I hope that they quote my initial effort and resulting formal paper which might be out in a few weeks.
Coming back to research work, the significant progress in yesterday's work was that we came to know that second hermite polynomial contains exactly the right quadratic variations for square of the diffusion but the form of solution I obtained might still not work as I came to know when I did numerical experiments today.
For the right solution form, we will have to look again at the simplest diffusion i.e the diffusion of brownian motion with drift. Previously I looked at brownian motion with centre at origin and that was very instructive and gave us the idea that FP equations can indeed be solved by this method.
Now we look at brownian motion that takes a starting value different from origin. I am not going to write the complete solution but it can only be solved in the form if we take squared difference from the mean.

Here is the second equation from top in yesterday's post.
${dw}[Z {(\frac{\partial w}{\partial Z})} +{(\frac{\partial w}{\partial Z})}^3 \frac{\partial Z^2}{\partial w^2}]$
$=\Big[-\mu_t(w)[Z {(\frac{\partial w}{\partial Z})} + {(\frac{\partial w}{\partial Z})}^3 \frac{\partial Z^2}{\partial w^2}] -\frac{\partial \mu_t(w)}{\partial w}{(\frac{\partial w}{\partial Z})}^2$
$+.5 {\sigma}^2 \big[ (Z^2-1)-3 Z{(\frac{\partial w}{\partial Z})}^2 {(\frac{\partial^2 Z}{\partial w^2})}+{(\frac{\partial w}{\partial Z})}^3 {(\frac{\partial Z^3}{\partial w^3})}] \Big] dt$

Please note that in following $\mu_{B}$ is mean of the density and $\mu_t$ is drift of the stochastic differential equation.
First thing to observe is that $Z {\frac{\partial w}{\partial Z}}$ on LHS of above equation which in case of brownian motion is equal to  $Z {\frac{\partial B}{\partial Z}}=\sigma \sqrt{t}Z =[B(t)-\mu_{B}]$.
And second observation is that there is no other possible way to solve for brownian motion that does not originate from zero (but has no time drift )than the following equations
$(B(t_1)-\mu_{B})^2 = (B(t_0)-\mu_{B})^2 + {\sigma}^2 Z^2 dt$
Please note that this very simple example perfectly agrees with Ito-hermite method where the above would be solved as
$B(t_1) = \mu_{B} + \sqrt{ (B(t_0)-\mu_{B})^2 + {\sigma}^2 Z^2 dt}$ and we had paid special care to signs.

And then we notice that Brownian motion with constand drift $\mu_t$ would have to be solved as
$(B(t_1)-\mu_{B}-\mu_t dt)^2 = (B(t_0)-\mu_{B})^2 + {\sigma}^2 Z^2 dt$

which perfectly agrees with out Ito-hermite method as the solution of brownian motion with drift there was
$B(t_1) = \mu_{B} +\mu_{t} dt + \sqrt{ (B(t_0)-\mu_{B})^2 + {\sigma}^2 Z^2 dt}$

Please note that these are initial informal notes that point towards the shape of the final formal solution but many of these things will slightly change in the coming days and in the final work. These are yet not the final equations or formal solution. There is a good possibility that $\mu_w$ in final formal solution will be slightly different and agree with the terms in Fokker-planck equation. These notes are only for motivation towards the formal solution.

Amin
Topic Author
Posts: 2421
Joined: July 14th, 2002, 3:00 am

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

In order to make a more informed attempt at the proper solution of fokker-planck equation, we look at the short term solution that comes from fourier transform.

Let us suppose we are given an SDE of the form
$dX(t)=\mu(X) dt + \sigma(X) dz(t)$
The short term solution to the associated fokker planck given by fourier transform is given as
refrence for the above is (  https://digital.library.adelaide.edu.au ... 02main.pdf ) : page 143
$P(t+\Delta t,y,t,x)=\frac{1}{\sqrt{4 Pi {\sigma(X)}^2 \Delta t }} \exp[-\frac{{(y-x-\mu(x) \Delta t +2 \frac{\partial {[\sigma(x)]}^2}{\partial x} \Delta t)}^2}{4{\sigma(x)}^2 \Delta t}]$
$* \exp[-\frac{\partial \mu(x)}{\partial x} \Delta t+ \frac{\partial^2 [{\sigma(x)}^2]}{\partial x^2}\Delta t]$

We can from the above formula easily infer which terms would have to go into a square-form and which terms would be added linearly. Every term inside the square on the above density would have to be added in a squared fashion with volatility. And the terms that are outside the square on second exponential would have to be added linearly. For people who are not familiar with previous work I did on Ito-hermite, please note that we are not in the above direct density framework. Informally you can know that by removing the normal density from the equation, we are dealing directly with terms on the exponential solution of FP equation, in a rough sense to take guidance about the proper formal solution in our framework.
We will have to be careful with expansion of second derivative term of fokker-planck (in Lamperti/Bessel form) since in the following term
${\sigma}^2 \Big[ Z^2-3 Z{(\frac{\partial w}{\partial Z})}^2 {(\frac{\partial^2 Z}{\partial w^2})}+{(\frac{\partial w}{\partial Z})}^3 {(\frac{\partial Z^3}{\partial w^3})}] \Big] dt$
it is only
${\sigma}^2 Z^2 dt$
that would probably appear in squared form  and other terms might have to be added later outside the squared expression. When a term in Fokker-planck of the SDE is expanded into multiple terms due to change of density derivatives, the above fourier transform guidance to include the terms inside squared complete differential expression will work for the first principal term only. just like in the above expansion of $.5 {\sigma}^2 \frac{d^2 P}{dw^2}$ the following terms are derived
${\sigma}^2 \Big[ Z^2-3 Z{(\frac{\partial w}{\partial Z})}^2 {(\frac{\partial^2 Z}{\partial w^2})}+{(\frac{\partial w}{\partial Z})}^3 {(\frac{\partial Z^3}{\partial w^3})}] \Big] dt$ while only ${\sigma}^2 Z^2 dt$ would probably appear in squared complete differential and rest of the terms will be manipulated and later added like our corrections in Ito-hermite algorithm. And I expect that case with drift is also similar.
I will come back with more experiments in another day or probably later today.

Amin
Topic Author
Posts: 2421
Joined: July 14th, 2002, 3:00 am

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

It seems that the terms that do not go inside the square and are not added in a squared fashion cannot be added directly. These terms have to be multiplied with
$\frac{1}{{(2 \pi)}^2}$ first and then added. I tried adding the following term
$-\frac{\partial \mu_t(w)}{\partial w}{(\frac{\partial w}{\partial Z})}^2$ after multiplying it with $\frac{1}{{(2 \pi)}^2}$
$-\frac{1}{{(2 \pi)}^2} \frac{\partial \mu_t(w)}{\partial w}{(\frac{\partial w}{\partial Z})}^2$
to the existing Ito-hermite program that I posted last time and it seemed to remarkably improve the mean without mean correction while preserving the original good shape of the density. Please note that previous corrections to Ito-hermite program were also multiplied by $\frac{1}{{(2 \pi)}^2}$
If you add any of the terms(that do not appear inside complete squared differential) directly in Ito-hermite density evolution, it will right away blow the program.
I am working on a new version of Ito-hermite density evolution program taking guidance from the Fokker-planck derivatives and hope to post it in a day.

Amin
Topic Author
Posts: 2421
Joined: July 14th, 2002, 3:00 am

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

In order to experiment with terms and to determine the exact solution of Fokker-planck, I tried to play with different terms that are not included in squared form. And without any detailed mathematical anlaysis, I tried some of those terms obtained from previous intermediate solution of fokker-planck equation. And I realized that simulation of densities of several SDEs especially improved in the tail and I got much better mean with mean-correction disabled. As I mentioned in previous post, every term that is not included in the squared form had to be scaled by $\frac{1}{{(2 \pi)}^2}$. I am including the program and I hope that it would be very helpful for other friends who are working on this or similar research.
function [TrueMean,ItoHermiteMean,MCMean] = ItoHermiteWilmottMeanCorrected080FP(x0,theta,kappa,gamma,sigma0,T)

%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%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.
%In this version, the mean correction has been disabled. If you are
%simulating standard mean-reverting SV type SDEs, please enable the mean
%correction by uncommenting the appropriate line in the body of code below.

dt=.125/16;   % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=128*T;%16*16;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
OrderA=4;  %
OrderM=4;  %
%dtM=.125/8;%Monte carlo time interval size dtM.
%TtM=8*8;%Monte carlo number of simulation intervals.
dtM=dt*1;
TtM=Tt/1;

dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=48;  % 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;

%theta=mu1/(-mu2);
%kappa=-mu2;

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

%you can specify any general mu1 and mu2 and beta1 and beta2.
%The present program is not perfect but quite reasonably good.
%if there is only one term in drift with beta=0, there will be minor
%errors. This case is yet not covered.

mu1=+1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
%sigma0=1.0;%Volatility value
alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;

w(1:Nn)=x0^(1-gamma)/(1-gamma);

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

for nn=2:Nn-1
ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1)-normcdf(.5*Z(nn)+.5*Z(nn-1),0,1);
%Above calculate probability mass in each probability subdivision.
end

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

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;%
d2wdZ2(1:Nn)=0;
dwdZ(1:Nn)=0;
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
tic

for tt=1:Tt

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

[wMu0dt,dwMu0dtdw,d2wMu0dtdw2,d3wMu0dtdw3,c1,c2,c3,c4] = CalculateDriftAndVolOrder4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);

dw(wnStart:Nn)=c1(wnStart:Nn).*Z(wnStart:Nn) ;% ...
dw2(wnStart:Nn)=dw(wnStart:Nn).^2;
% d2w(wnStart:Nn)=c2(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
% d3w(wnStart:Nn)=c3(wnStart:Nn).*(Z(wnStart:Nn).^3-3*Z(wnStart:Nn));
w(isnan(w)==1)=0;
wMu0dt(isnan(wMu0dt)==1)=0;
wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.

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

[w0,Sigma00,Sigma11,Sigma22,Sigma33,Sigma1,Sigma2,Sigma3,Sigma0] = CalculateHermiteRepresentationOfDensityOrderThreeNew(w,Z,wnStart,Nn,dNn,NnMidh,NnMidl);
d2wdZ2(wnStart:Nn)=d2wdZ2(wnStart:Nn)+ d2wMu0dtdw2(wnStart:Nn).*dwdZ(wnStart:Nn).^2+dwMu0dtdw(wnStart:Nn).*d2wdZ2(wnStart:Nn);%+ ...

dwdZ(wnStart:Nn)=(Sigma00)+Sigma1(wnStart:Nn).*2.*Z(wnStart:Nn);%.*ExpDamp(wnStart:Nn);

H2(wnStart:Nn)=(Z(wnStart:Nn).^2-1);
H3(wnStart:Nn)=(Z(wnStart:Nn).^3-3*Z(wnStart:Nn));
H4(wnStart:Nn)=(Z(wnStart:Nn).^4-6*Z(wnStart:Nn).^2+3);
CTerm(wnStart:Nn)=(.5*d2wMu0dtdw2(wnStart:Nn).*dwdZ(wnStart:Nn).^2 - ...
1*dwMu0dtdw(wnStart:Nn).*d2wdZ2(wnStart:Nn));

%Saved last
%Correction1(wnStart:Nn)=-1*sqrt(gamma)*(1+gamma*(1-gamma)*sigma0^2)*sqrt(8/8)/(22/7)^2.*CTerm(wnStart:Nn).*H2(wnStart:Nn);
%Correction0(wnStart:Nn)=-1*sqrt(gamma)*(1+(gamma)*(1-gamma)*sigma0^2)*sqrt(8/8)/(22/7)^2.*dt.*(d2wdZ2(wnStart:Nn)-.5*dwdZ(wnStart:Nn).^2).*H2(wnStart:Nn);

Correction1(wnStart:Nn)=-(1+gamma^2.*(1-gamma)*sigma0^2)/(22/7)^2.*CTerm(wnStart:Nn).*H2(wnStart:Nn);
Correction0(wnStart:Nn)=-(1+gamma^2.*(1-gamma)*sigma0^2)/(22/7)^2.*dt*(d2wdZ2(wnStart:Nn)-.5*dwdZ(wnStart:Nn).^2).*H2(wnStart:Nn);

if(tt==1)
w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+ ...
sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)));
end

%[dwdZ,d2wdZ2,d3wdZ3,d4wdZ4] = First4Derivatives2ndOrderEqSpaced(wnStart,Nn,dNn,w,Z);
dZdw(wnStart:Nn)=1.0./dwdZ(wnStart:Nn);
d2Zdw2(wnStart:Nn)=-1.*d2wdZ2(wnStart:Nn).*dZdw(wnStart:Nn).^3;
%d3Zdw3(1:Nn)=-d3wdZ3(1:Nn).*dZdw(1:Nn).^4+3*d2wdZ2(1:Nn).^2.*dZdw(1:Nn).^5;

if(tt>1)

fw0(wnStart:Nn)=exp(-.5*(1-gamma)*sigma0.^2.*Z(wnStart:Nn).^2);

%  dw(wnStart:Nn)=sigma0*sqrt(dt).*Z(wnStart:Nn);
%  dw2(wnStart:Nn)=sigma0^2*dt*Z(wnStart:Nn).^2;

w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+ ...
sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)))+ ...
+1/(2*pi)^2*dwMu0dtdw(wnStart:Nn).*dwdZ(wnStart:Nn).^2+ ...
-1/(2*pi)^2*wMu0dt(wnStart:Nn).*d2wdZ2(wnStart:Nn)+ ...
+1/(2*pi)^2*wMu0dt(wnStart:Nn).*Z(wnStart:Nn).*dwdZ(wnStart:Nn)+ ...
1*Correction1(wnStart:Nn).*fw0(wnStart:Nn) + ...
1*Correction0(wnStart:Nn).*fw0(wnStart:Nn) - ...
.5/(2*pi)^2*3*Z(wnStart:Nn).^1.*dZdw(wnStart:Nn).*d2wdZ2(wnStart:Nn)*sigma0^2.*dt;%.*fw0(wnStart:Nn);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Below are instability correction equations. There is instability in the
% SDE representation equation roughly between -.5 SD to +.5 SD in normal terms.
%Instability correction Algorithm takes two ends that are here given as
%lower end=w(NnMidl-2) and upper end=w(NnMidh+3).
%In the block below the interpolation starts from lower end.
%1st equation is
%And so on.
%There is a sixth equation that is not written below(since it takes us to
%the upper interpolation point which is already known) which is
%Sum of these six equations add up to the distance between upper and lower
%interpolation points.
%In the sum there are six instances of first derivative which is already
%calculated at the boundary, six instances of second derivative which is also
%already calculated at lower boundary and 18.0 instances of DeltaD1w which is
%unknown. This equation calculates the value of one unit of DeltaD1w addition
% I hope to find a better analytical method to tackle this issue. Friends
% are welcome to play with slightly changed related alternatives to the
%above interpolation.
CTheta=1.0;
D1wl=(11/6*w(NnMidl-2)-3*w(NnMidl-3)+3/2*w(NnMidl-4)-1/3*w(NnMidl-5))/dNn;
D2wl=(35/12*w(NnMidl-2) -26/3*w(NnMidl-3)+ 19/2*w(NnMidl-4)-14/3*w(NnMidl-5)+11/12*w(NnMidl-6))/dNn.^2;

%The last equation is not used since point w(NnMidh+3)is given but it is used in
%calculation of distance between two interpolation points. This distance is
%used to determine the value of DeltaD1w.

D1wh=(-11/6*w(NnMidh+3)+3*w(NnMidh+4)-3/2*w(NnMidh+5)+1/3*w(NnMidh+6))/dNn;
D2wh=(35/12*w(NnMidh+3) -26/3*w(NnMidh+4)+ 19/2*w(NnMidh+5)-14/3*w(NnMidh+6)+11/12*w(NnMidh+7))/dNn.^2;

%
w(NnMidh+2)=.5*whh(NnMidh+2)+.5*wll(NnMidh+2);
w(NnMidh+1)=.5*whh(NnMidh+1)+.5*wll(NnMidh+1);
w(NnMidh)=.5*whh(NnMidh)+.5*wll(NnMidh);
w(NnMidl)=.5*whh(NnMidl)+.5*wll(NnMidl);
w(NnMidl-1)=.5*whh(NnMidl-1)+.5*wll(NnMidl-1);

w1(1:Nn-1)=w(1:Nn-1);
w2(1:Nn-1)=w(2:Nn);

%         w1(1:Nn-2)=w(1:Nn-2);
%         w2(1:Nn-2)=w(2:Nn-1);
%         w3(1:Nn-2)=w(3:Nn);
%         w(w1(:)>w2(:))=w2(w1(:)>w2(:))-.25*(w3(w1(:)>w2(:))-w2(w1(:)>w2(:)));

%w(w1(:)>w2(:))=0;%Be careful;might not universally hold;

w(w<0)=0.0;
for nn=1:Nn
if(w(nn)<=0)
wnStart=nn+1;
end
end

%%1st order mean correction. We know the mean in original
%%coordinates so I shift the density into original coordinates,
%%apply the mean correction and then transform again into Lamperti
%%coordinates. Algebra solves two equations in two unknowns.
%%Equations are Sum_0^N{(Y_w(wnStart:Nn)-Y0)*W0}=1 and the second
%%equation is Sum_0^N{Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0}=u0
%%Two unknows are Y0 and W0. u0 is known mean.
tt;
u0=theta+(x0-theta)*exp(-kappa*(tt*dt)); %analytic mean of the density

%If you are not using stochastic volatility, replace above with
%true mean otherwise results would become garbage.

Y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));

Yn2Pn=sum(Y_w(wnStart:Nn).*Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
YnPn=sum(Y_w(wnStart:Nn).*ZProb(wnStart:Nn));
Pn=1.0;%%Sum(ZProb(1:Nn))
Y0=(Yn2Pn-u0*YnPn)/(YnPn-u0*Pn);

Y0Pn=Y0*Pn;
W0=1/(YnPn-Y0Pn);
YCorrected(wnStart:Nn)=Y_w(wnStart:Nn).*(Y_w(wnStart:Nn)-Y0)*W0;
wCorrected(wnStart:Nn)=(YCorrected(wnStart:Nn).^(1-gamma))./(1-gamma);
w(wnStart:Nn)=wCorrected(wnStart:Nn);

%I have disabled the mean correction. The above mean correction is only valid for
%standard mean reverting stochastic volatility type SDEs. To enable mean
%correction please uncomment the above last line in the block.

end

y_w(1:Nn)=0;
y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
Dfy_w(wnStart:Nn)=0;
Dfw(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
Dfy_w(nn) = (y_w(nn + 1) - y_w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
Dfw(nn) = (w(nn + 1) - w(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
fy_w(1:Nn)=0;
fw(1:Nn)=0;
for nn = wnStart:Nn-1
fy_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
fw(nn) = (normpdf(Z(nn),0, 1))/abs(Dfw(nn));
end

toc
%str=input('Analytic Values calculated; Press a key to start monte carlo');

rng(29079137, 'twister')
paths=100000;
YY(1:paths)=x0;  %Original process monte carlo.
Random1(1:paths)=0;
YYMean(1:TtM)=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);

%Uncomment for fourth order monte carlo

%
%   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,4,1).*YY(1:paths).^Fp(1,1,4,1)+YCoeff0(1,2,3,1).*YY(1:paths).^Fp(1,2,3,1)+ ...
%       YCoeff0(2,1,3,1).*YY(1:paths).^Fp(2,1,3,1)+YCoeff0(1,3,2,1).*YY(1:paths).^Fp(1,3,2,1)+ ...
%       YCoeff0(2,2,2,1).*YY(1:paths).^Fp(2,2,2,1)+YCoeff0(3,1,2,1).*YY(1:paths).^Fp(3,1,2,1)+ ...
%       YCoeff0(1,4,1,1).*YY(1:paths).^Fp(1,4,1,1)+YCoeff0(2,3,1,1).*YY(1:paths).^Fp(2,3,1,1)+ ...
%       YCoeff0(3,2,1,1).*YY(1:paths).^Fp(3,2,1,1)+YCoeff0(4,1,1,1).*YY(1:paths).^Fp(4,1,1,1))*dtM^3 + ...
%        (YCoeff0(1,1,5,1).*YY(1:paths).^Fp(1,1,5,1)+YCoeff0(1,2,4,1).*yy(1:paths).^Fp(1,2,4,1)+ ...
%     YCoeff0(2,1,4,1).*YY(1:paths).^Fp(2,1,4,1)+YCoeff0(2,1,1,1).*yy(1:paths).^Fp(2,1,1,1)+ ...
%      YCoeff0(2,2,3,1).*YY(1:paths).^Fp(2,2,3,1)+ ...
%      YCoeff0(3,1,3,1).*YY(1:paths).^Fp(3,1,3,1)+YCoeff0(1,4,2,1).*YY(1:paths).^Fp(1,4,2,1)+ ...
%      YCoeff0(2,3,2,1).*YY(1:paths).^Fp(2,3,2,1)+YCoeff0(3,2,2,1).*YY(1:paths).^Fp(3,2,2,1)+ ...
%      YCoeff0(4,1,2,1).*YY(1:paths).^Fp(4,1,2,1)+YCoeff0(1,5,1,1).*YY(1:paths).^Fp(1,5,1,1)+ ...
%      YCoeff0(2,4,1,1).*YY(1:paths).^Fp(2,4,1,1)+ ...
%      YCoeff0(3,3,1,1).*YY(1:paths).^Fp(3,3,1,1)+YCoeff0(4,2,1,1).*YY(1:paths).^Fp(4,2,1,1)+ ...
%      YCoeff0(5,1,1,1).*YY(1:paths).^Fp(5,1,1,1))*dtM^4+ ...
%       ((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+ ...
%     (YCoeff0(1,1,3,2).*YY(1:paths).^Fp(1,1,3,2)+YCoeff0(1,2,2,2).*YY(1:paths).^Fp(1,2,2,2)+ ...
%     YCoeff0(2,1,2,2).*YY(1:paths).^Fp(2,1,2,2)+YCoeff0(1,3,1,2).*YY(1:paths).^Fp(1,3,1,2)+ ...
%     YCoeff0(2,2,1,2).*YY(1:paths).^Fp(2,2,1,2)+YCoeff0(3,1,1,2).*YY(1:paths).^Fp(3,1,1,2)).*dtM^2.5+ ...
%     (YCoeff0(1,1,4,2).*YY(1:paths).^Fp(1,1,4,2)+YCoeff0(1,2,3,2).*YY(1:paths).^Fp(1,2,3,2)+ ...
%     YCoeff0(2,1,3,2).*YY(1:paths).^Fp(2,1,3,2)+YCoeff0(1,3,2,2).*YY(1:paths).^Fp(1,3,2,2)+ ...
%     YCoeff0(2,2,2,2).*YY(1:paths).^Fp(2,2,2,2)+ YCoeff0(3,1,2,2).*YY(1:paths).^Fp(3,1,2,2)+ ...
%     YCoeff0(1,4,1,2).*YY(1:paths).^Fp(1,4,1,2)+YCoeff0(2,3,1,2).*YY(1:paths).^Fp(2,3,1,2)+ ...
%     YCoeff0(3,2,1,2).*YY(1:paths).^Fp(3,2,1,2)+YCoeff0(4,1,1,2).*YY(1:paths).^Fp(4,1,1,2)).*dtM^3.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+ ...
%     (YCoeff0(1,1,3,3).*YY(1:paths).^Fp(1,1,3,3)+YCoeff0(1,2,2,3).*YY(1:paths).^Fp(1,2,2,3)+ ...
%     YCoeff0(2,1,2,3).*YY(1:paths).^Fp(2,1,2,3) + YCoeff0(1,3,1,3).*YY(1:paths).^Fp(1,3,1,3)+ ...
%     YCoeff0(2,2,1,3).*YY(1:paths).^Fp(2,2,1,3)+YCoeff0(3,1,1,3).*YY(1:paths).^Fp(3,1,1,3)).*dtM^3).*HermiteP1(3,1:paths) + ...
%     ((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )+ ...
%     (YCoeff0(1,1,2,4).*YY(1:paths).^Fp(1,1,2,4)+YCoeff0(1,2,1,4).*YY(1:paths).^Fp(1,2,1,4)+ ...
%     YCoeff0(2,1,1,4).*YY(1:paths).^Fp(2,1,1,4))*dtM^2.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);
%

YYMean(tt)=YYMean(tt)+sum(YY(1:paths))/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(y_w(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(500*gamma^2*1*sigma0);%Decrease the number of bins if the graph is too rough and increase if too straight lines.
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
plot(y_w(wnStart+1:Nn-1),fy_w(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('Ito-Hermite Method VS Monte Carlo Density, x0 = %.2f,theta=%.2f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f', x0,theta,kappa,gamma,sigma0,T));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

%title('text',sprintf('x0 = %f', x0),sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T))
%xlabel(a)
%title('Graphical Comparison of Ito-Hermite Method and Monte Carlo Method Density')
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

when you run the program with the following command,
ItoHermiteWilmottMeanCorrected080FP(1,1,1,.85,1,2)
here is the result from the program
ItoHermiteMean =0.999932367147095

One last thing about graphs in my program. When the monte carlo simulation density graphs are too rough, please decrease the variable NoOfBins  in the above program and when the graphs are too straight lines, increase the variable NoOfBins. I have also indicated it in the program.

Amin
Topic Author
Posts: 2421
Joined: July 14th, 2002, 3:00 am

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

Sorry friends, a quick update about yesterday's post. I believe that multiplier in terms that are not in the squared form is not $\frac{1}{{(2 \pi)}^2}$. It should rather probably be another  $4 \Delta t$  on top of first time integration. Here is the main update equation changed so can just plug it in the previous program and run it.

C_f=4*dt;
w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+ ...
sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)))+ ...
+C_f*dwMu0dtdw(wnStart:Nn).*dwdZ(wnStart:Nn).^2+ ...
-C_f*wMu0dt(wnStart:Nn).*d2wdZ2(wnStart:Nn)+ ...
+C_f*wMu0dt(wnStart:Nn).*Z(wnStart:Nn).*dwdZ(wnStart:Nn)+ ...
1*Correction1(wnStart:Nn).*fw0(wnStart:Nn) + ...
1*Correction0(wnStart:Nn).*fw0(wnStart:Nn) - ...
C_f*.5*3*Z(wnStart:Nn).^1.*dZdw(wnStart:Nn).*d2wdZ2(wnStart:Nn)*sigma0^2.*dt;

also attached as code
C_f=4*dt;
w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+ ...
sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).* ...
sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...
sign(dw(wnStart:Nn)).*dw2(wnStart:Nn)))+ ...
+C_f*dwMu0dtdw(wnStart:Nn).*dwdZ(wnStart:Nn).^2+ ...
-C_f*wMu0dt(wnStart:Nn).*d2wdZ2(wnStart:Nn)+ ...
+C_f*wMu0dt(wnStart:Nn).*Z(wnStart:Nn).*dwdZ(wnStart:Nn)+ ...
1*Correction1(wnStart:Nn).*fw0(wnStart:Nn) + ...
1*Correction0(wnStart:Nn).*fw0(wnStart:Nn) - ...
C_f*.5*3*Z(wnStart:Nn).^1.*dZdw(wnStart:Nn).*d2wdZ2(wnStart:Nn)*sigma0^2.*dt;

Amin
Topic Author
Posts: 2421
Joined: July 14th, 2002, 3:00 am

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

Friends, I had been working with the older program and tried to improve it and came back to a new version of the density simulation program based totally on fokker-planck derivatives. I completed the square for the solution of FP equation in bessel/Lamperti coordinates and I am trying to work out all the details of the numerical algorithm based totally on the derivatives derived from the solution of the FP equations and hope to post the new program in another day or two if things go well.

Amin
Topic Author
Posts: 2421
Joined: July 14th, 2002, 3:00 am

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

Though densities are different and more complicated, consider this following simple expression and its application to more complex SDEs.
If $B_t$ is brownian motion with volatility $\sigma$ the interesting expression is

$d\big[ \frac{{(B_t - B_0)}^2}{{\sigma}^2 \Delta t} - Z^2 \big] =0$

and consider its applications to more complicated SDEs.

Amin
Topic Author
Posts: 2421
Joined: July 14th, 2002, 3:00 am

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

I have previously asked friends on this forum to let me take this work on fokker-planck and densities of SDEs towards conclusion but I think there are many friends who can do this work far better than I can. And therefore I will ask those friends who are experts in this domain or have strong interest in this area to openly share and publish any related research and not wait for anything like a formal research paper from me as I had earlier asked. I would still request friends to please give me credit for whatever work I have presented on Wilmott forum. I will continue to do further research, keep posting it on Wilmott forum and might even write a paper but I would still want friends to go ahead with their research and not wait for me to present any formal research work.

Amin
Topic Author
Posts: 2421
Joined: July 14th, 2002, 3:00 am

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

Dear friends I made some notes about simulation of SDEs using Ito-hermite method. In order to simulate the SDEs with this method, first we have to calculate the complete differential of square of terms as they appear inside the normal exponential. I will start with a general SDE in Lamperti form. Here is the form of the SDE
$dw_s=\mu(w_s) ds + \sigma dZ(s)$   Equation(1)

Let us suppose$\zeta_w$ is the mean of density when we want to calculate the one time step ahead value of $w_t$ to calculate  the density.
First of all we want to calculate the differential and later we will go to fokker-planck

$\int_0^t d\Big[ {(w_s -\zeta_w - \int_0^s \mu(w_v) dv )}^2 \Big] = {\sigma}^2 Z^2 \Delta t$     Equation(2)

Let us suppose that we indicate the integral of drift by uppercase letter as
$M_t=\int_0^t \mu(w_s) ds$   Equation(3)

Expanding the left side in the equation(2), we get

$\int_0^t d\Big[ {({w_s}^2+ {\zeta_w}^2 + {M_s}^2 - 2 w_s \zeta_w -2 w_s M_s + 2 \zeta_w M_s )} \Big]$

applying the differential we get

$\int_0^t \Big[ 2w_s dw_s +{\sigma}^2 ds +2 M_s dM_s - 2 \zeta_w dw_s -2 d[w_s M_s] + 2 \zeta_w dM_s \Big]$
The first integral given as
$\int_0^t \Big[ 2w_s dw_s +{\sigma}^2 ds \Big]$ is straightforward. The only two integrals that are relatively difficult are calculated below

$\int_0^t M_s dM_s$
In order to calculate the above integral, we first calculate
$M_s=\int_0^s \mu(w_v) dv$
$dM_s= \mu(w_s) ds$
$M_s= \mu(w_0) \int_0^s dv + \frac{\partial \mu(w_0)}{\partial w} \mu(w_0) \int_0^s \int_0^v du dv + \frac{\partial \mu(w_0)}{\partial w} \sigma \int_0^s \int_0^v dz(u) dv + .5 \frac{\partial^2 \mu(w_0)}{\partial w^2} {\sigma}^2 \int_0^s \int_0^v du dv$
$=\mu(w_0) s + \frac{\partial \mu(w_0)}{\partial w} \mu(w_0) \frac{s^2}{2} + \frac{\partial \mu(w_0)}{\partial w} \sigma (1-\frac{1}{\sqrt{3}}) s z(s)$
$+ .5 \frac{\partial^2 \mu(w_0)}{\partial w^2} {\sigma}^2 \frac{s^2}{2}$

Now
$\int_0^t M_s dM_s= \int_0^t\Big[ \mu(w_0) s + \frac{\partial \mu(w_0)}{\partial w} \mu(w_0) \frac{s^2}{2} + \frac{\partial \mu(w_0)}{\partial w} \sigma (1-\frac{1}{\sqrt{3}}) s z(s)$
$+ .5 \frac{\partial^2 \mu(w_0)}{\partial w^2} {\sigma}^2 \frac{s^2}{2} \Big ] dMs$
$= \int_0^t\Big[ {\mu(w_0)}^2 s ds+ \frac{\partial \mu(w_0)}{\partial w} {\mu(w_0)}^2 \frac{s^2}{2} ds+ \frac{\partial \mu(w_0)}{\partial w} \mu(w_0) \sigma (1-\frac{1}{\sqrt{3}}) s z(s) ds$
$+ .5 \frac{\partial^2 \mu(w_0)}{\partial w^2} \mu(w_0) {\sigma}^2 \frac{s^2}{2} ds \Big ]$
$= {\mu(w_0)}^2 \frac{t^2}{2}+ \frac{\partial \mu(w_0)}{\partial w} {\mu(w_0)}^2 \frac{t^3}{6}+ \frac{\partial \mu(w_0)}{\partial w} \mu(w_0) \sigma (1-\frac{1}{\sqrt{3}}) \frac{1}{2} (1-\frac{1}{\sqrt{5}}) t^2 z(s)$
$+ .5 \frac{\partial^2 \mu(w_0)}{\partial w^2} \mu(w_0) {\sigma}^2 \frac{t^3}{6}$

The other difficult and more important integral is calculated as
$= \int_0^t d[w_s M_s]= \int_0^t M_s dw_s +w_s dM_s$
We substitute the value of $M_s$ from above equation No.  into the above equation as
$= \int_0^t \Big[ \big[ \mu(w_0) s + \frac{\partial \mu(w_0)}{\partial w} \mu(w_0) \frac{s^2}{2} + \frac{\partial \mu(w_0)}{\partial w} \sigma (1-\frac{1}{\sqrt{3}}) s z(s)$
$+ .5 \frac{\partial^2 \mu(w_0)}{\partial w^2} {\sigma}^2 \frac{s^2}{2} \big] \big[\mu(w_0) ds + \sigma dz(s) \big]$
$+w_0 \mu(w_0) ds + w_0 \frac{\partial \mu(w_0)}{\partial w} \mu(w_0) s ds + w_0 \frac{\partial \mu(w_0)}{\partial w} \sigma z(s) ds$
$+w_0 .5 \frac{\partial^2 \mu(w_0)}{\partial w^2} {\sigma}^2 s ds \Big]$
$= \Big[ \big[ {\mu(w_0)}^2 \frac {t^2}{2} + \frac{\partial \mu(w_0)}{\partial w} {\mu(w_0)}^2 \frac{t^3}{6}$
$+\frac{\partial \mu(w_0)}{\partial w} \sigma (1-\frac{1}{\sqrt{3}}) \frac{1}{2} (1-\frac{1}{\sqrt{5}}) t^2 z(t)$
$+ .5 \frac{\partial^2 \mu(w_0)}{\partial w^2} \mu(w_0) {\sigma}^2 \frac{t^3}{6} \big]$
$+ \big[ {\mu(w_0)} \sigma (1-\frac {1}{\sqrt{3}}) t z(t) + \frac{\partial \mu(w_0)}{\partial w} {\mu(w_0)} .5 t^2 z(t)/\sqrt{2}/\sqrt{6}$
$+\frac{\partial \mu(w_0)}{\partial w} \sigma (1-\frac{1}{\sqrt{3}}) \frac{t ({z(t)}^2-1)}{2\sqrt{2}}$
$+ .5 \frac{\partial^2 \mu(w_0)}{\partial w^2} \mu(w_0) {\sigma}^2 \frac{t^2 z(t)}{\sqrt{2}}{\sqrt{6}} \big]$
$+w_0 \mu(w_0) t + w_0 \frac{\partial \mu(w_0)}{\partial w} \mu(w_0) \frac{t^2}{2} + w_0 \frac{\partial \mu(w_0)}{\partial w} \sigma (1-\frac{1}{\sqrt{3}}) t z(t)$
$+w_0 .5 \frac{\partial^2 \mu(w_0)}{\partial w^2} {\sigma}^2 \frac{t^2}{2} \Big]$

Amin
Topic Author
Posts: 2421
Joined: July 14th, 2002, 3:00 am

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

Dear friends, when I wrote the post last night, it was already very late and by the time I completed the post I was barely awake and I had already slept briefly multiple times on my desk so I was unable to correct every mistake. In the expansion of second term in the integral  $\int_0^t d[w_s M_s]= \int_0^t \big[ M_s dw_s +w_s dM_s \big]$ , I did not properly expand $w_s$ and just used $w_0$ in its place and a proper expansion would require replacing $w_s$ by
$w_0+ \mu(w_0) s + \sigma dz(s)$  and then solving the resulting integrals after multiplication with $dM_s$
This would not be important to mention but there are some second order terms that we would miss if we do not expand this integral properly.

The other difficult and more important integral is calculated as
$= \int_0^t d[w_s M_s]= \int_0^t M_s dw_s +w_s dM_s$
We substitute the value of $M_s$ from above equation No.  into the above equation as
$= \int_0^t \Big[ \big[ \mu(w_0) s + \frac{\partial \mu(w_0)}{\partial w} \mu(w_0) \frac{s^2}{2} + \frac{\partial \mu(w_0)}{\partial w} \sigma (1-\frac{1}{\sqrt{3}}) s z(s)$
$+ .5 \frac{\partial^2 \mu(w_0)}{\partial w^2} {\sigma}^2 \frac{s^2}{2} \big] \big[\mu(w_0) ds + \sigma dz(s) \big]$
$+w_0 \mu(w_0) ds + w_0 \frac{\partial \mu(w_0)}{\partial w} \mu(w_0) s ds + w_0 \frac{\partial \mu(w_0)}{\partial w} \sigma z(s) ds$
$+w_0 .5 \frac{\partial^2 \mu(w_0)}{\partial w^2} {\sigma}^2 s ds \Big]$
$= \Big[ \big[ {\mu(w_0)}^2 \frac {t^2}{2} + \frac{\partial \mu(w_0)}{\partial w} {\mu(w_0)}^2 \frac{t^3}{6}$
$+\frac{\partial \mu(w_0)}{\partial w} \sigma (1-\frac{1}{\sqrt{3}}) \frac{1}{2} (1-\frac{1}{\sqrt{5}}) t^2 z(t)$
$+ .5 \frac{\partial^2 \mu(w_0)}{\partial w^2} \mu(w_0) {\sigma}^2 \frac{t^3}{6} \big]$
$+ \big[ {\mu(w_0)} \sigma (1-\frac {1}{\sqrt{3}}) t z(t) + \frac{\partial \mu(w_0)}{\partial w} {\mu(w_0)} .5 t^2 z(t)/\sqrt{2}/\sqrt{6}$
$+\frac{\partial \mu(w_0)}{\partial w} \sigma (1-\frac{1}{\sqrt{3}}) \frac{t ({z(t)}^2-1)}{2\sqrt{2}}$
$+ .5 \frac{\partial^2 \mu(w_0)}{\partial w^2} \mu(w_0) {\sigma}^2 \frac{t^2 z(t)}{\sqrt{2}}{\sqrt{6}} \big]$
$+w_0 \mu(w_0) t + w_0 \frac{\partial \mu(w_0)}{\partial w} \mu(w_0) \frac{t^2}{2} + w_0 \frac{\partial \mu(w_0)}{\partial w} \sigma (1-\frac{1}{\sqrt{3}}) t z(t)$
$+w_0 .5 \frac{\partial^2 \mu(w_0)}{\partial w^2} {\sigma}^2 \frac{t^2}{2} \Big]$