Serving the Quantitative Finance Community

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

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

Friends, I want to explain the series solution method for constant CDF lines of the density and relate it to my program and then discuss improvements I have made.
Previously we had derived the differential equation of constant CDF lines of solution of an SDE from Fokker-planck equation of the SDE. Here (Post 1202) are the details of that derivation:  https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1200#p867718

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

From method of iterated integrals for ordinary differential equations, For reference see: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2872598
we know that for a first order differential equation of the form
$\frac{dw}{dt}=f(w)$ ,

the solution to second order accuracy in time is given as
$w(t_1) \, =\, w(t_0)+ f(w(t_0)) \, \int_{t_0}^{t_1} \, dt +\, f'(w(t_0)) f(w(t_0)) \, \int_{t_0}^{t_1} \int_{t_0}^t \, ds \, dt \,$

or
$w(t_1) \, =\, w(t_0)+ f(w(t_0)) \, (t_1-t_0)+\, f'(w(t_0)) f(w(t_0)) \,\frac{ {(t_1-t_0)}^2}{2} \,$

incidentally after the next higher order term, the third order accurate in time solution can be written as
$w(t_1) \, =\, w(t_0)+ f(w(t_0)) \, (t_1-t_0)+\, f'(w(t_0)) f(w(t_0)) \,\frac{ {(t_1-t_0)}^2}{2} \,+\, \Big[ f''(w(t_0)) {f(w(t_0))}^2 + {f'(w(t_0))}^2 f(w(t_0)) \Big] \,\frac{ {(t_1-t_0)}^3}{6} \,$

But we have not added third order term in our series solution to constant CDF lines solution to Fokker-planck equation yet. In our case
$\frac{dw}{dt}=\mu(w) \, + \, .5 {\sigma}^2 Z \frac{\partial Z}{\partial w} - \, .5 {\sigma}^2 \frac{\partial^2 Z}{\partial w^2} [\frac{\partial Z}{\partial w}]^{-1}$
So
$f(w)=\mu(w) \, + \, .5 {\sigma}^2 Z \frac{\partial Z}{\partial w} - \, .5 {\sigma}^2 \frac{\partial^2 Z}{\partial w^2} [\frac{\partial Z}{\partial w}]^{-1}$
and
$f'(w)=\mu'(w) \, + \, .5 {\sigma}^2 \Big[{ (\frac{\partial Z}{\partial w})}^2 \,+Z \frac{\partial^2 Z}{\partial w^2}\, \Big] - \, .5 {\sigma}^2 \Big[ \frac{\partial^3 Z}{\partial w^3} [\frac{\partial Z}{\partial w}]^{-1} - \, {(\frac{\partial^2 Z}{\partial w^2})}^2 [\frac{\partial Z}{\partial w}]^{-2} \, \Big]$

Our second order solution to differential equation of constant CDF lines is given as
$w(t_1) \, =\, w(t_0)+ f(w(t_0)) \, \Delta t+\, f'(w(t_0)) f(w(t_0)) \,\frac{ {\Delta t}^2}{2} \,$
which becomes
$w(t_1) \, =\, w(t_0)+ \Bigg[ \mu(w) \, + \, .5 {\sigma}^2 Z \frac{\partial Z}{\partial w} - \, .5 {\sigma}^2 \frac{\partial^2 Z}{\partial w^2} [\frac{\partial Z}{\partial w}]^{-1}\Bigg] \, \Delta t \,$
$+\, \Bigg[ \mu'(w) \, + \, .5 {\sigma}^2 \Big[{ (\frac{\partial Z}{\partial w})}^2 \,+Z \frac{\partial^2 Z}{\partial w^2}\, \Big] - \, .5 {\sigma}^2 \Big[ \frac{\partial^3 Z}{\partial w^3} [\frac{\partial Z}{\partial w}]^{-1} - \, {(\frac{\partial^2 Z}{\partial w^2})}^2 [\frac{\partial Z}{\partial w}]^{-2} \, \Big]\Bigg] \, \Bigg[\mu(w) \, + \, .5 {\sigma}^2 Z \frac{\partial Z}{\partial w} - \, .5 {\sigma}^2 \frac{\partial^2 Z}{\partial w^2} [\frac{\partial Z}{\partial w}]^{-1} \Bigg] \,\frac{ {\Delta t}^2}{2} \,$

If wanted, we could extend the solution to third order in time for well behaved SDEs using the third order expansion of the differential equation that I described earlier.
In the next post, I will describe the series solution and then follow with another post to explain the program.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

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

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

To find the solution of density of SDE along constant CDF (or Constant Z-lines) lines,  we suppose that we can express the solution of w(t) where w is the SDE/PDE variable as a power series in Z. After all the research we have done in past few months, it seems very natural that we express  w(t) as a series in Z for any given time along time evolution of the density. We find the coordinates of this series at median of the density where Z=0. If we can represent w(t) at any given time as a valid series in Z, we can find all its derivatives and other smooth functions and also represent them as a series in Z. For a given evolution equation of the constant CDF lines of the density, we can represent all the variables in right hand side as a series in Z. The products, reciprocals and integrals of theses series again result as answers in another series in Z. So we find power series in Z for all variables in the right hand side of evolution equation and then do all mathematical operations on them in the form of series. The resulting answer is based on our original series at time $t_0$ but generally  has different series coefficients than that of series at time $t_0$. This new series that follows the application of SDE evolution equation at series coefficients at time $t_0$ represents the solution of constant CDF lines of SDE density $w(t_1)$ at next time level $t_1$ as a power series in Z again. Now we take this power series at time $t_1$ and find the series representations of its derivatives and other smooth functions and substitute them again in SDE evolution equation to find the series representation at time $t_2$.
I will write the proposed power series for w(t) and its derivatives and then explain how to find the power series of smooth functions of w when power series for w is given. The series for w(t) is given as
$w(t)=a_0 \,+ a_1 \, Z + a_2 \, Z^2 + a_3 \, Z^3 + a_4 \, Z^4 +\, ...$

The first, second and third  derivatives of w(t) w.r.t Z are given as
$\frac{\partial w(t)}{\partial Z}= a_1 \, + 2 \,a_2 \, Z +3 \, a_3 \, Z^2 +4 \, a_4 \, Z^3 +\, ...$
$\frac{\partial^2 w(t)}{\partial Z^2}= 2 \,a_2 \, +6 \, a_3 \, Z +12 \, a_4 \, Z^2 +\, ...$
$\frac{\partial^3 w(t)}{\partial Z^3}= 6 \, a_3 \, +24 \, a_4 \, Z \,+60 \, a_5 \, Z^2 +\, ...$

We do all the calculations of series evolution equation at median where Z=0. here all the derivatives are given only by leading coefficient associated with zeroth power of Z. So at Z=0,

$w(t)=a_0 \,$
$\frac{\partial w(t)}{\partial Z}= a_1 \,$
$\frac{\partial^2 w(t)}{\partial Z^2}= 2 \,a_2 \,$
$\frac{\partial^3 w(t)}{\partial Z^3}= 6 \, a_3 \,$

In order to represent smooth functions of w(t) as a series in Z, we suppose a form of the series and equate its coefficients with analytic derivatives w.r.t Z and this calculation is done at median that is Z=0.
So let us take drift in SDE, $\mu(w)$ as a smooth function of w and represent it as a power series in Z. In what follows $\mu(w)$ is considered a function operating on w just like we use symbol $f(w)$ to indicate that f is a function operating on w. We suppose that functional form of $\mu(w)$ and its derivatives is known analytically.
Let us suppose we know the form of power series of $\mu(w)$ with the coefficients at yet unknown. This power series is
$\mu(w(t))=b_0 \,+ b_1 \, Z + b_2 \, Z^2 + b_3 \, Z^3 + b_4 \, Z^4 +\, ...$

The first, second and third  derivatives of $\mu(w(t))$ w.r.t Z are given as
$\frac{\partial \mu(w(t))}{\partial Z}= b_1 \, + 2 \,b_2 \, Z +3 \, b_3 \, Z^2 +4 \, b_4 \, Z^3 +\, ...$
$\frac{\partial^2 \mu(w(t))}{\partial Z^2}= 2 \,b_2 \, +6 \, b_3 \, Z +12 \, b_4 \, Z^2 +\, ...$
$\frac{\partial^3 \mu(w(t))}{\partial Z^3}= 6 \, b_3 \, +24 \, b4 \, Z \,+60 \, b_5 \, Z^2 +\, ...$

We know that values of above drift and its derivatives are given by leading coefficients of the series when Z=0. So at Z=0,

$\mu(w(t))=b_0 \,$
$\frac{\partial \mu(w(t))}{\partial Z}= b_1 \,$
$\frac{\partial^2 \mu(w(t))}{\partial Z^2}= 2 \,b_2 \,$
$\frac{\partial^3 \mu(w(t))}{\partial Z^3}= 6 \, b_3 \,$

We also know the value of $w(t)$ at median which is $a_0$
but $\mu(w(t))=b_0 \,$  which means that $b_0=\mu(a_0) \,$
To find the second coefficient, we equate first derivative of drift w.r.t Z at $w_0=a_0$ with its coefficient. We know that
$\frac{\partial \mu(w(t))}{\partial Z} \,= \, \frac{\partial \mu(w(t))}{\partial w} \,\frac{\partial w(t)}{\partial Z} \,$ which equals
$=\, \frac{\partial \mu(w(t))}{\partial w} \, a_1$
but $\frac{\partial \mu(w(t))}{\partial Z}= b_1 \,$
so
$b_1=\, \frac{\partial \mu(w(t))}{\partial w} \, a_1=\, \frac{\partial \mu(a_0)}{\partial w} \, a_1$

Similarly
$\frac{\partial^2 \mu(w(t))}{\partial Z^2} \,= \, \frac{\partial^2 \mu(w(t))}{\partial w^2} \, {(\frac{\partial w(t)}{\partial Z} )}^2+ \, \, \frac{\partial \mu(w(t))}{\partial w} \,\frac{\partial^2 w(t)}{\partial Z^2} \,$
$= \, \frac{\partial^2 \mu(w(t))}{\partial w^2} {a_1}^2 \, + \, \, \, \frac{\partial \mu(w(t))}{\partial w} (2a_2)$

We equate this with   $\frac{\partial^2 \mu(w(t))}{\partial Z^2}= 2 \,b_2 \,$  to find out that
$b_2= 1/2 \, \Big[\, \frac{\partial^2 \mu(w(t))}{\partial w^2} {a_1}^2 \, + \, \, \, \frac{\partial \mu(w(t))}{\partial w} (2a_2) \Big]$
since above values of derivatives of w are calculated at $w=a_0$, we can write
$b_2= 1/2 \, \Big[\, \frac{\partial^2 \mu(a_0)}{\partial w^2} {a_1}^2 \, + \, \, \, \frac{\partial \mu(a_0)}{\partial w} (2a_2) \Big]$

similarly equating third derivative we get
$6 b_3= \, \frac{\partial^3 \mu(w)}{\partial w^3} {(\frac{\partial w(t)}{\partial Z} )}^3 \, + \,3 \, \frac{\partial^2 \mu(w)}{\partial w^2} \frac{\partial w(t)}{\partial Z} \frac{\partial^2 w(t)}{\partial Z^2} \, + \, \frac{\partial \mu(w)}{\partial w} \frac{\partial^3 w(t)}{\partial Z^3} \,$

which means
$b_3= 1/6 \, \Big[\, \frac{\partial^3 \mu(a_0)}{\partial w^3} {a_1}^3 \, + \,6 \, \frac{\partial^2 \mu(a_0)}{\partial w^2}\, a_1 \, a_2 + \, \frac{\partial \mu(a_0)}{\partial w} \, 6a_3 \, \Big]$

So this way we continue to find the power series of smooth functions of w by equating the series coefficients for derivatives at Z=0 with their functional form there and using chain rule of derivatives.

I am very sure that similar series method can be used for a lot of other partial differential equations where solution remains smooth and well-behaved and series representation converges well.

In next relatively non-technical post, I will just give an idea how to relate to variables I have used in my matlab implementation to equations in this and previous post. In another later post, I will explain some enhancements I have used to get the series solution method to get to work better closer to zero.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

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

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

Friends, I want to explain the series solution method for constant CDF lines of the density and relate it to my program and then discuss improvements I have made.
Previously we had derived the differential equation of constant CDF lines of solution of an SDE from Fokker-planck equation of the SDE. Here (Post 1202) are the details of that derivation:  https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1200#p867718

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

From method of iterated integrals for ordinary differential equations, For reference see: https://papers.ssrn.com/sol3/papers.cfm?abstract_id=2872598
we know that for a first order differential equation of the form
$\frac{dw}{dt}=f(w)$ ,

the solution to second order accuracy in time is given as
$w(t_1) \, =\, w(t_0)+ f(w(t_0)) \, \int_{t_0}^{t_1} \, dt +\, f'(w(t_0)) f(w(t_0)) \, \int_{t_0}^{t_1} \int_{t_0}^t \, ds \, dt \,$

or
$w(t_1) \, =\, w(t_0)+ f(w(t_0)) \, (t_1-t_0)+\, f'(w(t_0)) f(w(t_0)) \,\frac{ {(t_1-t_0)}^2}{2} \,$

incidentally after the next higher order term, the third order accurate in time solution can be written as
$w(t_1) \, =\, w(t_0)+ f(w(t_0)) \, (t_1-t_0)+\, f'(w(t_0)) f(w(t_0)) \,\frac{ {(t_1-t_0)}^2}{2} \,+\, \Big[ f''(w(t_0)) {f(w(t_0))}^2 + {f'(w(t_0))}^2 f(w(t_0)) \Big] \,\frac{ {(t_1-t_0)}^3}{6} \,$

But we have not added third order term in our series solution to constant CDF lines solution to Fokker-planck equation yet. In our case
$\frac{dw}{dt}=\mu(w) \, + \, .5 {\sigma}^2 Z \frac{\partial Z}{\partial w} - \, .5 {\sigma}^2 \frac{\partial^2 Z}{\partial w^2} [\frac{\partial Z}{\partial w}]^{-1}$
So
$f(w)=\mu(w) \, + \, .5 {\sigma}^2 Z \frac{\partial Z}{\partial w} - \, .5 {\sigma}^2 \frac{\partial^2 Z}{\partial w^2} [\frac{\partial Z}{\partial w}]^{-1}$
and
$f'(w)=\mu'(w) \, + \, .5 {\sigma}^2 \Big[{ (\frac{\partial Z}{\partial w})}^2 \,+Z \frac{\partial^2 Z}{\partial w^2}\, \Big] - \, .5 {\sigma}^2 \Big[ \frac{\partial^3 Z}{\partial w^3} [\frac{\partial Z}{\partial w}]^{-1} - \, {(\frac{\partial^2 Z}{\partial w^2})}^2 [\frac{\partial Z}{\partial w}]^{-2} \, \Big]$

Our second order solution to differential equation of constant CDF lines is given as
$w(t_1) \, =\, w(t_0)+ f(w(t_0)) \, \Delta t+\, f'(w(t_0)) f(w(t_0)) \,\frac{ {\Delta t}^2}{2} \,$
which becomes
$w(t_1) \, =\, w(t_0)+ \Bigg[ \mu(w) \, + \, .5 {\sigma}^2 Z \frac{\partial Z}{\partial w} - \, .5 {\sigma}^2 \frac{\partial^2 Z}{\partial w^2} [\frac{\partial Z}{\partial w}]^{-1}\Bigg] \, \Delta t \,$
$+\, \Bigg[ \mu'(w) \, + \, .5 {\sigma}^2 \Big[{ (\frac{\partial Z}{\partial w})}^2 \,+Z \frac{\partial^2 Z}{\partial w^2}\, \Big] - \, .5 {\sigma}^2 \Big[ \frac{\partial^3 Z}{\partial w^3} [\frac{\partial Z}{\partial w}]^{-1} - \, {(\frac{\partial^2 Z}{\partial w^2})}^2 [\frac{\partial Z}{\partial w}]^{-2} \, \Big]\Bigg] \, \Bigg[\mu(w) \, + \, .5 {\sigma}^2 Z \frac{\partial Z}{\partial w} - \, .5 {\sigma}^2 \frac{\partial^2 Z}{\partial w^2} [\frac{\partial Z}{\partial w}]^{-1} \Bigg] \,\frac{ {\Delta t}^2}{2} \,$

If wanted, we could extend the solution to third order in time for well behaved SDEs using the third order expansion of the differential equation that I described earlier.
In the next post, I will describe the series solution and then follow with another post to explain the program.
.
Friends this is not a very technical post and what I am trying to do here is pretty much obvious but I still decided to do some more algebra to bring the equations in copied post in line with equations used in the program.
We have that in copied post
$f(w)=\mu(w) \, + \, .5 {\sigma}^2 Z \frac{\partial Z}{\partial w} - \, .5 {\sigma}^2 \frac{\partial^2 Z}{\partial w^2} [\frac{\partial Z}{\partial w}]^{-1}$
and
$f'(w)=\mu'(w) \, + \, .5 {\sigma}^2 \Big[{ (\frac{\partial Z}{\partial w})}^2 \,+Z \frac{\partial^2 Z}{\partial w^2}\, \Big] - \, .5 {\sigma}^2 \Big[ \frac{\partial^3 Z}{\partial w^3} [\frac{\partial Z}{\partial w}]^{-1} - \, {(\frac{\partial^2 Z}{\partial w^2})}^2 [\frac{\partial Z}{\partial w}]^{-2} \, \Big]$

These equations use the form of derivatives with respect to w while our series derivatives are with respect to Z. We use the following relationships

$\frac{\partial Z}{\partial w}= {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-1}$
We can perform the above operation by series reciprocal function given with my matlab program. Higher derivatives cannot be inverted like that and we have to use relationships below
$\frac{\partial^2 Z}{\partial w^2}=- \, \frac{\partial^2 w}{\partial Z^2} {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-3}$
and
$\frac{\partial^3 Z}{\partial w^3}=- \, \frac{\partial^3 w}{\partial Z^3} {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-4}+ 3 \, {(\frac{\partial^2 w}{\partial Z^2})}^2 \, {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-5}$
So all higher derivatives are with respect to Z and first derivative and its negative powers can be found by series reciprocal program(which gives another power series for reciprocal of a variable).
Now we substitute the above relationships in $f(w)$  and $f'(w)$ to bring them in line with expressions used in my matlab program.
Now
$f(w)=\mu(w) \, + \, .5 {\sigma}^2 Z \frac{\partial Z}{\partial w} - \, .5 {\sigma}^2 \Big[- \, \frac{\partial^2 w}{\partial Z^2} {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-3} \, \Big] [\frac{\partial Z}{\partial w}]^{-1}$
$=\mu(w) \, + \, .5 {\sigma}^2 Z \big[\frac{\partial w}{\partial Z}\big]^{-1} + \, .5 {\sigma}^2 \, \frac{\partial^2 w}{\partial Z^2} {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-2} \,$
while
$f'(w)=\mu'(w) \, + \, .5 {\sigma}^2 \Bigg[{ (\frac{\partial w}{\partial Z})}^{-2} \,+Z \Big[- \, \frac{\partial^2 w}{\partial Z^2}\,{\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-3} \,\Big] \Bigg]$
$- \, .5 {\sigma}^2 \Bigg[ \, \Big[- \, \frac{\partial^3 w}{\partial Z^3} {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-4}+ 3 \, {(\frac{\partial^2 w}{\partial Z^2})}^2 \, {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-5} \Big][\frac{\partial Z}{\partial w}]^{-1} - \, \Big[{(- \, \frac{\partial^2 w}{\partial Z^2} {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-3})}^2 \Big]\, [\frac{\partial Z}{\partial w}]^{-2} \, \Bigg]$
$=\mu'(w) \, + \, .5 {\sigma}^2 \Bigg[{ (\frac{\partial w}{\partial Z})}^{-2} \,-Z \, \frac{\partial^2 w}{\partial Z^2}\,{\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-3} \, \Bigg] - \, .5 {\sigma}^2 \Bigg[ \, \Big[- \, \frac{\partial^3 w}{\partial Z^3} {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-3}+ 3 \, {(\frac{\partial^2 w}{\partial Z^2})}^2 \, {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-4} \Big] - \, \Big[ \, {(\frac{\partial^2 w}{\partial Z^2})}^2 {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-4}) \Big] \, \Bigg]$
$=\mu'(w) \, + \, .5 {\sigma}^2 \Bigg[{ (\frac{\partial w}{\partial Z})}^{-2} \,-Z \, \frac{\partial^2 w}{\partial Z^2}\,{\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-3} \, \Bigg] + \, .5 {\sigma}^2 \Bigg[ \, \Big[ \, \frac{\partial^3 w}{\partial Z^3} {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-3}- 2 \, {(\frac{\partial^2 w}{\partial Z^2})}^2 \, {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-4} \Big] \, \Bigg]$

So form of equations I have used in my program are
$f(w)=\mu(w) \, + \, .5 {\sigma}^2 Z \big[\frac{\partial w}{\partial Z}\big]^{-1} + \, .5 {\sigma}^2 \, \frac{\partial^2 w}{\partial Z^2} {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-2} \,$
and
$f'(w)=\mu'(w) \, + \, .5 {\sigma}^2 \Bigg[{ (\frac{\partial w}{\partial Z})}^{-2} \,-Z \, \frac{\partial^2 w}{\partial Z^2}\,{\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-3} \, \Bigg] + \, .5 {\sigma}^2 \Bigg[ \, \Big[ \, \frac{\partial^3 w}{\partial Z^3} {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-3}- 2 \, {(\frac{\partial^2 w}{\partial Z^2})}^2 \, {\big[ \, \frac{\partial w}{\partial Z} \, \big]}^{-4} \Big] \, \Bigg]$
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

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

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

Now I will explain the program and try to get the readers oriented with it.
In the program all leading coefficients of zeroth power are followed by a zero while coefficients of higher power are stored in an array.
In the expansion of w(t) given as
$w(t)=a_0 \,+ a_1 \, Z + a_2 \, Z^2 + a_3 \, Z^3 + a_4 \, Z^4 +\, ...$
In our program leading coefficient a0=$a_0$ is handled and input to functions separately while in the program array a has coefficients of  higher powers of power series of w(t) so
a(1)=$a_1$
a(2)=$a_2$
a(3)=$a_3$
and so on.
In the matlab program expansion of $\mu(w(t))$ is  given as
$\mu(w(t))=b_0 \,+ b_1 \, Z + b_2 \, Z^2 + b_3 \, Z^3 + b_4 \, Z^4 +\, ...$
In our program leading coefficient of drift b0=$b_0$ is handled and input to functions separately while in the program array b has coefficients of  higher powers of power series of $\mu(w(t))$ so
b(1)=$b_1$
b(2)=$b_2$
b(3)=$b_3$
and so on.

Similarly expansion of first derivative of drift w.r.t w(t) $\frac{\partial \mu(w)}{\partial w}$ is stored in array as b1 while its leading coefficient is denoted as b10 just like that of drift. Please note that first derivative of drift w.r.t w has its own power series and it is treated as a different function from drift.

Main function of the program that evolves the series solution is given as

function [c0,c] = EvolveSDESeriesSolution08New(a0,a,b0,b,b10,b1,sigma0,dt,ExpnOrder)
Here all the input arguments should be obvious now. These are coefficients of w, its drift and first derivtive of its drift followed by  by volatility coefficient, time step size and expansion order.
In output arguments c0 represents the leading coefficient of resulting power series at next time level and array c contains coefficients of first and higher powers of the power series.
In this function  d1=$\frac{\partial w}{\partial Z}$
d10 is leading coefficient of power series for first derivative and d1 is the array that stores coefficients of first and higher powers. d1(n) contains the coefficients of Z^n.
In this function  d2=$\frac{\partial^2 w}{\partial Z^2}$
d20 is leading coefficient of power series for second derivative derivative and d2 is the array that stores coefficients of first and higher powers. d2(n) contains the coefficients of Z^n.
Similarly d3 represents the third derivative of w with respect to Z and d30 is leading coefficient while d3 array contains coefficients of nth power of Z.
d1Sq array contains coefficients of power series of ${(\frac{\partial w}{\partial Z})}^2$ while d10Sq is its leading coefficient.
d1Cube array contains coefficients of power series of ${(\frac{\partial w}{\partial Z})}^3$ while d10Cube is its leading coefficient.
d1quad array contains coefficients of power series of ${(\frac{\partial w}{\partial Z})}^4$ while d10Quad is its leading coefficient.
d1SqInv array contains coefficients of power series of ${(\frac{\partial w}{\partial Z})}^{-2}$ while d10SqInv is its leading coefficient.
d1CubeInv array contains coefficients of power series of ${(\frac{\partial w}{\partial Z})}^{-3}$ while d10CubeInv is its leading coefficient.
d1quadInv array contains coefficients of power series of ${(\frac{\partial w}{\partial Z})}^{-4}$ while d10QuadInv is its leading coefficient.
d2Sq array contains coefficients of power series of ${(\frac{\partial^2 w}{\partial Z^2})}^2$ while d20Sq is its leading coefficient.
Z_d1Inv array contains coefficients of power series of $Z \,{(\frac{\partial w}{\partial Z})}^{-1}$ while Z_d1Inv0 is its leading coefficient.
There are a few other series with similarly obvious nomenclature.
T_Odt array contains coefficients of power series of $f(w)$ while T_Odt0 is its leading coefficient.
T_Odt2 array contains coefficients of power series of $f'(w)$ while T_Odt20 is its leading coefficient.
Term2 array contains coefficients of power series of $f(w) \, f'(w)$ while Term20 is its leading coefficient.

With previous explanation working of this function should be obvious.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

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

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

I still have to explain the functions related to drift and calculation of its power series.
There is a  function that converts derivatives of drift with respect to w into power series of Z given the coefficients of power series of w.
This function is
[b0,b] = CalculateDriftbCoeffs04(wMu0dt,dwMu0dtdw,a)
Its input is wMu0dt which is point value of drift at $w(Z_0)=a_0$ This is a scalar variable and not an array. while dwMu0dtdw is an array whose nth term contains nth derivative of drift with respect to w calculated at median. While a is an array whose nth term contains coefficient of nth power in power series of w.
In output, b0 is a scalar that is the leading coefficient of power series while b is an arry whose nth term contains coefficient of nth power of power series of drift.
CalculateDriftbCoeffs04 is a very general function and can be used to convert any function of w into power series of w when inputs that are derivatives of function with respect to w at median and an array that contains coefficients of power series of w.
In light of my previous posts working of this function should become obvious.

For the two term drift SDEs that we have been working with, the drift term in bessel coordinates is given as
$\mu(w) \,= \mu_1 [(1-\gamma)w]^{\frac{\beta_1-\gamma}{1-\gamma}} + \mu_2 [(1-\gamma)w]^{\frac{\beta_2-\gamma}{1-\gamma}} -.5 {\sigma}^2 \gamma [(1-\gamma)w]^{-1}$
The above equation is the bessel drift function $\mu(w)$ that we have used in our program. The expansion of this drift is used in calculation of $f(w)$ while expansion of its first derivative $\mu'(w)$ is used in the calculation of $f'(w)$
This is done in function
[wMu0dt,dwMu0dtdw,wMu1dt,dwMu1dtdw] = BesselDriftAndDerivatives08(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,ExpnOrder)
In the input side w0 is value of w at median while rest of the parameters are from the above drift function that I have latexed.
In the output  wMu0dt is the point value of drift at median and dwMu0dtdw is an array whose nth element contains nth derivative(wrt w) of the drift at median.
while  wMu1dt is the point value of first derivative of drift and dwMu1dtdw is an array that contains nth further derivative of the first derivative of drift calculated at median.
Incidentally  I want to mention that  wMu1dt=dwMu0dtdw(1) and  dwMu1dtdw(n)=dwMu0dtdw(n+1).
This should make working of this function obvious.
These point estimates and derivatives of drift are then fed to program
[b0,b] = CalculateDriftbCoeffs04(wMu0dt,dwMu0dtdw,a)
which converts them into power series of drift and its first derivative. This is done by calling the above program twice as
[b0,b] = CalculateDriftbCoeffs04(wMu0dta,dwMu0dtdwa,a);
[b10,b1] = CalculateDriftbCoeffs04(wMu1dt,dwMu1dtdw,a);

I am travelling tomorrow and will come back (in a day or two) with a new program that contains further  enhancements (that I did not have when I posted the last program yesterday) to get the series solution to work close to zero in a better way.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

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

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

Though I am sure everything about my series solution to constant CDF lines of density of SDEs would have become obvious to friends, some very basic things still might not be clear to some people who are still not oriented with the program so I will use a toy example to relate my addition of arrays that contain coefficients of power series of SDE to the addition of series coefficients.
Suppose we want to solve the following toy example with series solution
$w(t_1)=\mu(w(t_0)) + w(t_0)$

given that
$w(t_0)=a_0 \,+ a_1 \, Z + a_2 \, Z^2 + a_3 \, Z^3 + a_4 \, Z^4 +\, ...$
$\mu(w(t_0))=b_0 \,+ b_1 \, Z + b_2 \, Z^2 + b_3 \, Z^3 + b_4 \, Z^4 +\, ...$
$w(t_1)=c_0 \,+ c_1 \, Z + c_2 \, Z^2 + c_3 \, Z^3 + c_4 \, Z^4 +\, ...$

We have to find c's the coefficients associated with $w(t_1)$ given the known a's and b's that are coefficients of $w(t_0)$ and $\mu(w(t_0))$
Though we have found the series solution at Z=0, nowhere do we apply this zero value of Z(while calculating evolution equation to find c's), we simply do arithmetic on coefficients of the series.
$c_0=b_0 \, + \, a_0$
and coefficient of nth power of the resulting series
$c_n=b_n \, + \, a_n$

relating this to my matlab program where the convention is that leading coefficients have a variable name followed by a 0 and higher coefficients are stored in arrays with nth element of the array representing coefficient of nth power of Z.
So in the program we simply do series addition by adding leading coefficients resulting in leading coefficient of the new series. And then arrays containing higher coefficients are added to give the array whose nth element is coefficient of nth power of Z.
Throughout these calculations to get the coefficients of resulting series, we do arithmetic on series coefficients and never apply the value Z=0.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

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

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

Friends, I was in Kot Addu for past one week. I drove to Kot Addu last Thursday and drove back to Lahore yesterday. I am writing this post to describe the tactics of mind control agency to control me in the meantime. Last Thursday while going to Kot Addu, I took the route of Lahore-Multan motorway and got off the motorway at Khanewal exit. I knew that petty crooks of Pakistan army would try to drug the food on the way and I had previously thought that I would buy food and drinks from Lahore night before the journey but I recall that I was very tired at night before the travel and just decided to sleep. On motorway there are special areas where one could get petrol and buy food and drinks. I stopped at two of these recreation areas. At first place I just bought a diet coke and a sandwich from a small shop at the back. Both of these were good. That encouraged me to think that they might not be drugging the food on motorway and I later took another exit to a different recreation area some 120 kms ahead and bought water and red bull and another energy drink. Solid food I bought was good but all drinks including bottled water were badly drugged. Later I bought some diet cokes from a large store in Multan that lied on our path and found that all of them were drugged with mind control chemicals. I did not buy any food or drinks after that  during my drive to Kot Addu.
My first three days in in Kot Addu were very good. Everyday in the morning I would take the motorbike and drive on roads around the city with lush green fields and large trees. It was a very pleasant weather and I truly enjoyed driving motorcycle around the countryside and on small roads on banks of streams in the pleasant weather. Food and drinks in the Kot Addu city market were also very good. I was going to Kot Addu after two years and  there was far smaller mind control since it would be difficult to move all the mind control infrastructure quickly from Lahore. There were very few laser gadgets in my room since I had not lived in that room since more than eight years.
Probably since I was connecting more parts of my brain, I realized some urgency from mind control agencies to control me after first three days. They planted some device in motorcycle's seat (that was earlier not there for first few days) and it would constantly charge my anus as I drove the motorcycle and it was slightly painful ad disturbing. Then there was constant and irritating  gas in my room for last three days that was earlier not there. I also noticed that my family started to drug the food heavily in last three days. On Wednesday, My family ordered large amount of rice with chicken  from a professional cooking business to distribute among family friends and neighbors and all the rice was extremely heavily drugged. I had several cups of very strong tea after eating little rice and felt better. I went upstairs to my room and wrote a post on Wilmott and decided to go down to kitchen and have tea again. I had just a few sips of tea that I realized that I was fast losing my brain and I knew that somebody had drugged the tea in the meantime. That day I baked large amount of sweet potatoes in microwave oven and lived on them. So in the  last two days there was strong urgency among mind control agents to control me.
When I drove to Lahore yesterday from G.T.Road, I was able to buy some cheap energy drinks from rural surroundings of Kot Addu that were good. Later I bought two red bulls and a sting from  from chicha watni  and they were all badly drugged. I tried again to buy red bull and diet coke from Sahiwal and they were also drugged.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

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

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

Yesterday, I realized that mind control agents  had started drugging the city in a big way to target me. But Lahore is a very big city and it is very hard to drug food and water in a few days and first few days are relatively easy if the victim is smart enough. I had rice biryani from a place where I had never gone and they would least suspect I could go (As I said it is relatively easy at start). Then I bought Nestle bottled water from inner streets of an area a few kilometers away from where I live and I was surprised to find that water had already been drugged. I knew that Nestle water had been very good for past few months and the water had absolutely not been drugged at the manufacturing source and it must have been drugged in past two days. I knew that I would have to be very careful about water. I threw bottled water and I was lucky to find a water-cooler on unpopulated side of the canal (It was outside a marble tiles business) and I knew that water there would be good(at least in early days when drugging of water has just started) and was able to fill the bottles with water. I continued to work in the evening yesterday. I had some mild suspicion that they were changing my programs but I could not be certain therefore not emphasizing it.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

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

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

My past two days were very bad. Day before yesterday, I became a bit lax and had chicken biryani from a Biryani shop. It was drugged. I could not tell very well until I had taken all of the biryani. And later again yesterday I was hit when I had breakfast. I was returning from my morning walk and it would be around 7:00 am. I did not intend to eat anything but saw chicken kachoris and had one. It was perfectly good and that prompted me to think that other food at that place might be good as well. I had some chickpeas and they were drugged so I remained very down and sluggish all day yesterday. They also used gases on me almost all night in my room and I was feeling sick when I woke up. Really please ask them to stop using gases.
This morning I did not go for a walk ( I try to walk every other day) and I bought bottled Nestle water from a far off place in Wapda town but the water was drugged. I threw water in the bottles and filled them with tap water at a society in the maze of housing societies next to Valencia. This water had a slightly off taste and it also seemed very lightly drugged. I did not fill all bottles with water and tried to fill them from a different society(probably UET employers housing society) and the water was far more drugged with mind control chemicals there. I decided to use more drugged water for washing my face and body and thought that I would drink very lightly drugged water as such. I also realized that Pakistan army crooks had again drugged underwater reservoirs with mind control drugs in large parts of Lahore city.
I continue to take large amount of tea every day to lift myself and be able to do my research. If it were not for tea, I would be very down all the time also because of antipsychotics. I cannot take coffee since every coffee brand I tried in past one year was thoroughly drugged. I do not know how they do it but coffee beans are also drugged with mind control chemicals. So taking coffee to stimulate the mind is not a possibility. I therefore started taking tea as a (obviously less effective) substitute for coffee but now I am afraid they want to ask all tea manufacturers in Pakistan to drug the tea with mind control chemicals. And of course, this is all arranged by Pakistan army since many of their generals take huge bribes from mind control agencies in return for victimizing their own compatriots.
Finally I hope to post new version of the program with slightly better simulation of densities of SDEs close to zero. I will try to post it later today or tomorrow. I would have already done it if it were not for highly increased mind control activity in past few days after I came back to Lahore since I was hit by drugged food several times and remained very sluggish for past two days.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

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

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

Friends, very sorry for delay in the new program that works somewhat better close to zero. I am copying the condition I have imposed on coefficients from my matlab program.

for nn=1:ExpnOrder-1
if(abs(a(nn+1))>abs(a(nn))*2^(-nn+1))
a(nn+1)=sign(a(nn+1))*abs(a(nn))*2^(-nn+1);
end
end

here a(nn) is nn'th coefficient of power series that multiplies  Z^nn, Please feel free to play around with this growth restriction on absolute value of the coefficients.

I am posting the new program here.
function [] = FPESeriesCoeffsSolutionWilmottOdt2Downloaded()

%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.

dt=.125/16*4;   % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=128*5/4;%128*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/1;
TtM=T/dtM;

dNn=.1/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=46*2;  % No of normal density subdivisions
%NnMidl=18;%One half density Subdivision left from mid of normal density(low)
%NnMidh=19;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;
%NnT=46;
%NnD=(NnT-Nn)/2;
%NnD=5;

x0=.250;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.65;%50;   % volatility power.
kappa=1.0;%.950;   %mean reversion parameter.
theta=.05;%mean reversion target
sigma0=.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;
w0=x0^(1-gamma)/(1-gamma);

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

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

sigma11(1:OrderA+1)=0;
mu11(1:OrderA+1)=0;
mu22(1:OrderA+1)=0;
sigma22(1:OrderA+1)=0;
% index 1 correponds to zero level since matlab indexing starts at one.
sigma11(1)=1;
mu11(1)=1;
mu22(1)=1;
sigma22(1)=1;

for k=1:(OrderA+1)
if sigma0~=0
sigma11(k)=sigma0^(k-1);
end
if mu1 ~= 0
mu11(k)=mu1^(k-1);
end
if mu2 ~= 0
mu22(k)=mu2^(k-1);
end
if sigma0~=0
sigma22(k)=sigma0^(2*(k-1));
end
end
%Ft(1:TtM+1,1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0; %General time powers on hermite polynomials
Fp(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers on coefficients of hermite polynomials.
Fp1(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;%General x powers for bessel transformed coordinates.

%YCoeff0 and YCoeff are coefficents for original coordinates monte carlo.
%YqCoeff0 and YqCoeff are bessel/lamperti version monte carlo.

YCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
YqCoeff0(1:(OrderA+1),1:(OrderA+1),1:(OrderA+1),1:(OrderA+1))=0;
%Pre-compute the time and power exponent values in small multi-dimensional arrays
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma); %expand y^alpha where alpha=1;
YqCoeff = ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);%expand y^alpha1 where alpha1=(1-gamma)
YqCoeff=YqCoeff/(1-gamma); %Transformed coordinates coefficients have to be
%further divided by (1-gamma)

for k = 0 : (OrderA)
for m = 0:k
l4 = k - m + 1;
for n = 0 : m
l3 = m - n + 1;
for j = 0:n
l2 = n - j + 1;
l1 = j + 1;
%Ft(l1,l2,l3,l4) = dtM^((l1-1) + (l2-1) + (l3-1) + .5* (l4-1));
Fp(l1,l2,l3,l4) = (alpha + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
- (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));
Fp1(l1,l2,l3,l4) = (alpha1 + (l1-1) * beta1 + (l2-1) * beta2 + (l3-1) * 2* gamma + (l4-1) * gamma ...
- (l1-1) - (l2-1) - 2* (l3-1) - (l4-1));

YCoeff0(l1,l2,l3,l4) =YCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
YqCoeff0(l1,l2,l3,l4) =YqCoeff(l1,l2,l3,l4).*mu11(l1).*mu22(l2).*sigma22(l3).*sigma11(l4);
end
end
end
end

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

ExpnOrder=6;

wnStart=1;
tic

for tt=1:Tt
if(tt==1)
[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivatives(w0,YqCoeff0,Fp1,gamma,dt,ExpnOrder);
c0=wMu0dt+w0;
c(1)=sigma0*sqrt(dt)
c(2:10)=0.0;
w0=c0;
end

%The program intended to find solution of Z-series to 10th power but
%higher order coefficients would make some SDEs unstable so I just
%turned the coefficients c(7:10) and b(7:10) into zero as soon as they
%are calculated and just go with solution to 6th power of Z.

if(tt>1)

%below wMu0dta is drift and dwMu0dtdwa are its derivatives.
%wMu1dt is its first derivative and dwMu1dtdw are derivatives of first
%derivative.

%[wMu0dta,dwMu0dtdwa,wMu1dt,dwMu1dtdw] = BesselDriftAndDerivatives08(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,10);
[wMu0dta,dwMu0dtdwa,wMu1dt,dwMu1dtdw] = BesselDriftAndDerivatives08Bounded(w0,mu1,mu2,beta1,beta2,sigma0,gamma,dt,10);

%The program CalculateDriftbCoeffs04 takes value of drift and its
%derivatives wrt w and also takes Z-series of w and converts them into
%Z-series of drift. It is used to calculate Z-series of drift and Then
%to calculate Z_series of first derivative of drift wrt w.
[b0,b] = CalculateDriftbCoeffs08(wMu0dta,dwMu0dtdwa,a,gamma);
[b10,b1] = CalculateDriftbCoeffs08(wMu1dt,dwMu1dtdw,a,gamma);

b(ExpnOrder+1:10)=0.0;

b1(ExpnOrder+1:10)=0.0;
%The function below evolves the previous solution in coefficients of Z
%into next time step solution in coefficietns of Z-series.
[c0,c] = EvolveSDESeriesSolution08New(a0,a,b0,b,b10,b1,sigma0,dt,ExpnOrder,gamma);
c(ExpnOrder+1:10)=0.0;

w0=c0;

end

a0=c0;

a(1:ExpnOrder)=c(1:ExpnOrder);
a(ExpnOrder+1:10)=0;

Bound=abs(a(1)/a0)*gamma/(1-gamma);

for nn=1:ExpnOrder-1

%if(abs(a(nn+1))>abs(a(nn))*Bound)
%    a(nn+1)=sign(a(nn+1))*abs(a(nn))*Bound;
%end
if(abs(a(nn+1))>abs(a(nn))*2^(-nn+1))
a(nn+1)=sign(a(nn+1))*abs(a(nn))*2^(-nn+1);
end

end

end
c0=a0;
c=a;
wnStart=1;

w(wnStart:Nn)=c0;
for nn=1:ExpnOrder
w(wnStart:Nn)=w(wnStart:Nn)+c(nn)*Z(wnStart:Nn).^nn;
end
c0
c

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

Dfyy(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1

Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));

%Change of variable derivative for densities
end
Dfyy(Nn)=Dfyy(Nn-1);
Dfyy(1)=Dfyy(2);

pyy(1:Nn)=0;
for nn = wnStart:Nn

pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
end

toc

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

theta1=1;
rng(29079137, 'twister')
paths=50000;
YY(1:paths)=x0;  %Original process monte carlo.
Random1(1:paths)=0;
for tt=1:TtM

Random1=randn(size(Random1));
HermiteP1(1,1:paths)=1;
HermiteP1(2,1:paths)=Random1(1:paths);
HermiteP1(3,1:paths)=Random1(1:paths).^2-1;
HermiteP1(4,1:paths)=Random1(1:paths).^3-3*Random1(1:paths);
HermiteP1(5,1:paths)=Random1(1:paths).^4-6*Random1(1:paths).^2+3;

YY(1:paths)=YY(1:paths) + ...
(YCoeff0(1,1,2,1).*YY(1:paths).^Fp(1,1,2,1)+ ...
YCoeff0(1,2,1,1).*YY(1:paths).^Fp(1,2,1,1)+ ...
YCoeff0(2,1,1,1).*YY(1:paths).^Fp(2,1,1,1))*dtM + ...
(YCoeff0(1,1,3,1).*YY(1:paths).^Fp(1,1,3,1)+ ...
YCoeff0(1,2,2,1).*YY(1:paths).^Fp(1,2,2,1)+ ...
YCoeff0(2,1,2,1).*YY(1:paths).^Fp(2,1,2,1)+ ...
YCoeff0(1,3,1,1).*YY(1:paths).^Fp(1,3,1,1)+ ...
YCoeff0(2,2,1,1).*YY(1:paths).^Fp(2,2,1,1)+ ...
YCoeff0(3,1,1,1).*YY(1:paths).^Fp(3,1,1,1))*dtM^2 + ...
((YCoeff0(1,1,1,2).*YY(1:paths).^Fp(1,1,1,2).*sqrt(dtM))+ ...
(YCoeff0(1,1,2,2).*YY(1:paths).^Fp(1,1,2,2)+ ...
YCoeff0(1,2,1,2).*YY(1:paths).^Fp(1,2,1,2)+ ...
YCoeff0(2,1,1,2).*YY(1:paths).^Fp(2,1,1,2)).*dtM^1.5) .*HermiteP1(2,1:paths) + ...
((YCoeff0(1,1,1,3).*YY(1:paths).^Fp(1,1,1,3) *dtM) + ...
(YCoeff0(1,1,2,3).*YY(1:paths).^Fp(1,1,2,3)+ ...
YCoeff0(1,2,1,3).*YY(1:paths).^Fp(1,2,1,3)+ ...
YCoeff0(2,1,1,3).*YY(1:paths).^Fp(2,1,1,3)).*dtM^2).*HermiteP1(3,1:paths) + ...
((YCoeff0(1,1,1,4).*YY(1:paths).^Fp(1,1,1,4)*dtM^1.5 )).*HermiteP1(4,1:paths) + ...
(YCoeff0(1,1,1,5).*YY(1:paths).^Fp(1,1,1,5)*dtM^2.0).*HermiteP1(5,1:paths);

end

YY(YY<0)=0;
disp('Original process average from monte carlo');
MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
disp('Original process average from our simulation');
ItoHermiteMean=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates

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

MaxCutOff=30;
NoOfBins=round(1*900*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );
plot(yy(wnStart+1:Nn-1),pyy(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');

title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMean,TrueMean));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));

legend({'Ito-Hermite Density','Monte Carlo Density'},'Location','northeast')

str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');

end

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

NoOfTerms=3;

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

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

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

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

for mm=1:NoOfTerms

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

end

Bound=1/abs(wMu0dt/dwMu0dtdw(1))*gamma^1;

for nn=1:ExpnOrder

if(abs(dwMu0dtdw(nn+1))>abs(dwMu0dtdw(nn)*Bound))
dwMu0dtdw(nn+1)=sign(dwMu0dtdw(nn+1))*abs(dwMu0dtdw(nn))*Bound;
end
end

wMu1dt=dwMu0dtdw(1);
dwMu1dtdw(1:ExpnOrder)=dwMu0dtdw(2:ExpnOrder+1);

%wMu0dt
%dwMu0dtdw

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

end
Rest of the functions are old. When you run this program, you will get the output as

ItoHermiteMean =
0.052494659280937
true Mean only applicable to standard SV mean reverting type models otherwise disregard
TrueMean =
0.051347589399817
yy0 =
0.033647208005418
Original process average from monte carlo
MCMean =
0.051499886999017 - 0.000000442861801i
Original process average from our simulation
ItoHermiteMean =
0.052494659280937
true Mean only applicble to standard SV mean reverting type models otherwise disregard
TrueMean =
0.051347589399817
IndexMax =
1677
red line is density of SDE from Ito-Hermite method, green is monte carlo.

and the following graph.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

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

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

Friends, very sorry for delay in the new program that works somewhat better close to zero. I am copying the condition I have imposed on coefficients from my matlab program.

for nn=1:ExpnOrder-1
if(abs(a(nn+1))>abs(a(nn))*2^(-nn+1))
a(nn+1)=sign(a(nn+1))*abs(a(nn))*2^(-nn+1);
end
end

here a(nn) is nn'th coefficient of power series that multiplies  Z^nn, Please feel free to play around with this growth restriction on absolute value of the coefficients.
Please also try a related growth condition given as
if(abs(a(nn+1))>abs(a(nn))*2^(-nn/2))
a(nn+1)=sign(a(nn+1))*abs(a(nn))*2^(-nn/2);
end
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

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

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

Friends, Mind control agencies have been incessantly trying to control me since my last few days at Kot Addu more than a week ago. Whenever I am creative, it starts to ring bells in minds of people who want to fail me and they become active again to persecute me. And mind control officials who want to enter Jewish Godfather's paradise start a new round of extreme persecution.
Today I was hit again when I had a sting (which is a cheap energy drink) from Model Town and my condition remains deteriorating after that. They forced feeling of anxiety on me after I took the drink and highly increased my back ache. I have constant and increasing back ache for past 3-4 days which is caused since they pull neurotransmitters out of my brain.
Everyday since past several months, they cause dryness in my mouth when I sleep and something settles on upper part of my tongue and my upper jaw. When I wash my mouth with water, very large amount of foam comes out. I wash my mouth 2-4 times every night when I wake up during sleep.
Every time I sit in my car, they charge my head and face and I remain lightly tense after that. However when I wash my head and face after getting out of my car or after I return home, all of a sudden I feel fresh and lively again.
They also continue to release similar gas from time to time in my room to charge my face, head, nose and inside of my mouth. When I breath in gas that is when my nose and inside of my mouth gets charged.
I want to request friends to please help end my persecution since it never ends and continues to increase every day now.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

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

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

Friends, I had a sleep after writing the previous post and felt a bit better and was able to do some work.
Tomorrow I intend to write a post with the equations for solution of constant CDF lines of dz-integrals and dt-integrals of functions of SDE through their fokker-planck equation overlaid on the fokker-planck  equation of original SDE. Both fokker-planck equations are related by a common Z and we can solve the dependent fokker-planck equation of dz-integrals and dt-integrals  after solving the FP equation of independent SDE. I am working out the details and will try to share it with a nice post in the form of equations tomorrow. Though a worked out program might be about a week away.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

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

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

Friends we will try to find out how to find the series expansion and density of dt-integrals.

Suppose we are given an SDE of the form in bessel coordinates

$dw(t) \, = \, \mu(w) \, dt + \sigma_0 \,dz(t)$

please note that bessel transformation is given as $w(t)=\frac{{X(t)}^{(1-\gamma)}}{(1-\gamma)}$

So if our path integral in original coordinates is given by integral $\int_0^T \, X(t)^{\alpha} \, dt$

this integral in bessel coordinates would be given as  $\int_0^T \, \big[(1-\gamma) \, w(t) \big]^{\frac{\alpha}{(1-\gamma)}} \, dt$

We suppose $f_w(t)=\big[(1-\gamma) \, w(t) \big]^{\frac{\alpha}{(1-\gamma)}}$

and we want to find the integral  $\int_0^T \, f_w(t) \, dt$

We know that  $d\big[t \, f_w(t)]=\, t \, df_w(t) \, + \, f_w(t) \, dt$

giving us the formula  $\, f_w(t) \, dt \, =\, d\big[t \, f_w(t)]\,- t \, df_w(t) \,$

We know that $t \, df_w(t) \,=\, t\, \mu_f(w) \, dt \,+\, t\, \sigma_f(w) \, dz(t) \,$

We introduce another function $dg_w(t)\, =\, t \, df_w(t) \,=\, t\, \mu_f(w) \, dt \,+\, t\, \sigma_f(w) \, dz(t) \,$

Since  $\, f_w(t) \, dt \, =\, d\big[t \, f_w(t)]\,- t \, df_w(t) \,$

It can be written as  $\, f_w(t) \, dt \, =\, d\big[t \, f_w(t)]\,- \, dg_w(t)\, \,$

Taking integral on the above formula, we get  $\int_0^T \, f_w(t) \, dt \, =\, \int_0^T\, d\big[t \, f_w(t)]\,-\, \int_0^T \, dg_w(t)\, \,$
$=\int_0^T \, f_w(t) \, dt \, =\, T \, f_w(T)\,-\, \int_0^T \, dg_w(t)\, \,$

The value of term $\, T \, f_w(T)\,$ can be found directly from the transformation of solution of the original fokker-planck equation. It will be a power series in terms of Z in our case. We can write the fokker-planck equation of $g_w(t)\,$ and once we solve this fokker-planck equation overlaid on the Fokker-planck equation of $w(t)$, we would be able to find the series solution in terms of Z given by the equation

$\int_0^T \, f_w(t) \, dt \, =\, T \, f_w(T)\,-\, \, g_w(T)\, -\, g_w(0)$

The above equation will be an expression in terms of a series in Z and for that we would have to solve the dependent fokker-planck equation of  $g_w(t)$.  $g_w(0)$ in above equation would usually be zero.

The above equation could also be written for a small step as

$\int_{t_1}^{t_2} \, f_w(t) \, dt \, =\, t_2 \, f_w(t_2)\,- \, t_1 \, f_w(t_1)\, + \, g_w(t_2)\, -\, g_w(t_1)$

In next post, I will write equations that I have used to solve the fokker-planck equation of  $g_w(t)$ overlaid on fokker-planck equation of $w(t)$
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal

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

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

In previous post, we introduced a function by its differential equation as

$dg_w(t)\, =\, t \, df_w(t) \,=\, t\, \mu_f(w) \, dt \,+\, t\, \sigma_f(w) \, dz(t) \,$

it is $g_w(t)$ whose overlaid fokker-planck we want to find out. In some equations that follow, I have simply written $g$ instead of $g_w$.

We write this fokker-planck equation as

$\frac{\partial P(g_w,t)}{\partial t} = -\frac{\partial \big[t \, \mu_f(w) \, P(g_w,t) \big]}{\partial g}\, +.5 \, \frac{\partial^2 \big[t^2 \, {\sigma_f(w)}^2 \, P(g_w,t) \big]}{\partial g^2}$

To find out the constant CDF lines of the above equation, We follow the steps given in the post https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1200#p867718

After appropriate manipulations, we get the solution of above equation as

$\frac{\partial g_b}{\partial t} =t \mu_f(w_b) - .5 \, t^2 \, \frac{\partial \big[{\sigma_f(w_b)}^2 \big]}{\partial g}+ .5\, \, t^2 \,{\sigma_f(w_b)}^2 \, Z_b \, \frac{\partial Z}{\partial g} + \, .5\, t^2 \, {\sigma_f(w_b)}^2 \, {(\frac{\partial Z}{\partial g})}^2 \frac{\partial^2 g}{\partial Z^2}$

In the above equation suffix b means that this equation is valid for moving along a particular CDF line or equivalently by a particular value of Z indicated by $Z_b$

When finding out the series solution of $g_w(t)$  we would have to copy the power series (in terms of Z) of  $\mu_f(w)$ and $\sigma_f(w)$ from the fokker-planck of original SDE w(t) at every step in time and apply it to dependent fokker-planck of $g_w(t)$ . This is done at median of both fokker-planck equations.

Another problem is that we need to find the derivative  $\, \frac{\partial \big[{\sigma_f(w_b)}^2 \big]}{\partial g}$ which is not straightforward since  ${\sigma_f(w_b)}$  is not in terms of g. We would have to calculate this derivative with chain rule as

$\, \frac{\partial \big[{\sigma_f(w_b)}^2 \big]}{\partial g} \,= \, \frac{\partial \big[{\sigma_f(w_b)}^2 \big]}{\partial w} \, \frac{\partial w}{\partial Z} \, \frac{\partial Z}{\partial g}$

please note the in above equation all of the three terms  $\frac{\partial \big[{\sigma_f(w_b)}^2 \big]}{\partial w}$ , $\frac{\partial w}{\partial Z}$ ,  $\frac{\partial Z}{\partial g}$  would be different power series in Z evaluated at median Z=0. The first two terms would be calculated at median of base fokker-planck equation and last third term would be calculated at median of dependent fokker-planck equation.

So in order to solve the overlaid fokker-planck equations in a power series,  all we need to do is copy the series of w-terms (that are expressions in w) from their power series at median from the original fokker-planck equation and plug them in the dependent fokker-planck equation at median. We will also use chain rule to convert the derivatives of these w-terms with respect to $g_w$ and all involved terms are expressed as power-series in Z and are calculated at median of the both fokker-planck equations.

I have not yet calculated 2nd order expansion of the fokker-planck in terms of a Z power series. In next few days I will be working on this equation and code it in matlab and post the matlab program here.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal