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

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

In my previous post, I suggested to expand the SDE variable $X(t)$ as a function of hermite polynomials of the brownain motion but then I realized that drift and other functions of X(t) are basically non-linear functions and they non-linearly affect various hermite polynomials in the hermite polynomial expansion of X(t). And then I thought that if we can expand the drift and other functions of X(t) in terms of hermite polynomials, we would probably not even need to solve the fokker planck equation at all.  Here are some of my more developed thoughts that I want to share with friends but friends need to think analytically on their own since I have not tested my ideas with experiments and there might be mistakes in untested ideas. But I hope that there is enough food for thought to share the ideas with friends at this stage.

The main theme here is to find an analytical expansion of X(t) in terms of hermite polynomials and then find hermite expressions of $\mu(X)$ and $\sigma(X)$ and then analytically evolve the coefficients of hermite expansion of $X(t)$ over each time step
First some basics, we are dealing with hermite polynomials in standard normal $Z$ but most of it is directly applicable to hermite polynomials of broawnian motion $B(t)$ as well.

Differentiation of hermite polynomials      $\frac{d H_{n+1}(Z)}{dZ}=(n+1) H_{n}(Z)$
repeated integrations of integral over normal random variable as
$\int_0^Z dZ_1=Z=H_1(Z)$
$\int_0^Z \int_0^{Z_1} dZ_2 dZ_1=\frac{Z^2-1}{2}=\frac{1}{2} H_2(Z)$
$\int_0^Z \int_0^{Z_1} \int_0^{Z_2} dZ_3 dZ_2 dZ_1=\frac{1}{6} (Z^3-3Z)=\frac{1}{6}H_3(Z)$
and so on.
First of all we will deal with the case of Lamperti form of SDEs since it is lot easier to explain with very little latexing and lesser equations.
Suppose SDE under consideration is $X(t)=\mu(X) dt + \sigma dB(t)$
When we start the evolution of a density from delta source, the very initial small time step form is essentially a normal. We can denote the first step evolution as a Hermite form as given by equation  $X(Z)=X_0+\mu(X_0) \Delta t + \sigma \sqrt{\Delta t} Z=X_0+\mu(X_0) \Delta t + \sigma \sqrt{\Delta t} H_1(Z)$
We can write it as $Z=a_0+a_1 H_1(Z)$ where $a_0=X_0+\mu(X_0) \Delta t$ and $a_1=\sigma \sqrt{\Delta t}$

Now when we take the second step in the evolution of SDE, it is more involved since we are not starting from a delta origin anymore. Obviously we have to add the volatility in a squared fashion to second hermite polynomial but we also have to represent the drift as a function of hermite polynomials and add it linearly to the existing one step evolved hermite representation of X(t).
Next we have to learn how to expand functions of X(t) in terms of hermite polynomials once we are given a hermite representation of X(t) itself. We take $\mu(X)$ as a test function and expand it in terms of hermite polynomials. This has to follow iterated integrals form of Taylor series as

$\mu(X(Z))=\mu(X(Z_0)) + \int_0^Z \frac{d \mu(X)}{dX} \frac{dX}{dZ} dZ$
Above we have expanded the $\mu(X)$ to one term and evaluating the integral at $X_0=X(Z_0)$ and expanding one step further, we get
$\mu(X(Z))=\mu(X(Z_0)) + \frac{d \mu(X_0)}{dX} \frac{dX_0}{dZ} \int_0^Z dZ_1$
$+ \int_0^Z \int_0^{Z_1} \big[ \frac{d^2 \mu(X)}{dX^2} {(\frac{dX}{dZ})}^2 +\frac{d \mu(X)}{dX} \frac{dX^2}{dZ^2} \big] dZ_2 dZ_1$

When expanded to thirs hermite polynomial, drift can be written as
$\mu(X(Z))=\mu(X(Z_0)) + \frac{d \mu(X_0)}{dX} \frac{dX_0}{dZ} \int_0^Z dZ_1$
$+\big[ \frac{d^2 \mu(X_0)}{dX^2} {(\frac{dX}{dZ})}^2 +\frac{d \mu(X_0)}{dX} \frac{dX^2}{dZ^2} \big] \int_0^Z \int_0^{Z_1} dZ_2 dZ_1$
$+\big[ \frac{d^3 \mu(X_0)}{dX^3} {(\frac{dX}{dZ})}^3+3 \frac{d^2 \mu(X_0)}{dX^2} {(\frac{dX}{dZ})}^2 \frac{d^2 X}{dZ^2}+\frac{d \mu(X_0)}{dX} \frac{dX^3}{dZ^3} \big] \int_0^Z \int_0^{Z_1} \int_0^{Z_2} dZ_3 dZ_2 dZ_1 + .....$

which will be equivalent when written in form of hermite polynomials as
$\mu(X(Z))=\mu(X(Z_0)) + \frac{d \mu(X_0)}{dX} \frac{dX_0}{dZ} H_1(Z)$
$+\big[ \frac{d^2 \mu(X_0)}{dX^2} {(\frac{dX}{dZ})}^2 +\frac{d \mu(X_0)}{dX} \frac{dX^2}{dZ^2} \big] \frac{H_2(Z)}{2}$
$+\big[ \frac{d^3 \mu(X_0)}{dX^3} {(\frac{dX}{dZ})}^3+3 \frac{d^2 \mu(X_0)}{dX^2} {(\frac{dX}{dZ})}^2 \frac{d^2 X}{dZ^2}+\frac{d \mu(X_0)}{dX} \frac{dX^3}{dZ^3} \big] \frac{H_3(Z)}{6} + O (H_4(Z))$

All derivatives in the above are taken at $X=X(Z_0)$. We have to take enough terms in hermite expansion of functions of X(t) so as to be able to faithfully follow the true function over the range of Z=-4 to +4. In the above equation we can see that even though the hermite representation of X(t) can take just one derivative, the hermite representation of functions of X(t) could be far more complex and could take significantly many more hermite polynomials in expansion.
Continuing after one step evolution of Lamperti SDE given as $X(Z)=a_0+a_1 H_1(Z)$, we get a representation for drift as

$\mu(X(Z))=\mu(X(Z_0)) + \frac{d \mu(X_0)}{dX} a_1 H_1(Z)$
$+\frac{d^2 \mu(X_0)}{dX^2} {a_1}^2 \frac{H_2(Z)}{2}+\frac{d^3 \mu(X_0)}{dX^3} {a_1}^3 \frac{H_3(Z)}{6}$

The above expression for drift could be linearly added by each hermite polynomial to the first one step evolution of X(t) but first, in case of Lamperti, we would have to add the one step volatility in a squared fashion to first order hermite polynomial and then linearly add the drift and then move on step by step for new updated representations of drift and adding them to the previous stage hermite representation of X(t) and also adding volatility in a squared fashion to first hermite polynomial. Particularly after two steps from the delta origin, the hermite representation of X(t) would be

$X(Z)=X_0+ \mu(X_0) \Delta t + \mu(X(Z_0)) \Delta t + [ \sigma \sqrt{\Delta t +\Delta t} + \frac{d \mu(X_0)}{dX} \Delta t \sigma \sqrt{\Delta t}] H_1(Z)$
$+\frac{d^2 \mu(X_0)}{dX^2} {( \sigma \sqrt{\Delta t})}^2 \Delta t \frac{H_2(Z)}{2}+\frac{d^3 \mu(X_0)}{dX^3} {( \sigma \sqrt{\Delta t})}^3 \Delta t \frac{H_3(Z)}{6}$
where I have used $a_1=\sigma \sqrt{\Delta t}$
The third step would be even more complicated since third step would get a contribution in expansion of $\mu(X)$ in terms of hermite polynomials from large number of X(t) expansion hermite terms from step two.

Now we come towards a general evolution algorithm using hermite polynomials. First we take a one step evolution equation in terms of $\mu(X)$, $\sigma(X)$, and their derivatives with respect to X(t) and hermite polynomials of brownian motion$B(t)$. Secondly we find hermite representation of all functions of drift and volatility and their derivatives in terms of X(t) in the above terms in terms of hermites of normal random variable around $Z=0$. And then we add all volatility terms to existing previous step hermite representation of X(t) in a squared fashion and then we add all drift terms in hermite form to the hermite form of X(t) in a linear fashion.

Lamperti form is very attractive since it would take far smaller no of hermite polynomials in its expansion as its density remains reasonably similar to normal random variable density and we do not need to expand the volatility in this case.

I will soon be distributing the worked out code with this evolution scheme and then we move to full fledged solution of Fokker-Planck equation. Amin
Topic Author
Posts: 2152
Joined: July 14th, 2002, 3:00 am

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

Friends, first please apologize me for writing this post without numerically checking the results. I had some basic ideas today that I want to share with friends hence the reason for writing this post.
When we are doing an evolution or simulation scheme for the SDEs, we take an estimate of drift over a discrete period of time and also an estimate of volatility over the same discrete period of time. However many times the dynamics of drift are non-linear or mean reverting so that there is an instantaneous cancellation or addition to volatility effect and this instantaneous cancellation/addition is spread over the discrete time period chosen for simulation. Though we have considered the " discrete effects" because of drift and volatility, we also need to account for the instantaneous cancellation/addition and its cumulative effect on the "discrete effects" of volatility and drift. If we can calculate the first order instantaneous derivatives of the drift and volatility "discrete values" in the same currency, all we need is a correction to "discrete values simulation" equal to integrated values of the difference of first order instantaneous derivatives of the drift and volatility when calculated in the same currency.
Ok, for that purpose, I calculate a representation of drift and its first derivative in terms of derivatives with respect to brownian motion. In what follows below, B(t) is brownian motion and Z is standard normal.

We write $\frac{\partial \big [\mu(X) dt \big ]}{\partial B(t)} = \frac{\partial \big [\sqrt{(\mu(X) dt)} \sqrt{(\mu(X) dt)} \big ]}{\partial B(t)} = 2 \sqrt{(\mu(X) dt)} \frac{\partial \big [\sqrt{(\mu(X) dt)} \big ]}{\partial B(t)}$
$= 2 \sqrt{(\mu(X) dt)} \frac{\partial \big [\sqrt{(\mu(X) )} \big ]}{\partial Z}$
The above is the first derivative of drift in terms of brownian motion. Since drift starts from zero, this derivative is not very interesting. When multiplied with brownian motion (as an orthogonal basis) , it will retrieve the drift itself. Here, We are not particularly interested in drift but its first derivative. For that purpose, we take a second derivative with respect to brownian motion.
Taking a second derivative with respect to brownian motion ,we get

$\frac{\partial^2 \big [\mu(X) dt \big ]}{\partial {B(t)}^2} =2 {\Big[ \frac{\partial \big [\sqrt{(\mu(X) dt)} \big ]}{\partial B(t)}\Big] }^2+ 2 \sqrt{(\mu(X) dt)} \frac{\partial^2 \big [\sqrt{(\mu(X) dt)} \big ]}{\partial {B(t)}^2}$
which after a few manipulations and adjustment of $\sqrt{dt}$ yields
$\frac{\partial^2 \big [\mu(X) dt \big ]}{\partial {B(t)}^2} =2 {\Big[ \frac{\partial \big [\sqrt{(\mu(X))} \big ]}{\partial Z}\Big] }^2+ 2 \sqrt{(\mu(X))} \frac{\partial^2 \big [\sqrt{(\mu(X))} \big ]}{\partial {Z}^2}$
After expanding in terms of partial derivatives and cancelling we get the following expression for the derivative of first discrete effect of drift as
$\frac{\partial^2 \big [\mu(X) dt \big ]}{\partial {B(t)}^2} =\frac{\partial \mu(X)}{\partial X} \frac{\partial^2 X}{\partial N^2}+\frac{\partial^2 \mu(X)}{\partial X^2} {(\frac{\partial X}{\partial N})}^2$
Since this derivative is with respect to brownain motion, we can find a representation of its total effect in terms of hermite orthogonal basis of brownian motion as
$\Big[ \frac{\partial \mu(X)}{\partial X} \frac{\partial^2 X}{\partial N^2}+\frac{\partial^2 \mu(X)}{\partial X^2} {(\frac{\partial X}{\partial N})}^2 \Big] (\frac{{B(t)}^2-t}{2})$

When we are in a lamperti framework, the contribution from higher derivatives of brownian motion in volatility is almost zero and in that case I believe that it is the above term that has to be added/subtracted with second hermite basis in our discrete simulation framework to accommodate the integrated effect of integrated instantaneous first derivatives as I described earlier. I believe it is the above derivative effect that is being added from the mean-reverting/non-linear drift to the volatility effect instantaneously.

But I might be totally wrong as I have not numerically checked anything but I hope the ideas would be helpful.

I will be posting a full-fledged matlab program with the above effects included and describe here how it works. Amin
Topic Author
Posts: 2152
Joined: July 14th, 2002, 3:00 am

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

An informal update for time being
The right evolution equation along the CDF lines incorporating the effect of first and second mode (1st and 2nd equation post 815) is given as

$X(t+\Delta t)=X(t)+(B(t + \Delta t) -B(t)) {(\frac{\partial X}{\partial B})} + \mu(X) {\Delta t} -2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} {\Delta t}$
$-.5 {\sigma(X)}^2 (+3 {\frac{\partial B}{\partial X}}^{-1} \frac{\partial ^2 B}{\partial X^2}) {\Delta t} - ({(\frac{\partial X}{\partial B})} - {\sigma(X)}^2{(\frac{\partial B}{\partial X})}) (B(t + \Delta t) -B(t))$ (Equation(6) in post 815)

which can be written after simplifying it as

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

here in the above equation $(B(t + \Delta t) -B(t))=Z (\sqrt{t+\Delta t}-\sqrt{t})$
Now we are left with incorporating the change of probability effect from third mode(equation 3 in post 815).
After learning more about derivatives of the diffusion SDE variable X(t) with respect to Z and B(t), I am back to solution with fokker-planck equation. I think evolution of X(t) with respect to B(t) is given as in the above copied post is the only thing we need to do in a stable fashion. The third mode or change of probability derivative might not be needed since we do change of probability with respect to gaussian/BM probability distribution numerically in our previous framework. But a partial differential equation has to solve for this change of probability in its own calculations as it gives the final result in terms of X(t) probabilities so change of probability work calculations are embedded in the PDE. We may not need to deal with all the equations in what I called Mode 3/3rd equation or $M_3(t)$ in post 815. In our previous framework, we are only concerned with how X(t) moves with B(t) and analytical change of probability as done in many terms of mode 3 is not of much interest to us since we can very easily do a numerical change of probability derivative at the end of our calculations. In the fokker-planck equation below
$\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}$
The following
$- \frac{\partial \mu(X)}{\partial X} p(X,t)+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) p(X,t)$ are clearly change of probability derivatives and may be excluded in caclulations. I am hoping that I might be able to exclude all of the mode 3/equation 3 in post 815 from my calculations.

As we see this should not be much of a problem but I want a very stable algorithm in which we do not have to numerically differentiate for higher order derivatives on every step something that makes the algorithm unstable.
After all the understanding we have gained, I hope it will not take more than a week to finalize all this work to conclusion now. Amin
Topic Author
Posts: 2152
Joined: July 14th, 2002, 3:00 am

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

Friends, my work has been going rather slowly and there were several days that I was totally unable to work. I tried to solve the equations first by finite differences but my finite difference derivatives break down after a few steps at the boundary even though I try one-sided derivatives of the same order on the boundary. Since all my work is open, I do not want to stop other friends from doing work towards a very stable algorithm for fokker-planck equation. Numerical analysis is not my forte and there are very smart people in stochastics and mathematical finance who can probably very easily tackle these issues. I have therefore decided to make a brief research note on my work and upload it on SSRN within 7-10 days on the basic work I did towards solution of Fokker Planck equation and friends who want to do further work can cite my contribution by giving reference to the research note. In the meantime, I will continue to work towards a solution of SDEs on the pattern of work I did before for mean-reverting SDEs after incorporating the understanding I gained from the derivatives of Fokker-planck. But a full-fledged solution can be several weeks away.
I would want friends to do more innovative work but would request them to cite reference to my brief research note that I will upload on SSRN in next ten days. Amin
Topic Author
Posts: 2152
Joined: July 14th, 2002, 3:00 am

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

I was able to improve my finite differences a little bit to continue them for seven to eight time intervals of  $\Delta t=.125/8$ and then I tried all sort of extreme parameters at different values of SDE starting point very close to zero and close to one and several other places. I just played with Lamperti-transformed SDE. I was able to precisely calculate the coefficients and derivatives in the main diffusion term. I was also able to closely calculate the form of second term in the evolution. About the third term, I still cannot say anything.

$X(t+\Delta t)=X(t) + \mu(X) {\Delta t}$
$+ {\sigma(X)}^2 {(\frac{\partial B}{\partial X})}^3/ \big( {(\frac{\partial B}{\partial X})}^2-.25 \frac{\partial^2 B}{\partial X^2} \big) (B(t + \Delta t) -B(t))$   Diffusion term
$+\frac{1.5}{4} {\sigma(X)}^2 ({\frac{\partial B}{\partial X}} \frac{\partial ^2 B}{\partial X^2})/({\frac{\partial B}{\partial X}}^2) \frac{(Z^2-1)}{2} {\Delta t}$  second term
$+ .5 {\sigma(X)}^2 {(\frac{\partial^3 B}{\partial X^3})} / \big( {(\frac{\partial B}{\partial X})}^2-.25 \frac{\partial^2 B}{\partial X^2} \big) \Delta t$   Third term

In case of Lamperti transformed SDE, I am confident that diffusion term has been 100% perfectly approximated by the form I have given.
I am reasonably sure that second term is also right as given. It really needs to take a $\frac{(Z^2-1)}{2}$. There is absolutely no way the derivatives in the second term can appear in the solution to SDE without multiplication with $\frac{(Z^2-1)}{2}$. It has a definitely positive sign.
About the third term, anything I have written is tentative and I cannot say anything at all with certainty since the third derivative was unstable in my finite differences but I would be confident that this term also takes a second or third hermite multiplication and possibly cannot appear just with a simple $\Delta t$
I was able to perfectly replicate the true solution to density of SDE in Lamperti coordinates for a small time upto eight time steps each of $\Delta t=.125/8$ with extreme parameters everywhere (including quite close to zero) with the given form of second term and diffusion term and neglecting the third term. By extreme parameters, I imply mean reversion greater than 2.5 and vol greater than 1.5 in original coordinates. I tried these parameters even close to zero. The results were perfect with benign parameters as well.
However general non-Lamperti form of the SDE has no simple resemblance to this solution form given above. The diffusion term in Lamperti coordinates just does not simply carry over to the general case in original coordinates as the diffusion term in original coordinates does not take the second order correction as in the case of Lamperti transformed SDE given above. Amin
Topic Author
Posts: 2152
Joined: July 14th, 2002, 3:00 am

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

I have made some changes to above post and edited it.  