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

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: 2421
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: 2421
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: 2421
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: 2421
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: 2421
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.

Cuchulainn
Posts: 62167
Joined: July 16th, 2004, 7:38 am
Location: Amsterdam
Contact:

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

I don't know if the following is interesting, but here goes:

In the early we discussed Picard iteration for ODEs. This method is also used to prove existence of SDEs and it leads to a computational,form using iterated integrals (which I now understand where they come from).

https://arxiv.org/abs/1009.5556

Don't know if it related to your approach, Amin.

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

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

Thank you Cuch for the paper. Yes it is definitely related to my work and it was a good read. I will give a more detailed review later.
I want to tell the friends about progress on my work towards the solution to Fokker-Planck equation. I keep the solution in a hermite form as following
$X(Z)=a_0(Z) + a_1(Z) H_1(Z) +1/2 a_2(Z) H_2(Z)$
$+1/6 a_3(Z) H_3(Z) +1/24 a_4(Z) H_4(Z)$  Equation(1)
The coarse derivatives can be found as
$\frac{\partial X}{\partial Z} =a_1(Z) +a_2(Z) H_1(Z)+1/2 a_3(Z) H_2(Z) +1/6 a_4(Z) H_3(Z)$
$\frac{\partial^2 X}{\partial Z^2} = a_2(Z)+a_3(Z) H_2(Z) +1/2 a_4(Z) H_2(Z)$
$\frac{\partial^3 X}{\partial Z^3} = a_3(Z) +a_4(Z) H_1(Z)$

To find precise derivatives, we will have to differentiate the coefficients $a_n$ as well. Here is the precise first derivative
$\frac{\partial X}{\partial Z} =a_1(Z) +a_2(Z) H_1(Z)+1/2 a_3(Z) H_2(Z) +1/6 a_4(Z) H_3(Z)$
$+\frac{\partial a_1(Z)}{\partial Z}H_1(Z)+.5 \frac{\partial a_2(Z)}{\partial Z} H_2(Z)$
$+1/6 \frac{\partial a_3(Z)}{\partial Z} H_3(Z) +1/24 \frac{\partial a_4(Z)}{\partial Z} H_4(Z)$

and so on for higher derivatives.
I am now doing a relatively much limited role of finite differences in this new approach. I am very sure that we will actually be able to get a worked out program for the simulation of the SDE in next week to ten days. I will write more detailed equation early next week.

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

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

Tomorrow I will be posting my experimental programs for friends after writing explanation with the code. The first program is based solely on finite differences but I was able to make ghost points outside the boundary for finite-difference derivatives. For ghost points, I took divided differences of the grid points on the boundary until divided differences continued to decrease(akin to decreasing derivatives). Then I extend the divided difference on the bottom (bottom is chosen at the point where divided differences continued to decrease) and use it to recreate divided differences next to the divided differences from calculation by adding the bottom most point to divided differences on the right. This procedure matches upto five or six derivatives for the ghost point and is more robust than anything I tried but it still does not work well after eight or nine steps in time. But I am confident that people who played with finite differences and PDEs would certainly be able to extend it to find a valid way of smoothing the derivatives since really a lot of derivatives have been matched by the divided differences ghost point creation. I hope this program will be useful.
Second program is more interesting and in this program I have converted everything in the form of hermite polynomials like mentioned in the previous post. I convert the drift into hermite polynomials and I also convert other terms into respective hermite polynomials so that solution of the SDE density remains in a hermite form at any given time. But derivatives in this procedure are only approximate since taking accurate derivatives become unstable but I am sure some friends would like to find a way around it.
I will be posing these programs tomorrow after making minor changes and adding comments to the code.

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 attaching the program for solution of SDEs in bessel or Lamperti transformed coordinates based on finite differences. This program in the original form becomes unstable after 8-12 time steps because finite differences on the boundary become unstable. Otherwise the program works fine and is very accurate. Some approximations suggested continue to work for larger time intervals. I am attaching the program and its dependent functions. dependent functions for monte carlo part are not given as they have been uploaded in earlier posts. I will be posting other versions of the program after a few days.

function [] = FokkerPlanckNew02Wilmott()
%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
%Bessel process/Lamperti transformed version of the SDE and then changed coordinates to retreive
%the SDE in original coordinates.
%The simulation process is given as
%w(t+1)=w(t)+mu(w)dt+sigma0^2 dBdw^3/(dBdw^2-.25*d2Bdw2)*Z*(sqrt(t+dt)-sqrt(t))
%+1.5/4 dBdw * d2Bdw2* (Z^2-1)/2 *sigma0^2*dt.
%I have neglected the third hermite polynomial term with third derivative.
%That neglected term should probably be
%=.5/4*sigma0^2 *d3Bdw3*(Z^3-3*Z)/6 *dt;

dt=.125/8;   % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=8;     % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
%Tt2=0;%12*1-1;
OrderA=2;  %Keep at two.
OrderM=2;  %Keep at two.
dtM=dt;%Monte carlo time interval size dtM.
TtM=Tt;%Monte carlo number of simulation intervals.
x0=1.0;   % starting value of SDE
beta1=0.0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.750;   % volatility power.
kappa=1.0;%1.950;   %mean reversion parameter.
theta=1.0;   %mean reversion target
mu1=1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
sigma0=.840;%Volatility value
alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;

%For most stable grid, I choose 50 subdivisions with dZ=.2
dNn=.20;   % Normal density subdivisions width. would change with number of subdivisions
Nn=50;%120;%  % No of normal density subdivisions
NnMidl=25;%60;%One half left from mid of normal density(low)
NnMidh=26;%61;%One half right from the mid of normal density(high)
NnMid=4.0;
Z(1:Nn)=(((1:Nn)-5.50)*dNn-NnMid);

%In order to choose finer grid with 100 subdivisions with dZ=.1, comment
%the above block and decomment the block below
%dNn=.10;   % Normal density subdivisions width. would change with number of subdivisions
%Nn=100;%120;%  % No of normal density subdivisions
%NnMidl=50;;%One half left from mid of normal density(low)
%NnMidh=51;;%One half right from the mid of normal density(high)
%NnMid=4.0;
%Z(1:Nn)=(((1:Nn)-10.50)*dNn-NnMid);

%In order to choose even finer grid with 200 subdivisions with dZ=.05, comment
%the above block and decomment the block below
%dNn=.05;   % Normal density subdivisions width. would change with number of subdivisions
%Nn=200;%120;%  % No of normal density subdivisions
%NnMidl=100;%One half left from mid of normal density(low)
%NnMidh=101;%One half right from the mid of normal density(high)
%NnMid=4.0;
%Z(1:Nn)=(((1:Nn)-20.50)*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

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

wnStart=1;%\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\1;
tic

for tt=1:Tt

[dwdZ,d2wdZ2] = InterpolateDerivativesDivDiffs(wnStart,Nn,dNn,w,Z);

dZdw(1:Nn)=1.0./dwdZ(1:Nn);
d2Zdw2(1:Nn)=-1.*d2wdZ2(1:Nn).*dZdw(1:Nn).^3;%Inverse second derivative formula
dwdB(1:Nn)=dwdZ(1:Nn)./sqrt((tt-1)*dt);%convert to derivative w.r.t brownian motion

dBdw(1:Nn)=1.0./dwdB(1:Nn);

d2wdB2(1:Nn)=0;
d2Bdw2(1:Nn)=0;
if(tt>3)
d2wdB2(1:Nn)=d2wdZ2(1:Nn)./((tt-1)*dt);;%convert to derivative w.r.t brownian motion

d2Bdw2(1:Nn)=-1.*d2wdB2(1:Nn).*dBdw(1:Nn).^3; %Inverse second derivative formula
end

tt
plot((1:Nn),dBdw(1:Nn),'r');
str=input('Look at first derivative')

plot((1:Nn),d2Bdw2(1:Nn),'r');
str=input('Look at second derivative')

%end

[wMu0dt]=CalculateDrift(mu1,mu2,beta1,beta2,gamma,sigma0,w,wnStart,Nn,dt);

w(isnan(w)==1)=0;

wMu0dt(isnan(wMu0dt)==1)=0;

if(tt==1)
dBdw(wnStart:Nn)=1/sigma0;
end
if(tt>=1)

w(wnStart:Nn)=w(wnStart:Nn)+wMu0dt(wnStart:Nn)+  ...
+1*1.5/4*sigma0^2.*dBdw(wnStart:Nn).*d2Bdw2(wnStart:Nn)./(dBdw(wnStart:Nn).^2)*dt.*((Z(wnStart:Nn).^2-1)/2).^1+ ...
sigma0^2.*dBdw(wnStart:Nn).^3./(dBdw(wnStart:Nn).^2-1*.25*d2Bdw2(wnStart:Nn)).*Z(wnStart:Nn).*(sqrt(tt*dt)-sqrt((tt-1)*dt));% + ...

%%An approximation is below
%w(wnStart:Nn)=w(wnStart:Nn)+wMu0dt(wnStart:Nn)+  ...
%    sigma0^2.*dBdw(wnStart:Nn).^3./(dBdw(wnStart:Nn).^2-0*.25*d2Bdw2(wnStart:Nn)).*Z(wnStart:Nn).*(sqrt(tt*dt)-sqrt((tt-1)*dt));
%Another approximation is below
%w(wnStart:Nn)=w(wnStart:Nn)+wMu0dt(wnStart:Nn)+  ...
%    +1*1.5/4*sigma0^2.*dBdw(wnStart:Nn).*d2Bdw2(wnStart:Nn)./(dBdw(wnStart:Nn).^2)*dt.*((Z(wnStart:Nn).^2-1)/2).^1+ ...
%    sigma0^2.*dBdw(wnStart:Nn).^3./(dBdw(wnStart:Nn).^2-0*.25*d2Bdw2(wnStart:Nn)).*Z(wnStart:Nn).*(sqrt(tt*dt)-sqrt((tt-1)*dt));

end
%    if(tt==1)
%
%    w(wnStart:Nn)=w(wnStart:Nn)+wMu0dt(wnStart:Nn)+  ...
%        sigma0.*sqrt(tt*dt).*Z(wnStart:Nn);
%
%    end

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

%     w1(1:Nn-1)=w(1:Nn-1);
%     w2(1:Nn-1)=w(2:Nn);
%     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.

u0=theta+(x0-theta)*exp(-kappa*(tt*dt)); %analytic mean of the density
%If you are not using mean reverting 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);

%Uncomment the line below for mean correction to work.

%     w(wnStart:Nn)=wCorrected(wnStart:Nn);

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

%plot(w(1:Nn),fw(1:Nn),'r',wOld(1:Nn),fwOld(1:Nn),'b');
%str=input('Look at graphs');

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

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.
%Pre-compute the time and power exponent values in small multi-dimensional arrays
for k = 0 : (OrderA+1)
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));

end
end
end
end
%Below call the program that calculates the Ito-Taylor coeffs.
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);
%Yq=ItoTaylorCoeffsNew(alpha1,beta1,beta2,gamma);
%YqCoeff=Yq/(1-gamma);%These coordinates are divided by (1-gamma)
%Calculate the density with monte carlo below.
%rng(85109634, 'twister')
rng(96348510, 'twister')
paths=250000;
YY(1:paths)=x0;  %Original process monte carlo.
% YY2dt(1:paths)=0.0;
%YY4dt(1:paths)=0.0;
%X(1:paths)=x0^(1-gamma)/(1-gamma);
Random1(1:paths)=0;
for tt=1:TtM
%    t1=tt*dtM;
%    t0=(tt-1)*dtM;
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;
YY0=YY;
%   YYPrev=YY;
%   XX0=X;

for k = 0 : OrderM
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;
if(l4<=5)

%                    X(1:paths)=X(1:paths) + YqCoeff(l1,l2,l3,l4) * ...
%                   mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
%                   ((1-gamma)*XX0(1:paths)).^(Fp1(l1,l2,l3,l4)/(1-gamma)) .* ...
%                   HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);

YY(1:paths)=YY(1:paths) + YCoeff(l1,l2,l3,l4) * ...
mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
YY0(1:paths).^Fp(l1,l2,l3,l4) .* HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);
%Above is original coordinates monte carlo equation.
end
end
end
end
end
end
YY(YY<0)=0;

sum(YY(:))/paths %origianl coordinates monte carlo average.
sum(y_w(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
% sum(y_w(5:95).*ZProb(5:95)) %Original process average from coordinates
y_w(isnan(y_w)==1)=0.0;
sum(y_w(wnStart:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
%sum(y_w(4:55).*ZProb(4:55))
theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average

BinSize=.0075;
if(theta<=.26)
BinSize=.0075;%Please change this bin size accordingly. When data range is too small decrease it.\
end

if(theta<=.15)
BinSize=.0075/5;%Please change this bin size accordingly. When data range is too small decrease it.\
end
if(theta>=.95)
BinSize=.0075*2;
end

if(theta>=10)
BinSize=.075*1;
end
if(theta>=50)
BinSize=.075*5;
end
if(theta>=100)
BinSize=.075*10;
end
if(theta<=.025)
BinSize=.0075/12;
end
%Please change this bin size accordingly. When data range is too small decrease it.
%When density range is too large increase it.Roughly, When starting point of SDE is
%x0=1, set at .0075*2; when starting point is x0=.1, set at .0075/2.0;

MaxCutOff=3000;
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
function [wMu0dt]=CalculateDrift(mu1,mu2,beta1,beta2,gamma,sigma0,w,wnStart,Nn,dt)

dt2=dt^2/2;
dt3=dt^3/6.0;
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
theta1=mu1;
theta2=mu2;
theta3=.5*sigma0^2*(-gamma);
omega1=beta1-gamma;
omega2=beta2-gamma;
omega3=gamma-1;
GG(wnStart:Nn)=theta1*yy(wnStart:Nn).^omega1+theta2*yy(wnStart:Nn).^omega2+theta3*yy(wnStart:Nn).^omega3;
dGG(wnStart:Nn)=omega1*theta1*yy(wnStart:Nn).^(omega1-1)+omega2*theta2*yy(wnStart:Nn).^(omega2-1)+ ...
omega3*theta3*yy(wnStart:Nn).^(omega3-1);
DD(wnStart:Nn)=mu1*yy(wnStart:Nn).^beta1+mu2*yy(wnStart:Nn).^beta2;
dGG2(wnStart:Nn)=theta1*omega1*(omega1-1).*yy(wnStart:Nn).^(omega1-2)+ ...
theta2*omega2*(omega2-1).*yy(wnStart:Nn).^(omega2-2)+ ...
theta3*omega3*(omega3-1).*yy(wnStart:Nn).^(omega3-2);
QQ(wnStart:Nn)=sigma0^2*yy(wnStart:Nn).^(2*gamma);
ddGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*yy(wnStart:Nn).^(omega1+beta1-2)+ ...
mu1*omega2*theta2*(omega2+beta1-1).*yy(wnStart:Nn).^(omega2+beta1-2)+ ...
mu1*omega3*theta3*(omega3+beta1-1).*yy(wnStart:Nn).^(omega3+beta1-2)+ ...
mu2*omega1*theta1*(omega1+beta2-1).*yy(wnStart:Nn).^(omega1+beta2-2)+ ...
mu2*omega2*theta2*(omega2+beta2-1).*yy(wnStart:Nn).^(omega2+beta2-2)+ ...
mu2*omega3*theta3*(omega3+beta2-1).*yy(wnStart:Nn).^(omega3+beta2-2);
d2dGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*(omega1+beta1-2).*yy(wnStart:Nn).^(omega1+beta1-3)+ ...
mu1*omega2*theta2*(omega2+beta1-1).*(omega2+beta1-2).*yy(wnStart:Nn).^(omega2+beta1-3)+ ...
mu1*omega3*theta3*(omega3+beta1-1).*(omega3+beta1-2).*yy(wnStart:Nn).^(omega3+beta1-3)+ ...
mu2*omega1*theta1*(omega1+beta2-1).*(omega1+beta2-2).*yy(wnStart:Nn).^(omega1+beta2-3)+ ...
mu2*omega2*theta2*(omega2+beta2-1).*(omega2+beta2-2).*yy(wnStart:Nn).^(omega2+beta2-3)+ ...
mu2*omega3*theta3*(omega3+beta2-1).*(omega3+beta2-2).*yy(wnStart:Nn).^(omega3+beta2-3);
ddGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*yy(wnStart:Nn).^(omega1+2*gamma-3)+ ...
sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*yy(wnStart:Nn).^(omega2+2*gamma-3)+ ...
sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*yy(wnStart:Nn).^(omega3+2*gamma-3);
d2dGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*(omega1+2*gamma-3).*yy(wnStart:Nn).^(omega1+2*gamma-4)+ ...
sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*(omega2+2*gamma-3).*yy(wnStart:Nn).^(omega2+2*gamma-4)+ ...
sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*(omega3+2*gamma-3).*yy(wnStart:Nn).^(omega3+2*gamma-4);

VV(wnStart:Nn)=sigma0.*yy(wnStart:Nn).^gamma;

wMu0dt(wnStart:Nn)=GG(wnStart:Nn)*dt + dGG(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...
.5*dGG2(wnStart:Nn).*QQ(wnStart:Nn).*dt2+ ...
ddGGDD(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
.5*ddGG2QQ(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
.5*d2dGGDD(wnStart:Nn).*QQ(wnStart:Nn).*dt3 + ...
.25*d2dGG2QQ(wnStart:Nn).*QQ(wnStart:Nn).*dt3;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [dwdZ,d2wdZ2] = InterpolateDerivativesDivDiffs(wnStart,Nn,dNn,w,Z)
POrder=10;
%Below we invert the w array at the right boundary to find ghost point at the end boundary.
mm=0;
for nn=Nn:-1:Nn-POrder+1
mm=mm+1;
wInv(mm)=w(nn);
end

%wGhostPointE is ghost point at the end boundary.
%wGhostPointS is ghost point at the start boundary.
%Both use a slightly different function.
[wGhostPointE] = ExtrapolateGhostPointWithDivDiffs(wInv,dNn,POrder);
[wGhostPointS] = ExtrapolateGhostPointWithDivDiffs02(wnStart,Nn,dNn,w,POrder);

%Below is the calculation of first and second derivatives.

dwdZ(wnStart)=(-1*wGhostPointS+1*w(wnStart+1))/(2*dNn);
dwdZ(wnStart+1:Nn-1)=(-1*w(wnStart:Nn-2)+1*w(wnStart+2:Nn))/(2*dNn);
dwdZ(Nn)=(wGhostPointE-1*w(Nn-1))/(2*dNn);

d2wdZ2(wnStart)=(wGhostPointS+w(wnStart+1)-2*w(wnStart))/(dNn.^2);
d2wdZ2(wnStart+1:Nn-1)=(w(wnStart:Nn-2)-2*w(wnStart+1:Nn-1)+1*w(wnStart+2:Nn))/(dNn.^2);
d2wdZ2(Nn)=(wGhostPointE+w(Nn-1)-2*w(Nn))/(dNn.^2);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [GhostPoint] = ExtrapolateGhostPointWithDivDiffs(w,dNn,POrder)
%POrder is the order till which divided differences will continue.
DMat=zeros(POrder,POrder);
Dmat(1,1:POrder)=w(1:POrder);
Order2=POrder;
EndFlag=0;
for nn=1:POrder-1
Order2=Order2-1;
for mm=1:Order2
Dmat(nn+1,mm)=(Dmat(nn,mm+1)-Dmat(nn,mm))/dNn;
end

if((abs(Dmat(nn+1,1))<.000000001) && (EndFlag==0))
EndFlag=1;
FinalOrder=nn;
end
end
%FinalOrder is the order of divided differences used to extrapolate the ghost
%point on boundary.

if(EndFlag==0)
FinalOrder=POrder;
end

EndFlag2=0;
for nn=1:POrder-1
if((abs(Dmat(nn+1,1))>abs(Dmat(nn,1)))&&(EndFlag2==0))
FinalOrder=nn-1;
EndFlag2=1;
end
end
if(FinalOrder<2)
FinalOrder=2;
end
FinalOrder
%Below we construct the boundary ghost point.
XX(FinalOrder)=Dmat(FinalOrder,1)+Dmat(FinalOrder,1)-Dmat(FinalOrder,2);
for nn=FinalOrder-1:-1:1
XX(nn)=-XX(nn+1).*dNn+Dmat(nn,1);
end

GhostPoint=XX(1);

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [GhostPoint] = ExtrapolateGhostPointWithDivDiffs02(wnStart,Nn,dNn,w,POrder)
%POrder is the order till which divided differences will continue.
DMat=zeros(POrder,POrder);
Dmat(1,1:POrder)=w(wnStart:wnStart+POrder-1);
Order2=POrder;
EndFlag=0;
for nn=1:POrder-1
Order2=Order2-1;
for mm=1:Order2
Dmat(nn+1,mm)=(Dmat(nn,mm+1)-Dmat(nn,mm))/dNn;
end

if((abs(Dmat(nn+1,1))<.000000001) && (EndFlag==0))
EndFlag=1;
FinalOrder=nn;
end
end
%FinalOrder is the order of divided differences used to extrapolate the ghost
%point on boundary.

if(EndFlag==0)
FinalOrder=POrder;
end

EndFlag2=0;
for nn=1:POrder-1
if((abs(Dmat(nn+1,1))>abs(Dmat(nn,1)))&&(EndFlag2==0))
FinalOrder=nn-1;
EndFlag2=1;
end
end
if(FinalOrder<1)
FinalOrder=1;
end
FinalOrder
%Below we construct the boundary ghost point.
XX(FinalOrder)=Dmat(FinalOrder,1)+Dmat(FinalOrder,1)-Dmat(FinalOrder,2);
for nn=FinalOrder-1:-1:1
XX(nn)=-XX(nn+1).*dNn+Dmat(nn,1);
end

GhostPoint=XX(1);

end



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 written a program for two dimensional correlated stochastic volatility models in my framework which I used for 1D Lamperti transformed SDEs. I am testing the code and comparing it with two dimensional monte carlo. I hope to post a full-fledged and very general stochastic volatility pricing program within a few days. I hope most friends would like this new program. Should not be more than a few days to finesse everything.

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 sorry for the delay. But before we can precisely calculate two dimensional density in a stochastic volatility framework, we have to first write the program for precise calculation of dt-integrals and dz-integrals from the SDE evolution. I had earlier presented the idea of Ito-calculus on the SDE grid and even posted a program. Unfortunately I made several errors in that program. So despite that the original idea was sound, the resulting densities were widely off.. I have fixed the errors in the previous program. I had earlier made the wrong assumption that transition probabilities between any two points on the SDE can be copied from the Brownian motion transition probabilities. This was utterly wrong. I was led to this false conclusion because bessel process and brownian motion densities and their transition probabilities are very close to each other, at least, at the start of the processes. We have to calculate the transition densities from the bessel process in order to get precise results. Similarly I have fixed previous errors in the calculation of effective Z or normal increment between any two points on the SDE density to do Ito-calculus/change of variable on the SDE density grid. I have removed the errors and now I move towards calculation of all functionals of the SDE and their evolution on the same grid on which the SDE evolves in time. This is very interesting and helps us quickly calculate the dt-integrals and dz-integrals along the evolution of the SDE with great precision. Once the densities of these integrals are handily available, they can be converted into two dimensional process densities quite easily in several ways. In two to three days, I will be posting the program for calculation of the densities of functionals of SDE on the same grid on which the SDE is evolving. This is very interesting and I hope most friends would like this new 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

OK friends, I wrote a different program for solution of densities of 1D SDEs. I took a fixed grid with 700 cells and pre-calculated transition probabilities and then evolved the SDE. This way you can almost solve for the density of any SDE. zero remains a minor problem but I will fix that. This grid works for any SDE including SDEs other than mean reverting types. Now on this grid, we can calculate forward expectation of any function for which we can write Ito change of variable formula. We then convert expectations of the function into function values and calculate its density. When density itself is known there is no point in simple Ito change of variable since we can analytically transform the density at terminal point and convert it into the density of the function of SDE variable. But dt-integrals cannot be solved by any simple transformation and they are also precisely calculated in my framework on a fixed grid/lattice. Forward expectation or Ito change of variable in a grid is not previously known in mathematics and it is a new thing. I have used a large grid but since most of the variables have been pre-calculated for each grid point, there is a main loop in which a fast computer running in C++ will easily zap the whole thing in 150 ms for SDE and its functionals though it takes me 2-3 seconds in matlab on my slow laptop.
We can use this quadratic variation in calculation of solution of 2D SDEs. We can also calculate the quadratic variations related to variance swaps and other variance/vol instruments once the volatility/variance SDE has been specified. We can easily know what is the expected quadratic variation associated with each point on the terminal variance density. I also believe that we can find forward expectation of other European type derivatives like swaps on this lattice but for that we will have to specify proper transition variances. Another application will be finding the density of the integrated asset for Asian options on a 1D lattice/grid framework.  I will post my code tomorrow or on Saturday after writing comments and explanations.
Meanwhile, as I have been previously saying that my programs are transported to some people in US, if somebody recently mentions that they have done the same thing like Ito change of variable formula or forward expectation on a fixed lattice/grid with transition probabilities, please know that they have seen my programs and that is what they are talking about. So please beware of the thieves who want to steal my 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

In this program, I am simulating the SDE given as
dyy(t)=mu1 yy(t)^beta1 dt + mu2 yy(t)^beta2 dt +sigma yy(t)^gamma dz(t)
I have not directly simulated the SDE but simulated the transformed  Bessel process version of the SDE.
There is one caveat as I used the process for fast mean reverting SDEs and I specified the mean for such mean reverting SDEs which is commonly known
but If true mean is not known in original coordinates just comment the appropriate if loop in the program. The program will work perfectly well for all SDEs away from zero. Just comment the appropriate if loop for other non-mean reverting type SDEs.
Currently, it uses a mean from mean-reverting SV models. Current mean correction is helpful for processes that go closer to zero. In the next version, I will calculate true mean for every SDE and then generalize this mean correction aspect perfectly to all SDEs.
Disclaimer: The program is very general and can be used for any SDE that does not go to zero very fast. this program works perfectly well for mean reverting
diffusions with high volatility and high drift for x0>=.1 and theta>=.1 and gamma(vol exponent)>.55;I have not checked for all possible type of SDEs but the above guideline would apply. The program is extremely precise when the SDE process does not get to zero. Conditions on are more restrictive for the Ito process of the SDE.  I have already made some corrections to make it better close to zero but due to limited time for this first version, I did not have time to make it perfect close to zero. It is doable and in the next version, I will make it perfectly well very close to zero.

For the Ito functions on the transition probability grid, first Ito function is f=yy^omega1,and second Ito function being calculated is fdt=yy^omega1*dt,
So the current program simulates three processes, two Ito processes and one original SDE. Simulation of Ito processes takes more time than simulating the original SDE process.
For both the Ito processes, the value omega1 has to be greater than .75  when gamma is smaller than .65; when gamma is grerater than .75, you can choose omega1 close to .5. for gamma=.95 and close, you can choose any value of omega1 since the  process does not reach zero. The reason is that when gamma is smaller than .65 and omega1 is .5, the bessel process becomes unstable close to zero.It will be perfectly fixed in the next version but it requires some more work. Very small instability or inaccuracy becomes explosive in the transition probabilities framework and can make the whole process inaccurate. I did not work on all possible values of omega1 due to limited time for this first version. Next version will cover this issue.
I checked the program for one year only. Though I know from a few runs that  original SDE works perfectly fine once it remains stable for one year and remains away from zero or does not go to zero too fast, Ito processes may lose accuracy since in this version, I did not use higher order generators for Ito processes so low order generators lose slight accuracy with time. Next version will fix this problem.

I create a fixed grid of size Nn. The grid is based on bessel/lamperti version of the process.Bessel process very approximately remains closer to a gaussian process so it also helps in calculating a grid that is relatively uniform in Z which is the driving noise of the SDE. I calculate the drift and volatility of the SDE and then calculate a Zt(driving noise) which links the starting point of the SDE to neighboring grid points.
In the program the state dependent drift is called wMu0dt. if w1 is starting point on grid and w2 is final point and state dependent volatility(defined here as coefficient of Z in expansion of SDE)is dw0(w1),  then Zt(w1,w2)=(w2-w1-wMu0dt)/dw0(w1);This is the Zt(w1,w2) which will be used for transition probabilities calculation and later for Ito on the fixed grid as the SDE process moves from w1 to w2. This Zt(w1,w2) is the same as Zt(yy1,yy2) where w is the bessel process and yy is the process in the original coordinates so Zt remains the same across original and bessel coordinates.
For the SDE probabilities calculation, I calculate their evolutions from point w1 at time to to w2 at time t+dt as
P(w2)=Sum(over w1) Pt(w1,w2)*P(w1);
Adding time index explicitly to above equation, we get
P(w2(t+dt))=Sum(over w1(t)) Pt(w1(t),w2(t+dt))*P(w1(t));
Here P(w1) is the occupation probability of point w1 at time t.
And P(w2) is the occupation probability which is being calculated at time t+dt.
Please know that all of these probabilites are integrated probability mass and not raw probabilities. Transition probabilities and occupation probabilities are calculated to be integrated probability mass over the appropriate grid intervals. Calculating the drift from bessel version is very accurate but very closer to
zero the bessel version is very inaccurate and relatively unstable so there I use the original coordinates to determine Zt and associted transition
probabilities. Transition probabilities only depend on Zt (the driving noise). below(and also in the program code) I have used yy as the SDE process in original coordinates.
When evolving the Ito functionals of the diffusion, we do not evolve their probabilities but evolve thier expectations for the appropriate grid cells. Since the correct expectation is known at every grid cell, and the probability is known from the original process, we can calculate the appropriate value for the Ito-process at the terminal time by dividing the expectaation at that cell by integrated probability mass in that cell.
I have not calculated arbitrary payoffs along the process path since they require a specification of forward variance at each cell in forward point so you still cannot do swaps or arithmetic addition along the process. with the current program, You can only evolve a proper Ito process starting at initial time.
Next week I will release another version where forward variance will be specified and all forward coupons and swap type payoffs could be specified. And arithmetic asian type options would also be doable after proper specification of forward variance. I simply did not have the time to complete all of this in this initial version.
Here is the main program.
function [] = SDEandFunctionalsDensityWilmott()
%Disclaimer: The program is very general and can be used for any SDE that
%does not go to zero.
%this program works perfectly well for mean reverting
%diffusions with high volatility and high drift for x0>=.1
%and theta>=.1 and gamma(vol exponent)>.55;I have not checked for all
%possible type of SDEs but the above guideline would apply. The program is
%extremely precise when the SDE process does not get to zero.
%Conditions on are more restrictive for the Ito process of the SDE.
%I have already made some corrections to make it better close to zero but
%due to limited time for this first version, I did not have time to make it
%perfect close to zero. It is doable and in the next version, I will make
%it perfectly well very close to zero.
%In this program, I am simulating the SDE given as
%dyy(t)=mu1 yy(t)^beta1 dt + mu2 yy(t)^beta2 dt +sigma yy(t)^gamma dz(t)
%I have not directly simulated the SDE but simulated the transformed
%Bessel process version of the SDE and then changed coordinates to retreive
%the SDE in original coordinates.
%The present program will directly simulate only the Bessel Process version of the
%SDE in transformed coordinates.
%There is one caveat as I used the process for fast mean reverting SDEs
%I specified the mean for such mean reverting SDEs which is commonly known
%but If true mean is not known in original coordinates just comment
%the appropriate if loop in the program. The program will work perfectly well for all
%SDEs away from zero. Just comment the appropriate if loop for other SDEs.
%Currently, it uses a mean from mean-reverting SV models. Current mean
%correction is helpful for processes that go closer to zero. In the next
%version, I will calculate true mean for every SDE and then generalize this
%mean correction aspect perfectly to all SDEs.

%For the Ito functions, first Ito function is f=yy^omega1,and second Ito
%function being calculated is fdt=yy^omega1*dt,
%So the current program simulates three processes. Simulation of Ito
%processes takes more time than simulating the original process.
%For both the Ito processes, the value omega1 has to be greater than .75
%when gamma is smaller than .65;
%when gamma is grerater than .75, you can choose omega1 close to .5.
%for gamma=.95 and close, you can choose any value of omega1 since the
%process does not reach zero.
%The reason is that when gamma is smaller than .65 and omega1 is .5, the
%bessel process becomes unstable close to zero.It will be perfectly fixed
%in the next version. Very small instability or inaccuracy becomes explosive
%in the transition probabilities framework and can make the
%whole process inaccurate. I did not work on all possible values of omega1
%due to limited time for this first version. Next version will cover this
%issue.
%I checked the program for one year only. Though I am confident that
%original SDE will work perfectly fine once it remains stable for one year
%and remains away from zero, Ito processes may lose accuracy since in the
%original version, I did not use higher order generators for Ito processes.
%Next version will fix this problem.

%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%or skype ahsan.amin2999
%I create a fixed grid of size Nn. The grid is based on bessel/lamperti
%version of the process.Bessel process very approximately remains closer to
%a gaussian process so it also helps in calculating a grid that is
%relatively uniform in Z which is the driving noise of the SDE. I calculate
%the drift and volatility of the SDE and then calculate a Zt(driving noise) which links
%the starting point of the SDE to neighboring grid points.
%In the program the state dependent drift is called wMu0dt. if w1 is
%starting point on grid and w2 is final point and state dependent
%volatility(defined here as coefficient of Z in expansion of SDE)is dw0(w1),
%then Zt(w1,w2)=(w2-w1-wMu0dt)/dw0(w1);This is the Zt(w1,w2) which will be used for
%transition probabilities calculation and later for Ito on the fixed grid as the
%SDE process moves from w1 to w2. This Zt(w1,w2) is the same as Zt(yy1,yy2) where
%w is the bessel process and yy is the process in the original coordinates
%so Zt remains the same across original and bessel coordinates.
%For the SDE probabilities calculation, I calculate their evolutions from
%point w1 at time to to w2 at time t+dt as
%P(w2)=Sum(over w1) Pt(w1,w2)*P(w1);
%Adding time to above equation, we get
%P(w2(t+dt))=Sum(over w1(t)) Pt(w1(t),w2(t+dt))*P(w1(t));
%Here P(w1) is the occupation probability of point w1 at time t.
%And P(w2) is the occupation probability which is being calculated at time
%t+dt.
%Please know that all of these probabilites are integrated probability mass
%and not raw probabilities. Transition probabilities and occupation
%probabilities are calculated to be integrated probability mass over the
%appropriate grid intervals.
%Calculating the drift from bessel version is very accurate but very closer to
%zero the bessel version is very inaccurate and relatively unstable so
%there I use the original coordinates to determine Zt and associted transition
%probabilities. Transition probabilities only depend on Zt (the driving noise).
%below(and also in the program code) I have used yy as the SDE process in original coordinates.

%When evolving the Ito functionals of the diffusion, we do not evolve their
%probabilities but evolve thier expectations for the appropriate grid
%cells. Since the correct expectation is known at every grid cell, and the
%probability is known from the original process, we can calculate the
%appropriate value for the Ito-process at the terminal time by dividing the
%expectaation at that cell by integrated probability mass in that cell.
%I have not calculated arbitrary payoffs along the process path since they
%require a specification of forward variance at each cell in forward point
%so you still cannot do swaps or arithmetic addition along the process. You
%can only evolve a proper Ito process starting at initial time.
%Next week I will release another version where forward variance will be
%specified and all forward coupons and swap type payoffs could be
%specified. And arithmetic asian type options would also be doable after proper
%specification of forward variance. I simply did not have the time to
%complete all of this in this initial version.
%I have simulated the dt-integral given as an ito process as
% \int_0^t yy(s)^omega1 ds  %%%here omega1 is power on yy(s)
%for this I used the true dt-formula on each step given as
%\int_t0^t1 yy(s)^omega1 ds=t1* yy(t1)^pmega1- t0* yy(t0)^omega1 - ...
% omega1 * yy(t0)^(omega1-1)*(mu1*yy(t0)^beta1 + mu2*yy(t0)^beta2) *(t1^2-t0^2)/2 -...
% omega1 * yy(t0)^(omega1-1)*(sigma0 *yy(t0)^gamma)* N * (t1^1.5-t0^1.5)/sqrt(3) -...
%.5*omega1*(omega1-1)*sigma0^2*yy(t0)^(omega1+2*gamma-2)*(t1^2-t0^2)/2.0;

%I have used the same simulation formula in monte carlo and my
%transition probabilities framework. I have used original coordinates
%both in monte carlo and transition probabilities framework for calculation
% of dt-integrals.
% REference for the method for path-integrals/dt-integral I used in my
%program is my paper on
% SSRN titled "Solution of Stochastic Volatility Models Using Variance
%Transition Probabilities and Path Integrals"
%https://ssrn.com/abstract=2149231 or
%http://dx.doi.org/10.2139/ssrn.2149231
%For diffusions going into zero, you may see instabilities, I have not
%corrected for that yet in this initial version.

dt=.125/8;   % Simulation time interval.%For diffusions close to zero
%decrease dt for accuracy.
Tt=8*8;     % Number of simulation levels. Terminal time= Tt*dt; //.125/8*8*8=1 year;
OrderA=2;  %For monte carlo order, Keep at two.
OrderM=2;  %For monte corlo order,Keep at two.OrderM has to be equal to OrderA
%is the true constraint otherwise you can change the order upto 4.

x0=1;   % starting value of SDE in original coordinates
beta1=0.0;   % the first drift term power.
beta2=1.0;   % Second drift term power.
gamma=.850;   % volatility power.
kappa=1.0;%1.450;   %mean reversion parameter.
theta=1;   %mean reversion target
%Below you can get rid of kappa and theta and freely choose your own values
%for both drift coefficients.
mu1=1*theta*kappa;   % first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
sigma0=.950;%Volatility value
omega1=.750; %power on x(t) in dt-integral. This is also the power for the
%process which is evolved using Ito formula on the grid.
%I have tested powers between .65 and 1;
%Omega1 has to be greater than .75 when gamma is smaller than .65;
%when gamma is grerater than .75, you can choose omega1 close to .5.
%The reason is that when gamma is smaller than .65 and omega1 is .5, the
%bessel process becomes unstable close to zero.It will be perfectly fixed
%in the next version. Very small instability or inaccuracy can make the
%whole process inaccurate. I did not work on it due to limited time for
%this first version.

theta1=1;%Keep theta1 at one. It will be applicable in the next version.
%Changing theta1 will give erroneous results.
%I have used omega1 and theta1=1 for Ito process simulation.
%that is that I simulate yy(t)^omega1 as an Ito process of yy(t) of
%the SDE in original coordinates.
%may fail for extrme high or low powers.

alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
%Do not change above alpha.

tic
w0=x0^(1-gamma)/(1-gamma);
w0T=w0;
for tt=1:Tt
w0T=w0T+CalculateDrift(mu1,mu2,beta1,beta2,gamma,sigma0,w0T,1,1,dt);
end
%Above I made an approximate calculation of median of the density in bessel
%coordinates at the terminal time.

wS1=(w0T-sigma0*sqrt(Tt*dt)*3.5);%wS1 means starting point 1 for the bessel
%grid which is kept at 3.5 standard normal deviations below the projected
%median. The argument is that bessel standard deviations are roughly close
%to normal standard deviations.Obviously this does not apply to
%mean-reverting type processes and appropriate integrated volatility could
%be used.
wS2=(w0-sigma0*sqrt(Tt*dt)*3.5);%wS2 second candidate for the starting point of the grid.
if(wS1<wS2)
wS0=wS1;
else
wS0=wS2;
end
%Above wS0 is the true lower starting point of the grid chosen from the
%lesser of wS1 or wS2.
if(wS0<0)
wS0=.0001;
end
%wE1 and wE2 are the candidates for the end point of the grid and we choose
%one of them below as temporary end point of the grid given by wE0.
wE1=(w0T+sigma0*sqrt(Tt*dt)*3.5);
wE2=(w0+sigma0*sqrt(Tt*dt)*3.5);
if(wE1>wE2)
wE0=wE1;
else
wE0=wE2;
end

Nn=801; %Total numer of nodes in the grid. hereit is approximate.
dNn0=(wE0-wS0)/(Nn-1);%temporary calculation of grid interval.
Nn1=(w0-wS0)/dNn0;%temporary calculation of the grid size.
%Nn=Nn1;
dNn=(w0-wS0)/floor(1+Nn1);%True grid interval so that grid starting point and the
%SDE origin w0 both lie on the grid.
Nn=floor(1+(wE0-wS0)/dNn);%True grid size clculated again slightly
%different from above temporary grid size.
wE0=wS0+Nn*dNn; %True grid end point
w(1:Nn)=wS0+((1:Nn)-1)*dNn;%Calculate each grid cell centre.
wS(1:Nn)=w(1:Nn)-.5*dNn;%Grid starts;
wE(1:Nn)=w(1:Nn)+.5*dNn;%Grid Ends.

cN=1+ceil(sqrt(dt)*sigma0*5.5/dNn);
%Above cN calculates the 5.5 standard deviations distance for one step
%volatility. This distance will applied so that we do not have to do
%transition probabilities calculation for the entire grid but only from
%below 5.5 SDs to positive 5.5 SDs.
%w0
Nn0=round(1+(w0-wS0)/dNn);%This is the cell associated with delta origin
%that will have occupation probability one at the start.
w(Nn0)
wS0
wE0
str=input('Look at numbers')

yy(1:Nn)=((1-gamma)*w(1:Nn)).^(1/(1-gamma));%Original coordinates for each cell.
fyy(1:Nn)=theta1*yy(1:Nn).^omega1;
%w
b1(1:Nn)=w(1:Nn);
b2(1:Nn)=w(1:Nn);
wnStart=1;
Ns=65;  %Ns is used at the starting-N. This is where I calculate Zt from
%original coordinates for accuracy. 65 is arbitrary but accurate enough here.

%Below we calculate the drift wMu0dt1 and vol(inclusive of sqrt(dt)) in original coordiates
%This is more accurate for first Ns grid points.
[wMu0dt1,dw01]=CalculateDriftAndVolOriginalCoordinates(mu1,mu2,beta1,beta2,gamma,sigma0,w,1,Nn,dt);
%Below we calculate the drift and vol in bessel coordiates. This is
%precision accurate away from zero.
[wMu0dt2,dw02]=CalculateDriftAndVol(mu1,mu2,beta1,beta2,gamma,sigma0,w,1,Nn,dt);
% wMu0dt(1:Ns)=wMu0dt0(1:Ns);%First Ns points get a drift in original
% wMu0dt(Ns+1:Nn)=wMu0dt1(Ns+1:Nn);
% dw0(1:Ns)=dw01(1:Ns);
% dw0(Ns+1:Nn)=dw02(Ns+1:Nn);
%Below we calculate different coefficients for drift and H1=Z and H2=(Z^2-1).
%For Ito Process of f=theta1*yy(t)^omega1;
[fMu0dt,dfz0,df2z0]=CalculateFunctionDriftAndVol(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,w,wnStart,Nn,dt);
df2z0(1)=0;
dfz0(1)=dfz0(2);
fMu0dt(1)=fMu0dt(2);

%Below we calculate different coefficients for drift and H1=Z and H2=(Z^2-1).
%For Ito Process of fdt=theta1*yy(t)^omega1*dt;
[fdtMu0dt1,fdtMu0dt2,dfdtz01,df2dtz0]=CalculateFdtDriftAndVol(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,w,wnStart,Nn,dt);

Mm=Nn;
for mm=1:Mm

if((mm-cN)>=1)
c1N(mm)=mm-cN;
else
c1N(mm)=1;
end
%c1N(mm) is the starting point on t+dt grid for
%point mm on the first grid.This is 5.5 SDs less than w(mm)
if((mm+cN)<=Nn-1)
c2N(mm)=mm+cN;
else
c2N(mm)=Nn-1;
end
%c2N(mm) is the ending point on t+dt grid for
%point mm on the first grid.This is 5.5 SDs more than w(mm)
%so using c1N(mm) and c2N(mm) as +-5.5 SDs, we iterate the time step
%to only relevant values and our algorithm has close to order(Nn) complexity
% rather than O(Nn^2) complexity if we had to iterate over all the
% points on the next grid.
if(mm<=Ns)  %Here for first Ns grid points, we used Zt derived from original coordinates
Zt(mm,c1N(mm):c2N(mm))=(yy(c1N(mm):c2N(mm))-yy(mm)-1*wMu0dt1(mm)+.5*sigma0^2*yy(mm).^(2*gamma-1)*dt)/dw01(mm);
pt(mm,c1N(mm):c2N(mm))=exp(-.5* Zt(mm,c1N(mm):c2N(mm)).^2)./sqrt(2*pi*sigma0.^2*dt).*dNn;
%Please not that I have used a scaling factor of sigma0^2*dt since we are integrating the transition probability
%over the target grid cell to calculate the integrated probability
%mass.So even though Zt was derived from original coordinates, the
%final integration is done in bessel coordinates to calculate the
%integrated transition probability.
end

if(mm>Ns)
Zt(mm,c1N(mm):c2N(mm))=(b2(c1N(mm):c2N(mm))-b1(mm)-1*wMu0dt2(mm))/dw02(mm);
pt(mm,c1N(mm):c2N(mm))=exp(-.5* Zt(mm,c1N(mm):c2N(mm)).^2)./sqrt(2*pi*dw02(mm).^2).*dNn;
%In the above we calculate Zt and integrated transition
%probabilities from bessel coordinates.
end

Pn(mm)=sum(pt(mm,c1N(mm):c2N(mm)));%Integrate probability mass on the next target grid that
%evolved from grid point mm.
%If true mean is not known in original coordinates just comment
%the if loop below. Currently, it uses a mean-reverting mean SV models.
if( (b1(mm)<wE0-.5*sigma0*sqrt(dt))&& (Pn(mm)>0) )

pt0(1:c2N(mm))=pt(mm,1:c2N(mm));
TrueMean=theta+(yy(mm)-theta)*exp(-kappa*(dt));
ptCorrected0 = CorrectProbAndMean(TrueMean,c1N(mm),c2N(mm),yy,pt0,Pn(mm));
pt(mm,1:c2N(mm))=ptCorrected0(1:c2N(mm));

%In this if loop, the transition probabilities are corrected so
%that sum of next stage target probabilities sum to one and
%mean equals the process mean at next stage given its evolution
%from point w(mm).
%If true mean is not known in original coordinates just comment
%this if loop.

%     P2n(mm)= sum(pt(mm,c1N(mm):c2N(mm)));
%     P2n(mm)
%    sum(yy(c1N(mm):c2N(mm)).*pt(mm,c1N(mm):c2N(mm)))

%    yy(mm)
%    theta+(yy(mm)-theta)*exp(-kappa*dt)
%    mm
%    str=input('Look at values') ;
end

dfGenerator(mm,c1N(mm):c2N(mm))=fMu0dt(mm)+dfz0(mm).*Zt(mm,c1N(mm):c2N(mm))+df2z0(mm).*(Zt(mm,c1N(mm):c2N(mm)).^2-1);
EdfGenerator(mm,c1N(mm):c2N(mm))=dfGenerator(mm,c1N(mm):c2N(mm)).*pt(mm,c1N(mm):c2N(mm));

dfdtz01H1(mm,c1N(mm):c2N(mm))=dfdtz01(mm).*Zt(mm,c1N(mm):c2N(mm));%*(t2^1.5-t1^1.5)/sqrt(3) ...
d2fdtz01H2(mm,c1N(mm):c2N(mm))=df2dtz0(mm).*(Zt(mm,c1N(mm):c2N(mm)).^2-1);%*(t2^2-t1^2)/2).*Pprev(mm).*pt(mm,c1N(mm):c2N(mm));

end
Pprev(1:Nn)=0;
Pprev(Nn0)=1;
Pnew(1:Nn)=0.0;
Efprev(1:Nn)=0;
Efprev(Nn0)=theta1*((1-gamma).*w(Nn0)).^(omega1/(1-gamma));
Efnew(1:Nn)=0;
Efdtprev(1:Nn)=0;
Efdtprev(Nn0)=0;%theta1*((1-gamma).*w(Nn0)).^(omega1/(1-gamma));
Efdtnew(1:Nn)=0;
for tt=1:Tt
t1=(tt-1)*dt;
t2=tt*dt;
dt02=(t2^2-t1^2)/2;
dt03=(t2^1.5-t1^1.5)/sqrt(3);
for mm=1:Mm-1
Pnew(c1N(mm):c2N(mm))=Pnew(c1N(mm):c2N(mm))+Pprev(mm).*pt(mm,c1N(mm):c2N(mm));
%Above is SDE process density simulation that evolves integraed probability.
%Efnew(c1N(mm):c2N(mm))=Efnew(c1N(mm):c2N(mm))+Efprev(mm).*pt(mm,c1N(mm):c2N(mm))+fGenerator(mm,c1N(mm):c2N(mm)).*Pprev(mm).*pt(mm,c1N(mm):c2N(mm));
Efnew(c1N(mm):c2N(mm))=Efnew(c1N(mm):c2N(mm))+Efprev(mm).*pt(mm,c1N(mm):c2N(mm))+EdfGenerator(mm,c1N(mm):c2N(mm)).*Pprev(mm);
%Above is simulation of expected value of first Ito-process.

Efdtnew(c1N(mm):c2N(mm))=Efdtnew(c1N(mm):c2N(mm))+Efdtprev(mm).*pt(mm,c1N(mm):c2N(mm))+ ...
((-t1*fyy(mm)+t2*fyy(c1N(mm):c2N(mm)))-fdtMu0dt1(mm)*dt02 ...
-dfdtz01H1(mm,c1N(mm):c2N(mm))*dt03 ...
-d2fdtz01H2(mm,c1N(mm):c2N(mm))*dt02).*Pprev(mm).*pt(mm,c1N(mm):c2N(mm));
%Above is simulation of expected value of second dt-ito-integral process.
%this takes longer since time values are changing on every time
%level because there is a t in Ito expansion.

end
Pprev(1:Nn)=Pnew(1:Nn);
Efprev(1:Nn)=Efnew(1:Nn);
Efdtprev(1:Nn)=Efdtnew(1:Nn);
Pnew(1:Nn)=0;
Efnew(1:Nn)=0.0;
Efdtnew(1:Nn)=0.0;
end
Pnew=Pprev;
Efnew=Efprev;
Efdtnew=Efdtprev;

%Below we calculate Z(0,1) i.e standard gaussian density associated with
%the probability mass in the grid. It assigns a Z, standard normal, value
%to each point in the grid using cumultive probability mapping

Cprob=0;
for mm=1:Mm
Cprob=Cprob+Pnew(mm);
Zb(mm)=norminv(Cprob,0,1);
end
for mm=Nn0:-1:1
if(Zb(mm)==-Inf)
Zb(mm)=Zb(mm+1);
end
end
for mm=Nn0:Mm
if(Zb(mm)==Inf)
Zb(mm)=Zb(mm-1);
end
end
wmStart=1;
wmStartFlag=1;
for mm=1:Mm-1
if((Zb(mm+1)-Zb(mm))>0  && (wmStartFlag==1))
wmStart=mm;
wmStartFlag=0;
end
end
wmEnd=Nn;
wmEndFlag=1;
for mm=Mm:-1:2
if((Zb(mm)-Zb(mm-1))>0  && (wmEndFlag==1))
wmEnd=mm;
wmEndFlag=0;
end
end

for mm=wmStart:wmEnd-1
Z(mm)=.5*Zb(mm+1)+.5*Zb(mm);
end

fy(wmStart:wmEnd-1)=Efnew(wmStart:wmEnd-1)./Pnew(wmStart:wmEnd-1);
%Convert the expected values in each cell to standard Values after
%associated with the cell after removing the
%integrated probability in the cell.Above for fy ito process.
fydt(wmStart:wmEnd-1)=Efdtnew(wmStart:wmEnd-1)./Pnew(wmStart:wmEnd-1);
%Convert the expected values in each cell to standard Values after
%associated with the cell after removing the
%integrated probability in the cell.Above for fydt ito process.

Dfy_w(wmStart:wmEnd)=0;

Dfy(1:Mm)=0;
Dfydt(1:Mm)=0;
for mm=wmStart+1:wmEnd-2
Dfy(mm) = (fy(mm + 1) - fy(mm - 1))/(Z(mm + 1) - Z(mm - 1));
Dfydt(mm) = (fydt(mm + 1) - fydt(mm - 1))/(Z(mm + 1) - Z(mm - 1));
end

ffy(1:Mm)=0;
ffydt(1:Mm)=0;
for mm = wmStart+1:wmEnd-2
ffy(mm) = (normpdf(Z(mm),0, 1))/abs(Dfy(mm));
ffydt(mm) = (normpdf(Z(mm),0, 1))/abs(Dfydt(mm));%Origianl coordinates density
end
%fx2dt
%xx2dt
for mm=wmStart+1:wmEnd-2
Dfy_w(mm) = (yy(mm + 1) - yy(mm - 1))/(Z(mm + 1) - Z(mm - 1));
%Change of variable derivative for densities
end
fy_w(1:Mm)=0;
for mm = wmStart:wmEnd-1
fy_w(mm) = (normpdf(Z(mm),0, 1))/abs(Dfy_w(mm));%Origianl coordinates density
end

%P2n
%Zt(1,:)
%pt(1,:)
sum(yy(wmStart+1:wmEnd-2).*Pnew(wmStart+1:wmEnd-2)) %Original process average from coordinates
theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average

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

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:Tt+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.
%Pre-compute the time and power exponent values in small multi-dimensional arrays
for k = 0 : (OrderA+1)
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) = dt^((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));
end
end
end
end
%Below call the program that calculates the Ito-Taylor coeffs.
YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);
%Calculate the density with monte carlo below.
rng(85109634, 'twister')
paths=200000;
YY(1:paths)=x0;  %Original process monte carlo.
YY2dt(1:paths)=0.0;
Random1(1:paths)=0;
for tt=1:Tt
t1=tt*dt;
t0=(tt-1)*dt;
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;
YY0=YY;
YYPrev=YY;
%     YY2dt(1:paths)=YY2dt(1:paths)+.5*theta1*YY(1:paths).^omega1*dt;
%     YY2dt(1:paths)=YY2dt(1:paths)+.5*dt*YYPrev(1:paths).^mp;
for k = 0 : OrderM
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;
if(l4<=5)

YY(1:paths)=YY(1:paths) + YCoeff(l1,l2,l3,l4) * ...
mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...
YY0(1:paths).^Fp(l1,l2,l3,l4) .* HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);
%Above is original coordinates monte carlo equation.
end
end
end
end
end

%     YY2dt(1:paths)=YY2dt(1:paths)+.5*theta1*YY(1:paths).^omega1*dt;

YY2dt(1:paths)=YY2dt(1:paths)+t1*YY(1:paths).^omega1-t0*YYPrev(1:paths).^omega1- ...
omega1.*(mu1*YYPrev(1:paths).^(omega1+beta1-1)+mu2*YYPrev(1:paths).^(omega1+beta2-1)).*(t1^2-t0^2)/2.0 - ...
.5*omega1*(omega1-1).*sigma0^2*YYPrev(1:paths).^(omega1+2*gamma-2).*(t1^2-t0^2)/2 - ...
omega1*sigma0*YYPrev(1:paths).^(omega1+gamma-1) .*HermiteP1(2,1:paths).*(t1^1.5-t0^1.5)/sqrt(3);

end
YY(YY<0)=0;

fYY2(1:paths)=theta1*YY(1:paths).^omega1;

sum(fYY2(:))/paths   %dt-integral average from monte carlo
sum(fy(wmStart+1:wmEnd-2).*Pnew(wmStart+1:wmEnd-2)) %dt-integral from transition
%probabilities framework.

sum(YY2dt(:))/paths   %dt-integral average from monte carlo
sum(fydt(wmStart+1:wmEnd-2).*Pnew(wmStart+1:wmEnd-2)) %dt-integral from transition
%probabilities framework.

sum(YY(:))/paths %origianl coordinates monte carlo average.
sum(yy(wmStart+1:wmEnd-2).*Pnew(wmStart+1:wmEnd-2)) %Original process average from coordinates
theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average

BinSize=.0075/1;%Please change this bin size accordingly. When data range is too small decrease it.
%When density range is too large increase it.
MaxCutOff=30;
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);
plot(yy(wmStart+1:wmEnd-2),fy_w(wmStart+1:wmEnd-2),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
str=input('red line is density of SDE from squared orthogonal incrments method, green is monte carlo.');

BinSize=.0075*2/1;%Please change this bin size accordingly. When data range is too small decrease it.

MaxCutOff=40;
[fY2Density,IndexOutfY2,IndexMaxfY2] = MakeDensityFromSimulation_Infiniti(fYY2,paths,BinSize,MaxCutOff);
plot(fy(wmStart+1:wmEnd-2),ffy(wmStart+1:wmEnd-2),'r',IndexOutfY2(1:IndexMaxfY2),fY2Density(1:IndexMaxfY2),'g');
str=input('red line is Ito transformed Process from transition probability framework , green is monte carlo.');

BinSize=.0075*2/1;%Please change this bin size accordingly. When data range is too small decrease it.
%When density range is too large increase it.
MaxCutOff=40;
[Y2dtDensity,IndexOutY2dt,IndexMaxY2dt] = MakeDensityFromSimulation_Infiniti(YY2dt,paths,BinSize,MaxCutOff);
plot(fydt(wmStart+1:wmEnd-2),ffydt(wmStart+1:wmEnd-2),'r',IndexOutY2dt(1:IndexMaxY2dt),Y2dtDensity(1:IndexMaxY2dt),'g');
str=input('red line is dt-integral from transition probability framework , green is monte carlo.');

end

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 are the dependencies functions of the above program.

function [ptCorrected] = CorrectProbAndMean(TrueMean,c1N,c2N,yy,pt,Pn)
%  u0=theta+(yy(mm)-theta)*exp(-kappa*(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));
%
u0=TrueMean;
Yn2Pn=sum(yy(c1N:c2N).*yy(c1N:c2N).*pt(c1N:c2N));
YnPn=sum(yy(c1N:c2N).*pt(c1N:c2N));
Pn0=Pn;%%Sum(ZProb(1:Nn))
Y0=(Yn2Pn-u0*YnPn)/(YnPn-u0*Pn0);
%
Y0Pn=Y0*Pn0;
W0=1/(YnPn-Y0Pn);
%      W0
%      Y0
%      yy(c1N)
%      yy(c2N)
if((Y0>=yy(c1N))&&(Y0<yy(c2N)))
if(W0<0)
Y0=yy(c1N)-1;
else
Y0=yy(c2N)+.01;
end
end
Y0Pn=Y0*Pn0;
W0=1/(YnPn-Y0Pn);
ptCorrected(c1N:c2N)=pt(c1N:c2N).*(yy(c1N:c2N)-Y0)*W0;

%sum( ptCorrected(c1N:c2N));
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
function [YMu0dt,dw0]=CalculateDriftAndVolOriginalCoordinates(mu1,mu2,beta1,beta2,gamma,sigma0,w,wnStart,Nn,dt)

dt0=dt;
dtsqrt=sqrt(dt);
dtonehalf=dt^1.5;
dt2=dt^2/2;
dt3=dt^3/6.0;

theta1=mu1;
theta2=mu2;
theta3=0;%.5*sigma0^2*(-gamma);
omega1=beta1;
omega2=beta2;
omega3=0;
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));

DD(wnStart:Nn)=mu1*yy(wnStart:Nn).^beta1+mu2*yy(wnStart:Nn).^beta2;
dDD(wnStart:Nn)=mu1*beta1*yy(wnStart:Nn).^(beta1-1)+mu2*beta2*yy(wnStart:Nn).^(beta2-1);

f(wnStart:Nn)=theta1*yy(wnStart:Nn).^(omega1)+theta2*yy(wnStart:Nn).^(omega2)+ ...
theta3*yy(wnStart:Nn).^(omega3);

df(wnStart:Nn)=omega1*theta1*yy(wnStart:Nn).^(omega1-1)+omega2*theta2*yy(wnStart:Nn).^(omega2-1)+ ...
omega3*theta3*yy(wnStart:Nn).^(omega3-1);

d2f(wnStart:Nn)=theta1*omega1*(omega1-1).*yy(wnStart:Nn).^(omega1-2)+ ...
theta2*omega2*(omega2-1).*yy(wnStart:Nn).^(omega2-2)+ ...
theta3*omega3*(omega3-1).*yy(wnStart:Nn).^(omega3-2);
QQ(wnStart:Nn)=sigma0^2*yy(wnStart:Nn).^(2*gamma);
ddGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*yy(wnStart:Nn).^(omega1+beta1-2)+ ...
mu1*omega2*theta2*(omega2+beta1-1).*yy(wnStart:Nn).^(omega2+beta1-2)+ ...
mu1*omega3*theta3*(omega3+beta1-1).*yy(wnStart:Nn).^(omega3+beta1-2)+ ...
mu2*omega1*theta1*(omega1+beta2-1).*yy(wnStart:Nn).^(omega1+beta2-2)+ ...
mu2*omega2*theta2*(omega2+beta2-1).*yy(wnStart:Nn).^(omega2+beta2-2)+ ...
mu2*omega3*theta3*(omega3+beta2-1).*yy(wnStart:Nn).^(omega3+beta2-2);

d_dfDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*yy(wnStart:Nn).^(omega1+beta1-2)+ ...
mu2*omega1*theta1*(omega1+beta2-1).*yy(wnStart:Nn).^(omega1+beta2-2)+ ...
mu1*omega2*theta2*(omega2+beta1-1).*yy(wnStart:Nn).^(omega2+beta1-2)+ ...
mu2*omega2*theta2*(omega2+beta2-1).*yy(wnStart:Nn).^(omega2+beta2-2)+ ...
mu1*omega3*theta3*(omega3+beta1-1).*yy(wnStart:Nn).^(omega3+beta1-2)+ ...
mu2*omega3*theta3*(omega3+beta2-1).*yy(wnStart:Nn).^(omega3+beta2-2);

d2_dfDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*(omega1+beta1-2).*yy(wnStart:Nn).^(omega1+beta1-3)+ ...
mu1*omega2*theta2*(omega2+beta1-1).*(omega2+beta1-2).*yy(wnStart:Nn).^(omega2+beta1-3)+ ...
mu1*omega3*theta3*(omega3+beta1-1).*(omega3+beta1-2).*yy(wnStart:Nn).^(omega3+beta1-3)+ ...
mu2*omega1*theta1*(omega1+beta2-1).*(omega1+beta2-2).*yy(wnStart:Nn).^(omega1+beta2-3)+ ...
mu2*omega2*theta2*(omega2+beta2-1).*(omega2+beta2-2).*yy(wnStart:Nn).^(omega2+beta2-3)+ ...
mu2*omega3*theta3*(omega3+beta2-1).*(omega3+beta2-2).*yy(wnStart:Nn).^(omega3+beta2-3);
d_d2fQQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*yy(wnStart:Nn).^(omega1+2*gamma-3)+ ...
sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*yy(wnStart:Nn).^(omega2+2*gamma-3)+ ...
sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*yy(wnStart:Nn).^(omega3+2*gamma-3);
d2_d2fQQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*(omega1+2*gamma-3).*yy(wnStart:Nn).^(omega1+2*gamma-4)+ ...
sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*(omega2+2*gamma-3).*yy(wnStart:Nn).^(omega2+2*gamma-4)+ ...
sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*(omega3+2*gamma-3).*yy(wnStart:Nn).^(omega3+2*gamma-4);

VV(wnStart:Nn)=sigma0.*yy(wnStart:Nn).^gamma;

dVV(wnStart:Nn)=sigma0.*gamma*yy(wnStart:Nn).^(gamma-1);
d2VV(wnStart:Nn)=sigma0.*gamma*(gamma-1)*yy(wnStart:Nn).^(gamma-2);

YMu0dt(wnStart:Nn)=f(wnStart:Nn)*dt+ ...
df(wnStart:Nn).*DD(wnStart:Nn)*dt2 + ...
.5*d2f(wnStart:Nn).*QQ(wnStart:Nn)*dt2+ ...
d_dfDD(wnStart:Nn).*DD(wnStart:Nn)*dt3 + ...
.5*d_d2fQQ(wnStart:Nn).*DD(wnStart:Nn)*dt3 + ...
.5*d2_dfDD(wnStart:Nn).*QQ(wnStart:Nn)*dt3 + ...
.25*d2_d2fQQ(wnStart:Nn).*QQ(wnStart:Nn)*dt3;

dw0(1:Nn)=VV(wnStart:Nn).*dtsqrt+ ...
dVV(wnStart:Nn).*DD(wnStart:Nn)*1/sqrt(3).*dt^1.5+ ...
.5*d2VV(wnStart:Nn).*QQ(wnStart:Nn)*1/sqrt(3).*dt^1.5+ ...
d_dfDD(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3)).*dt^1.5 + ...
.5*d_d2fQQ(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3)).*dt^1.5;

%dVV(wnStart:Nn).*VV(wnStart:Nn).*dtsqrt+ ...

%dw0(1:Nn)=sigma0.*sqrt(dt0);%+dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*dt0^1.5;

end

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

function [fMu0dt1,fMu0dt2,dfz01,df2z0]=CalculateFdtDriftAndVol(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,w,wnStart,Nn,dt)
dtsqrt=sqrt(dt);
dtonehalf=dt^1.5;
dt2=dt^2/2;
dt3=dt^3/6.0;
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
f(wnStart:Nn)=theta1*yy(wnStart:Nn).^(omega1);
df(wnStart:Nn)=omega1*theta1*yy(wnStart:Nn).^(omega1-1);
DD(wnStart:Nn)=mu1*yy(wnStart:Nn).^beta1+mu2*yy(wnStart:Nn).^beta2;
d2f(wnStart:Nn)=theta1*omega1*(omega1-1).*yy(wnStart:Nn).^(omega1-2);

QQ(wnStart:Nn)=sigma0^2*yy(wnStart:Nn).^(2*gamma);
d_dfDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*yy(wnStart:Nn).^(omega1+beta1-2)+ ...
mu2*omega1*theta1*(omega1+beta2-1).*yy(wnStart:Nn).^(omega1+beta2-2);
d2_dfDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*(omega1+beta1-2).*yy(wnStart:Nn).^(omega1+beta1-3)+ ...
mu2*omega1*theta1*(omega1+beta2-1).*(omega1+beta2-2).*yy(wnStart:Nn).^(omega1+beta2-3);

d_d2fQQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*yy(wnStart:Nn).^(omega1+2*gamma-3);

d2_d2fQQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*(omega1+2*gamma-3).*yy(wnStart:Nn).^(omega1+2*gamma-4);

VV(wnStart:Nn)=sigma0.*yy(wnStart:Nn).^gamma;

d_dfVV(wnStart:Nn)= ...
sigma0.*omega1.*theta1.*(omega1+gamma-1).*yy(wnStart:Nn).^(omega1+gamma-2);

d2_dfVV(wnStart:Nn)= ...
sigma0.*omega1.*theta1.*(omega1+gamma-1).*(omega1+gamma-2).*yy(wnStart:Nn).^(omega1+gamma-3);

fMu0dt1(wnStart:Nn)=df(wnStart:Nn).*DD(wnStart:Nn) + ...
.5*d2f(wnStart:Nn).*QQ(wnStart:Nn);%dt+ ...
fMu0dt2(wnStart:Nn)= d_dfDD(wnStart:Nn).*DD(wnStart:Nn) + ...
.5*d_d2fQQ(wnStart:Nn).*DD(wnStart:Nn) + ...
.5*d2_dfDD(wnStart:Nn).*QQ(wnStart:Nn) + ...
.25*d2_d2fQQ(wnStart:Nn).*QQ(wnStart:Nn);%dt2

dfz01(1:Nn)=df(wnStart:Nn).*VV(wnStart:Nn);%.*dtsqrt+ ...
%dfz02(1:Nn)=d_dfVV(wnStart:Nn).*DD(wnStart:Nn)+ ...
%    .5*d2_dfVV(wnStart:Nn).*QQ(wnStart:Nn)+ ...
%    d_dfDD(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3)) + ...
%    .5*d_d2fQQ(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3));%.*dt^1.5/sqrt(3);

df2z0(1:Nn)=d_dfVV(wnStart:Nn).*VV(wnStart:Nn);%.*dt/2;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [wMu0dt0,dw0]=CalculateDriftAndVolSmall(mu1,mu2,beta1,beta2,gamma,sigma0,w,wnStart,Nn,dt)

dt0=dt;

theta1=mu1;
theta2=mu2;
theta3=.5*sigma0^2*(-gamma);
omega1=beta1-gamma;
omega2=beta2-gamma;
omega3=gamma-1;
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));

DD(wnStart:Nn)=mu1*yy(wnStart:Nn).^beta1+mu2*yy(wnStart:Nn).^beta2;
dDD(wnStart:Nn)=mu1*beta1*yy(wnStart:Nn).^(beta1-1)+mu2*beta2*yy(wnStart:Nn).^(beta2-1);
yyf(wnStart:Nn)=yy(wnStart:Nn)+DD(wnStart:Nn).*dt0+dDD(wnStart:Nn).*DD(wnStart:Nn)*dt0^2/2;

wf(wnStart:Nn)=yyf(wnStart:Nn).^(1-gamma)/(1-gamma);

wMu0dt0(wnStart:Nn)=wf(wnStart:Nn)-w(wnStart:Nn);

dGG(wnStart:Nn)=omega1*theta1*yy(wnStart:Nn).^(omega1-1)+omega2*theta2*yy(wnStart:Nn).^(omega2-1)+ ...
omega3*theta3*yy(wnStart:Nn).^(omega3-1);

VV(wnStart:Nn)=sigma0.*yy(wnStart:Nn).^gamma;

dw0(1:Nn)=sigma0.*sqrt(dt0);%+dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*dt0^1.5;

end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [fMu0dt,dfz0,df2z0]=CalculateFunctionDriftAndVol(theta1,omega1,mu1,mu2,beta1,beta2,gamma,sigma0,w,wnStart,Nn,dt)
dtsqrt=sqrt(dt);
dtonehalf=dt^1.5;
dt2=dt^2/2;
dt3=dt^3/6.0;
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
f(wnStart:Nn)=theta1*yy(wnStart:Nn).^(omega1);
df(wnStart:Nn)=omega1*theta1*yy(wnStart:Nn).^(omega1-1);
DD(wnStart:Nn)=mu1*yy(wnStart:Nn).^beta1+mu2*yy(wnStart:Nn).^beta2;
d2f(wnStart:Nn)=theta1*omega1*(omega1-1).*yy(wnStart:Nn).^(omega1-2);

QQ(wnStart:Nn)=sigma0^2*yy(wnStart:Nn).^(2*gamma);
d_dfDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*yy(wnStart:Nn).^(omega1+beta1-2)+ ...
mu2*omega1*theta1*(omega1+beta2-1).*yy(wnStart:Nn).^(omega1+beta2-2);
d2_dfDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*(omega1+beta1-2).*yy(wnStart:Nn).^(omega1+beta1-3)+ ...
mu2*omega1*theta1*(omega1+beta2-1).*(omega1+beta2-2).*yy(wnStart:Nn).^(omega1+beta2-3);

d_d2fQQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*yy(wnStart:Nn).^(omega1+2*gamma-3);

d2_d2fQQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*(omega1+2*gamma-3).*yy(wnStart:Nn).^(omega1+2*gamma-4);

VV(wnStart:Nn)=sigma0.*yy(wnStart:Nn).^gamma;

d_dfVV(wnStart:Nn)= ...
sigma0.*omega1.*theta1.*(omega1+gamma-1).*yy(wnStart:Nn).^(omega1+gamma-2);

d2_dfVV(wnStart:Nn)= ...
sigma0.*omega1.*theta1.*(omega1+gamma-1).*(omega1+gamma-2).*yy(wnStart:Nn).^(omega1+gamma-3);

fMu0dt(wnStart:Nn)=df(wnStart:Nn).*DD(wnStart:Nn).*dt + ...
.5*d2f(wnStart:Nn).*QQ(wnStart:Nn).*dt+ ...
d_dfDD(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...
.5*d_d2fQQ(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...
.5*d2_dfDD(wnStart:Nn).*QQ(wnStart:Nn).*dt2 + ...
.25*d2_d2fQQ(wnStart:Nn).*QQ(wnStart:Nn).*dt2;

dfz0(1:Nn)=df(wnStart:Nn).*VV(wnStart:Nn).*dtsqrt+ ...
d_dfVV(wnStart:Nn).*DD(wnStart:Nn)*1/sqrt(3).*dt^1.5+ ...
.5*d2_dfVV(wnStart:Nn).*QQ(wnStart:Nn)*1/sqrt(3).*dt^1.5+ ...
d_dfDD(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3)).*dt^1.5 + ...
.5*d_d2fQQ(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3)).*dt^1.5;

df2z0(1:Nn)=d_dfVV(wnStart:Nn).*VV(wnStart:Nn).*dt/2;

end

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

function [wMu0dt,dw0]=CalculateDriftAndVol(mu1,mu2,beta1,beta2,gamma,sigma0,w,wnStart,Nn,dt)
dtsqrt=sqrt(dt);
dtonehalf=dt^1.5;
dt2=dt^2/2;
dt3=dt^3/6.0;
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
theta1=mu1;
theta2=mu2;
theta3=.5*sigma0^2*(-gamma);
omega1=beta1-gamma;
omega2=beta2-gamma;
omega3=gamma-1;
GG(wnStart:Nn)=theta1*yy(wnStart:Nn).^omega1+theta2*yy(wnStart:Nn).^omega2+theta3*yy(wnStart:Nn).^omega3;
dGG(wnStart:Nn)=omega1*theta1*yy(wnStart:Nn).^(omega1-1)+omega2*theta2*yy(wnStart:Nn).^(omega2-1)+ ...
omega3*theta3*yy(wnStart:Nn).^(omega3-1);
DD(wnStart:Nn)=mu1*yy(wnStart:Nn).^beta1+mu2*yy(wnStart:Nn).^beta2;
dGG2(wnStart:Nn)=theta1*omega1*(omega1-1).*yy(wnStart:Nn).^(omega1-2)+ ...
theta2*omega2*(omega2-1).*yy(wnStart:Nn).^(omega2-2)+ ...
theta3*omega3*(omega3-1).*yy(wnStart:Nn).^(omega3-2);
QQ(wnStart:Nn)=sigma0^2*yy(wnStart:Nn).^(2*gamma);
ddGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*yy(wnStart:Nn).^(omega1+beta1-2)+ ...
mu1*omega2*theta2*(omega2+beta1-1).*yy(wnStart:Nn).^(omega2+beta1-2)+ ...
mu1*omega3*theta3*(omega3+beta1-1).*yy(wnStart:Nn).^(omega3+beta1-2)+ ...
mu2*omega1*theta1*(omega1+beta2-1).*yy(wnStart:Nn).^(omega1+beta2-2)+ ...
mu2*omega2*theta2*(omega2+beta2-1).*yy(wnStart:Nn).^(omega2+beta2-2)+ ...
mu2*omega3*theta3*(omega3+beta2-1).*yy(wnStart:Nn).^(omega3+beta2-2);
d2dGGDD(wnStart:Nn)=mu1*omega1*theta1*(omega1+beta1-1).*(omega1+beta1-2).*yy(wnStart:Nn).^(omega1+beta1-3)+ ...
mu1*omega2*theta2*(omega2+beta1-1).*(omega2+beta1-2).*yy(wnStart:Nn).^(omega2+beta1-3)+ ...
mu1*omega3*theta3*(omega3+beta1-1).*(omega3+beta1-2).*yy(wnStart:Nn).^(omega3+beta1-3)+ ...
mu2*omega1*theta1*(omega1+beta2-1).*(omega1+beta2-2).*yy(wnStart:Nn).^(omega1+beta2-3)+ ...
mu2*omega2*theta2*(omega2+beta2-1).*(omega2+beta2-2).*yy(wnStart:Nn).^(omega2+beta2-3)+ ...
mu2*omega3*theta3*(omega3+beta2-1).*(omega3+beta2-2).*yy(wnStart:Nn).^(omega3+beta2-3);
ddGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*yy(wnStart:Nn).^(omega1+2*gamma-3)+ ...
sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*yy(wnStart:Nn).^(omega2+2*gamma-3)+ ...
sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*yy(wnStart:Nn).^(omega3+2*gamma-3);
d2dGG2QQ(wnStart:Nn)=theta1*sigma0^2.*omega1*(omega1-1).*(omega1+2*gamma-2).*(omega1+2*gamma-3).*yy(wnStart:Nn).^(omega1+2*gamma-4)+ ...
sigma0^2*theta2*omega2*(omega2-1).*(omega2+2*gamma-2).*(omega2+2*gamma-3).*yy(wnStart:Nn).^(omega2+2*gamma-4)+ ...
sigma0^2*theta3*omega3*(omega3-1).*(omega3+2*gamma-2).*(omega3+2*gamma-3).*yy(wnStart:Nn).^(omega3+2*gamma-4);

VV(wnStart:Nn)=sigma0.*yy(wnStart:Nn).^gamma;

wMu0dt(wnStart:Nn)=GG(wnStart:Nn)*dt + dGG(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...
.5*dGG2(wnStart:Nn).*QQ(wnStart:Nn).*dt2+ ...
ddGGDD(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
.5*ddGG2QQ(wnStart:Nn).*DD(wnStart:Nn).*dt3 + ...
.5*d2dGGDD(wnStart:Nn).*QQ(wnStart:Nn).*dt3 + ...
.25*d2dGG2QQ(wnStart:Nn).*QQ(wnStart:Nn).*dt3;

dw0(1:Nn)=sigma0.*dtsqrt;%+dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*dtonehalf;

end