SERVING THE QUANTITATIVE FINANCE COMMUNITY

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

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

In the above post my definition of hermite polynomial might be slightly non-standard
and here are first few hermites according to my definition
H1(z(t))=z(t)
H2(z(t))=z(t)^2-t.
H3(Z(t))=z(t)^3-3*t*z(t)
and so on so we have t dependent coefficients in different terms based on variance. These are the kinds of hermite polynomials in integrals that arise in solution of Ito-Taylor stochastic integrals.
Even though I made a mistake when I earlier said that these stochastic integrals commute, they do take the same order of hermite polynomials as earlier but their coefficients change based on the order of stochastic integrals.

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

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

I am posting a new version of the old programs correcting the error. I have run the new program on density of square of a lognormal SDE with drift. And it can then be compared to analytical solution for the density of squared lognormal. Here is the new code with Ito-Taylor.
Clear[x,t,x0,Y,P1,P2,l1,l2,Z,mu,sigma,mu0,sigma0,alpha,beta,gamma,a0,k,k1,m,L1,L2,YAns,H,Hc,Hz,X1,fX1,DfX1,Nn,n,HCoeff,ZCoeff,DZCoeff,zz];
P1=24;
P2=24;
L1=0;
L2=0;
Array[Y,{P1+L1+1,P1+L1+1,P2+L2+1},{0,0,0}];
For[k=0,k<=P1+L1,k++,(For[m=0,m<=P1+L1,m++,(For[n=0,n<=P2+L2,n++,(Y[k,m,n]=0);]);]);];
a0=1;
mu0=.25;
sigma0=.75;
alpha=2;
gamma=1.0;
beta=1;
Y[L1,L1,L2]=a0;
For[k=0,k<P1,k++,(For[m=0,m<=If[k<=P2,k,P2],m++,(l0=L1+k-m;l3=L2+m;
For[n=0,n<=l0-L1,n++,(l1=l0-n;l2=L1+n;
Y[l1+1,l2,l3]=Y[l1+1,l2,l3]+SetPrecision[(1.0/(l1+l2+1)*(1.0-Sqrt[l3]/Sqrt[2*l1+2*l2+l3+2]))*Y[l1,l2,l3]*(alpha+(l1-L1)*beta+(l2-L1)*2*gamma+(l3-L2)*gamma-(l1-L1)-2*(l2-L1)-(l3-L2))*mu0,50];
Y[l1,l2+1,l3]=Y[l1,l2+1,l3]+SetPrecision[(1.0/(l1+l2+1)*(1.0-Sqrt[l3]/Sqrt[2*l1+2*l2+l3+2]))*Y[l1,l2,l3]*.5*(alpha+(l1-L1)*beta+(l2-L1)*2*gamma+(l3-L2)*gamma-(l1-L1)-2*(l2-L1)-(l3-L2))*(alpha+(l1-L1)*beta+(l2-L1)*2*gamma+(l3-L2)*gamma-(l1-L1)-2*(l2-L1)-(l3-L2)-1)*sigma0^2,50];
If[l3<L2+P2,(Y[l1,l2,l3+1]=Y[l1,l2,l3+1]+SetPrecision[(1.0/Sqrt[2*l1+2*l2+l3+1]/Sqrt[l3+1])*Y[l1,l2,l3]*(alpha+(l1-L1)*beta+(l2-L1)*2*gamma+(l3-L2)*gamma-(l1-L1)-2*(l2-L1)-(l3-L2))*sigma0,50]),0];);];)];);];
Array[HCoeff,P2+L2+1,0];
For[k=0,k<P2+L2+1,k++,HCoeff[k]=0;]
x=1;
t=2;
For[k=0,k<=P1,k++,(For[m=0,m<=If[k<=P2,k,P2],m++,(l0=L1+k-m;l3=L2+m;
For[n=0,n<=l0-L1,n++,(l1=l0-n;l2=L1+n;
HCoeff[l3]=HCoeff[l3]+SetPrecision[Y[l1,l2,l3]*x^(alpha+(l1-L1)*beta+(l2-L1)*2*gamma+(l3-L2)*gamma-(l1-L1)-2*(l2-L1)-(l3-L2))*(Sqrt[t^(2*l1+2*l2-2*L1+l3)]),50];)];)];);];
Array[Hz,{P2+L2+1,P2+L2+1},{0,0}];
For[k=0,k<P2+L2+1,k++,For[m=0,m<P2+L2+1,m++,Hz[k,m]=0;]];
Hz[0,0]=1.0;
For[k=0,k<P2+L2,k++,(For[m=0,m<=k,m++,(Hz[k+1,m+1]=SetPrecision[Hz[k,m],50];
If[m>=1,Hz[k+1,m-1]=SetPrecision[Hz[k+1,m-1]-m*Hz[k,m],50],0];)])];
Array[ZCoeff,P2+L2+1,0];
For[k=0,k<P2+L2+1,k++,ZCoeff[k]=0];
For[k=0,k<P2+L2+1,k++,(For[m=0,m<=k,m++,ZCoeff[m]=ZCoeff[m]+SetPrecision[Hz[k,m]*HCoeff[k],50]])];
Array[DZCoeff,P2+L2,0];
For[k=0,k<P2+L2,k++,DZCoeff[k]=0];
For[k=0,k<P2+L2,k++,DZCoeff[k]=SetPrecision[ZCoeff[k+1]*(k+1),50]];
Nn=81;
Array[X1,Nn,0];
Array[DfX1,Nn,0];
Array[fX1,Nn,0];
Array[zz,Nn,0];
For[n=0,n<Nn,n++,(X1[n]=0;DfX1[n]=0;fX1[n]=0;zz[n]=0;)];
For[n=0,n<Nn,n++,(zz[n]=n*.1-4;X1[n]=0;DfX1[n]=0;fX1[n]=0;
X1[n]=ZCoeff[0];DfX1[n]=DZCoeff[0];For[k=1,k<P2+L2+1,k++,(X1[n]=X1[n]+SetPrecision[ZCoeff[k]*zz[n]^k,50];If[k<P2+L2,DfX1[n]=DfX1[n]+SetPrecision[DZCoeff[k]*zz[n]^k,50],0];)];
fX1[n]=SetPrecision[PDF[NormalDistribution[0,1],zz[n]]/Abs[DfX1[n]],50];)];
data1=Table[{X1[n],fX1[n]},{n,0,58}];
p1=ListLinePlot[data1]
And here is the code for analytical density of square of lognormal SDE with drift.
For[ n = 0, n < 900, n++, ( Z[n] = n*.01 + .0001;
Logf[n] = .50/Sqrt[2*Pi]/Z[n]/(sigma0*Sqrt[t])* Exp[-.5*(Log[Sqrt[Z[n]]] - Log[Sqrt[x]] - (mu0 - .5*sigma0^2)*t)^2/(sigma0^2*t)];)];
data2 = Table[{Z[n], Logf[n]}, {n, 0, 800}];
p2 = ListLinePlot[data2, PlotRange -> All]
And then single line of code for overlaying both densities in a single graph
Show[p1, p2]
I hope that this new program will be helpful.

Here is the graph comparison image. The parameters are mu=.25, sigma=.75, X0=1,t=2 for geometric Brownian motion/lognormal SDE. And graphs show the density of      squared X(t) that is  X(t)^2.  The following graph shows the overlaid comparison for analytical solution with Ito-Taylor expansion. Both images are indistinguishable. I have removed a lot of the area from the tail where probability was miniscule to focus on the main area.

Here is graph comparison for the lognormal SDE density with drift but that is not squared. The parameters are mu=.25, sigma=.75, X0=1,t=2. It is again overlaid graph comparison with the same parameters as described earlier but it is the SDE variable and not its square. Here is the overlaid comparison graph. I have again tried to focus on the main graph body and not on the tails.

Here is another overlaid graph comparison of ito-taylor expansion and analytical formula for five year density of lognormal SDE with slightly decreased volatility. Here mu=.25, sigma=.65,t=5, X0=1.

Last edited by Amin on May 11th, 2018, 8:52 pm, edited 6 times in total.

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

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

Here I have latexed the update formulas for stochastic integrals

$\int_0^t s^m H_{n}(z(s)) ds = \frac {t^{m+1} H_{n}(z(t))}{(m+1)} (1 - \frac {\sqrt{n}}{\sqrt{(2m+n+2)}})$

$\int_0^t s^m H_{n}(z(s)) dz(s)= \frac {t^m H_{n+1}(z(t))}{ \sqrt{(n+1)} \sqrt{(2m+n+1)}}$

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

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

Doing stress test with lognormal SDE with high drift to make sure that approximation matches with closed form formula. here mu=.75, sigma=.75, x0=1, t=5. I have used 24 hermite polynomials.
And both analytical and Ito-Taylor graphs overlaid are shown below. I have taken out extreme tail so that focus is on the main body of graph where we can make a visual comparison.

The purpose here is not to use the method for lognormal SDE but to make sure that the method works with lognormal SDE where a closed form solution is known and to make sure that calculation of stochastic integrals is sound. Method has largest domain of validity with lognormal type SDEs where diffusion term is lognormal like. And it has very small domain of validity with CEV exponent=.5 or diffusion terms like in heston. And a crossover for domain of validity is diffusion terms with powers .75 or greater where domain of validity is reasonable and can easily be several years. And for lognormal like diffusion terms, I believe you can easily do Ito-Taylor expansion for more than ten or twenty years or even more. And for some SDEs of interest that take mean reversion with lognormal like diffusion terms, I am sure the method is very interesting. Even if domain of validity is small, the method can be used with other techniques to easily simulate(without monte carlo) the SDEs or their various stochastic integrals out to practically any time horizon.

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

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

For friends who did not follow the earlier posts where I gave some basic details. here is the general SDE whose density is found by Ito-Taylor.
$dX(t)=\mu X(t)^{\beta} dt + \sigma X(t)^{\gamma} dz(t)$
and relating to parameters used in mathematica
mu0=$\mu$
sigma0=$\sigma$
and density is found for $X^{\alpha}$
For example for lognormal density calculations beta=1, gamma=1 and also alpha=1.
For calculation of density of X^2 where dX is given by lognormal SDE, beta=1, gamma=1, and alpha=2.

and you can specify the order of dz and dt integrals by changing the values of L1 and L2.
If L1=0 and L2=0, it means that density would be calculated for $X^{\alpha}$
If L1=1 and L2=0, it means that order of dt integral is one and the density would be calculated for
$\int_0^t X(s)^{\alpha} ds$
if L1=0 and L2=1, it means that density would be calculated for
$\int_0^t X(s)^{\alpha} dz(s)$
At this stage I would not use the program for higher mixed integrals since I did not redo this part for mixed integrals when they do not commute and it would require slight tweaking.
and in the mid of first code snippet, I specify the initial conditions for x, the starting point of diffusion and time horizon t. They are shown in mathematica code as
x=1;
t=2;
And you can play around with these initial values and time horizon.

And Sorry, there was a minor error in my analytic density of $X^2$ where X is the lognormal density variable. This is regarding second code snippet in post 692. here I have to use log[x] instead of log[Sqrt[x]]. But it was working fine for particular initial condition since the initial condition was x=1. But results hold seamlessly for any initial value after this minor change. This is error in my analytic formula for density of X^2 of lognormal SDE and Ito-Taylor part was fine as earlier posted. I am reposting second code snippet of post 692 after this correction.
For[ n = 0, n < 900, n++, ( Z[n] = n*.01 + .0001; Logf[n] = .50/Sqrt[2*Pi]/Z[n]/(sigma0*Sqrt[t])* Exp[-.5*(Log[Sqrt[Z[n]]] - Log[x] - (mu0 - .5*sigma0^2)*t)^2/(sigma0^2*t)];)];
data2 = Table[{Z[n], Logf[n]}, {n, 0, 800}];
p2 = ListLinePlot[data2, PlotRange -> All]

And for those friends who would want to compare the lognormal density itself. Here is the code snippet for that. I did not share this previously.
For[ n = 0, n < 1000, n++, ( Z[n] = n*.01 + .0001; Logf[n] = 1/Sqrt[2*Pi]/Z[n]/(sigma0*Sqrt[t])* Exp[-.5*(Log[Z[n]] - Log[x] - (mu0 - .5*sigma0^2)*t)^2/(sigma0^2*t)];)];
data2 = Table[{Z[n], Logf[n]}, {n, 0, 500}];
p2 = ListLinePlot[data2, PlotRange -> All]

Please make sure to change the ranges for density calculations and showing graph according to the true density range and visuals.

Giving slight more detail about working of the code. I have calculated coefficients of Ito-Taylor in an array Y(l1,l2,l3). l1 corresponds to order of dt-integral with respect to drift term. l2 shows the order of dt-integral with respect to quadratic variations. l3 shows order of dz-integrals. And there are all sorts of combinations of l1,l2,l3 carrying Ito-Taylor coefficients of the respective iterated stochastic integrals. while calculating these coefficients based on differentiation, I also multiply them with weights for hermite polynomials required for hermite polynomial update. so basically updating from Y(l1,l2,l3) to Y(l1+1,l2,l3), I multiply them with
(1.0/(l1+l2+1)*(1.0-Sqrt[l3]/Sqrt[2*l1+2*l2+l3+2])) which are the coefficients of first hermite dt-integral update in post 693. here order of dt-integrals is l1+l2 which is the sum of order associated with drift and quadratic variation dt-integrals. Update to Y(l1,l2+1,l3) from Y(l1,l2,l3) is exactly the same in terms of l1,l2 and l3. And then Update to Y(l1,l2,l3+1) from Y(l1,l2,l3) requires dz-integral update coefficients which are given by second formula in 693 and coded as  (1.0/Sqrt[2*l1+2*l2+l3+1]/Sqrt[l3+1]) .

And also P1 and P2 are order of expansion of dt and dz terms. For lognormal I freely used P1=P2=24. Though lognormal-like diffusions can be easily expanded to arbitrary higher hermite polynomials without any problems, I have seen that in case of heston-like diffusions where gamma=.5 become unstable as the order of expansions increases. In order to play with expansion order, you have to change P1 and P2.

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

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

I have two questions so as to have clarity about my contribution to literature.

a) Has somebody systematically solved these stochastic integrals that I wrote in post 693  before? I was reading some recent research papers and I noticed that authors did first few integrals with hand with a lot of effort. Would like to know if there is any reference about these integrals. Can somebody help with the literature review here? Of course, would be good to claim ownership if I have done them first so as to obviate the possibility of others claiming the same thing. But surely some expert literature review would be helpful. I believe any stochastic integral of Brownian motion can be solved with the equations in post 693.

b) Also I mostly see applications of Ito-Taylor expansions to monte carlo simulations of SDEs in almost all of the research papers. I have presented a framework in this thread to analytically calculate the density of SDEs and their stochastic integrals  and relate these derived densities to density of standard normal random variable. Have some other people done that before?

I also see some contribution in writing a simple algorithm to calculate the Ito-Taylor coefficients in an efficient manner without resorting to symbolic computing but that is probably more basic. But I believe there is enough original material now to put together a very good research paper after throwing in some interesting , useful and practical examples that are based on analytic solution of densities of stochastic integrals of SDEs and their applications to mathematical finance.

And I made a lot of errors on the way. Thanks to all friends for pardoning and overlooking the mistakes.

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

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

Here I have latexed the update formulas for stochastic integrals

$\int_0^t s^m H_{n}(z(s)) ds = \frac {t^{m+1} H_{n}(z(t))}{(m+1)} (1 - \frac {\sqrt{n}}{\sqrt{(2m+n+2)}})$

$\int_0^t s^m H_{n}(z(s)) dz(s)= \frac {t^m H_{n+1}(z(t))}{ \sqrt{(n+1)} \sqrt{(2m+n+1)}}$
As I suggested, we can calculate most stochastic integrals using this update. Let us suppose, we want to calculate
$\int_0^t z(s)^4 ds$  Equation(1)
Here is the way to go. We have to do it in terms of hermite polynomials.
fourth hermite polynomial representation of Brownian motion goes like
$H_4(z(t)) = z(t)^4 -6 t z(t)^2 +3 t^2$  Eq(2)
and second hermite polynomial is
$H_2(z(t)) = z(t)^2 - t$  Eq(3)
and we know from the quoted post
$\int_0^t H_4(z(s)) ds = t H_4(z(t)) (1- \frac{2}{\sqrt{6}} )$  Eq(4)
so from Eq(2) and Eq(4)
$\int_0^t z(s)^4 ds = t H_4(z(t)) (1 - \frac{2}{\sqrt{6}} ) + 6 \int_0^t s z(s)^2 ds - 3 \int_0^t s^2 ds$ Eq(5)
In order to solve the second term, we take the second hermite polynomial and use the formula in quotes to calculate that
$\int_0^t s H_2(z(s)) ds = \frac{t^2 H_2(z(t))}{2} (1 - \frac{\sqrt{2}}{\sqrt{6}} )$  Eq(6)
so
$\int_0^t s z(s)^2 ds = \frac{t^2 H_2(z(t))}{2} ( 1 - \frac{\sqrt{2}}{\sqrt{6}} ) + \int_0^t s^2 ds$  Eq(7)
and we substitute Eq(7) in Eq(5) to get
$\int_0^t z(s)^4 ds = t H_4(z(t)) ( 1 - \frac{2}{\sqrt{6}} ) + 3 t^2 H_2(z(t)) ( 1 - \frac{\sqrt{2}}{\sqrt{6}} ) + 3 \int_0^t s^2 ds$
so
$\int_0^t z(s)^4 ds = t H_4(z(t)) (1- \frac{2}{\sqrt{6}} ) + 3 t^2 H_2(z(t)) ( 1 - \frac{\sqrt{2}}{\sqrt{6}} ) + t^3$
and we could take $t^3$ common and have a representation in terms of standard normal. And then you can easily simulate the density of $\int_0^t z(s)^4 ds$ or analytically calculate and plot its density.

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

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

Ok friends, there is possibility of a slight bug in the program I posted in post 692. It will work perfectly well with lognormal and calculation of stochastic integrals is perfectly fine. But I am slightly afraid that I might have multiplied stochastic integrals in reverse in the algorithm and that would not matter for lognormal SDE since it takes a constant multiple of initial values. Though the program seems to work fine for some other simple SDEs where I checked but I want to double check more carefully and if there is an error, I will be posting a new version in a few days. But again good thing is that calculation of stochastic integrals is good and possible error is in the algorithm that can be easily fixed and it is a programming issue and not a mathematical or research issue. I will be updating in a day or two.
Last edited by Amin on May 15th, 2018, 9:34 pm, edited 1 time in total.

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

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

Sorry Double.

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

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

It will be a few days before I post a new full fledged program since some of the program architecture has to be completely changed.

But some interesting things about hermite polynomials.

1. If we have two different stochastic integrals or dependent processes of a single SDE and we want to get a density representation of the sum/difference and we know their hermite representations, We can simply add/subtract the respective hermite coefficients as hermite is an orthogonal basis. By coefficients of hermite polynomials, I mean different variable weights that we have assigned to each hermite polynomial when we found the hermite representation. And these weights are not to be confused with the fixed coefficients assigned to different powers in each hermite polynomial.

2. The same applies if we have two different stochastic integrals or dependent processes of two different SDEs driven by a single normal random variable and we want to get a density representation of the sum/difference of stochastic processes of different SDEs driven by the same normal noise and we know their hermite representations, again we can simply add/subtract the respective hermite coefficients/weights as hermite is an orthogonal basis.

3. If we have two different stochastic integrals of totally different SDEs driven by independent normal noises, and we want to add them and get a single density representation for the sum, we could easily do that if we have respective hermite representations of each stochastic integral since there is a map from two different hermite polynomials to a single hermite polynomial of the sum of two processes.

$H_{n}(x+y) = 2^{\frac{-n}{2}} \sum_{k=0}^{n} \binom{n}{k} H_{n-k}(x \sqrt{2}) H_{k}(y \sqrt{2})$

I copied the formula from Wikipedia and I think it is based on a combination of Taylor expansion and Fourier transforms. This is a formula (and its fourier counterparts) that relates hermite representation of a sum of two independent processes in terms of their hermite representation and coefficients of the resulting process can be found from the same formula in terms of coefficients/weights of hermite polynomials of two independent processes. I am sure there can be vast applications of this formula and one very simple application can be finding a single hermite representation of sum of two processes driven by two independent lognormal SDEs or sum of various stochastic integrals of two different lognormal SDEs. And once we have a single hermite representation of a sum of two stochastic processes, we can easily calculate density of such a stochastic process. I am sure there will also be many other interesting applications of this formula.

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

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

I am posting this interim version of the new program. Works perfectly with lognormal and driftless CEV but with  lognormal drift + CEV diffusion, the results are very slightly off when the Ito-Taylor results have converged. As of this time, You cannot do dt-integrals or dz-integrals of diffusions with this program but I will be adding this capability in the next version. I will be posting a new version at the end of next week after carefully checking for the possible error and reason why lognormal drift + CEV SDEs are not perfectly exact as their might be some error in my algorithm. This is just an interim version.
Clear[x, t, x0, Y, P1, P2, l1, l2, Z, mu, sigma, mu0, sigma0, alpha, beta, gamma, a0, k, k1, m, L1, L2, YAns, H, Hc, Hz, X1, fX1, DfX1, Nn, n, HCoeff, ZCoeff, DZCoeff, zz, KIndex, MIndex, MaxIndex, Y1,Y2, CoeffDX1, CoeffDX2, CoeffdtIntegral, CoeffdzIntegral, qIndex1, qIndex21, qIndex22, qIndex23, qIndex0, qqIndex1, qqIndex2, n1, m1, l3, P, k0, m0];
P1 = 12;
P2 = 12;
P3 = 12;
L1 = 0;
L2 = 0;
P = If[P1 > P2, P1, P2];
Array[Y, {P + L1 + 1, P + L1 + 1, P + L2 + 1}, {0, 0, 0}];
For[k = 0, k <= P + L1, k++, (For[m = 0, m <= P + L1, m++, (For[n = 0, n <= P + L2, n++, (Y[k, m, n] = 0);]);]);];
Array[KIndex, {P + 1}, {0}];
For[k = 0, k <= P + 1, k++, If[k > 0, (k0 = k - 1; KIndex[k] = (k0*(k0 + 1)*(2*k0 + 1) + 9*k0*(k0 + 1) + 12*k0)/12 + 1), KIndex[k] = 0]];
Array[MIndex, {P1 + 1}, {0}];
For[m = 0, m <= P + 1, m++, If[m > 0, (m0 = m - 1; MIndex[m] = m0*(m0 + 1)/2 + m0 + 1), MIndex[m] = 0]];
MaxIndex = 1 + 2*MIndex[Round[P/2 + 1]]*KIndex[Round[P/2 + 1]];
Array[Y1, MaxIndex, 0];
Array[Y2, MaxIndex, 0];
For[m = 0, m < MaxIndex, m++, (Y1[m] = 1; Y2[m] = 0;)];
Array[CoeffdtIntegral, {P + 1, P + 1}, {0, 0}];
Array[CoeffdzIntegral, {P + 1, P + 1}, {0, 0}];
For[k = 0, k < P, k++, (For[m = 0, m <= k, m++, (CoeffdtIntegral[k - m, m] = 1/(m + 1)*(1 - Sqrt[k - m]/Sqrt[2*m + (k - m) + 2]);
CoeffdzIntegral[k - m, m] = 1/Sqrt[(2*m + (k - m) + 1)*(k - m + 1)];)])];
a0 = 1;
mu0 = .25;
sigma0 = .75;
alpha = 1;
gamma = 1.0;
beta = 1;
Y[L1, L1, L2] = a0;
For[k = 0, k < P, k++, (For[m = 0, m <= k, m++, (l3 = L2 + k - m;
For[n = 0, n <= m, n++, (l1 = L1 + m - n; l2 = L1 + n; qIndex1 =
MIndex[m] + n + 1; QIndex1 = (qIndex1 - 1)*KIndex[P - k + 1];
qIndex21 = MIndex[m + 1] + n + 1;
QIndex21 = (qIndex21 - 1)*KIndex[P - k];
qIndex22 = MIndex[m + 1] + (n + 1) + 1;
QIndex22 = (qIndex22 - 1)*KIndex[P - k];
qIndex23 = MIndex[m] + n + 1;
QIndex23 = (qIndex23 - 1)*KIndex[P - k];
CoeffDX1 = (alpha + (l1 - L1)*beta + (l2 - L1)*2*gamma + (l3 - L2)*gamma - (l1 - L1) - 2*(l2 - L1) - (l3 - L2));
CoeffDX2 = CoeffDX1 - 1; MuDX1 = CoeffDX1*mu0;
Sigma2DX2 = .5*CoeffDX1*CoeffDX2*sigma0^2;
For[k1 = 0, k1 < P - k, k1++, (For[m1 = 0, m1 <= k1, m1++, (For[n1 = 0, n1 <= m1, n1++, (qqIndex2 = KIndex[k1] + MIndex[m1] + n1 + 1;
qqIndex1 = KIndex[k1 + 1] + MIndex[m1 + 1] + n1 + 1;
Y2[qqIndex2 + QIndex21] = Y2[qqIndex2 + QIndex21] + Y1[qqIndex1 + QIndex1]*CoeffdtIntegral[k1 - m1, m1]*MuDX1;
qqIndex1 = KIndex[k1 + 1] + MIndex[m1 + 1] + (n1 + 1) + 1;
Y2[qqIndex2 + QIndex22] = Y2[qqIndex2 + QIndex22] + Y1[qqIndex1 + QIndex1]*CoeffdtIntegral[k1 - m1, m1]*Sigma2DX2;
qqIndex1 = KIndex[k1 + 1] + MIndex[m1] + n1 + 1;
Y2[qqIndex2 + QIndex23] = Y2[qqIndex2 + QIndex23] + Y1[qqIndex1 + QIndex1]*CoeffdzIntegral[k1 - m1, m1]*SigmaDX1;)];)])];
Y[l1 + 1, l2, l3] = Y2[1 + QIndex21];
Y[l1, l2 + 1, l3] = Y2[1 + QIndex22];
Y[l1, l2, l3 + 1] = Y2[1 + QIndex23];)];)];
For[qIndex0 = 0, qIndex0 < MaxIndex, qIndex0++, (Y1[qIndex0] = Y2[qIndex0]; Y2[qIndex0] = 0;)];)];
Array[HCoeff, P2 + L2 + 1, 0];
For[k = 0, k < P2 + L2 + 1, k++, HCoeff[k] = 0;]
x = 1.0;
t = 1.0;
For[k = 0, k <= P3, k++, (For[m = 0, m <= If[k < P1, k, P1], m++, (l0 = L1 + m; l3 = L2 + k - m;For[n = 0, n <= If[m < P2, m, P2], n++, (l1 = L1 + m - n; l2 = L1 + n;
HCoeff[l3] = HCoeff[l3] + SetPrecision[Y[l1, l2, l3]*x^(alpha + (l1 - L1)*beta + (l2 - L1)*2*gamma + (l3 - L2)*gamma - (l1 - L1) - 2*(l2 - L1) - (l3 - L2))*(Sqrt[t^(2*l1 + 2*l2 - 2*L1 + l3)]), 50];)];)];);];
Array[Hz, {P2 + L2 + 1, P2 + L2 + 1}, {0, 0}];
For[k = 0, k < P2 + L2 + 1, k++, For[m = 0, m < P2 + L2 + 1, m++, Hz[k, m] = 0;]];
Hz[0, 0] = 1.0;
For[k = 0, k < P2 + L2, k++, (For[m = 0, m <= k, m++, (Hz[k + 1, m + 1] = SetPrecision[Hz[k, m], 50];
If[m >= 1, Hz[k + 1, m - 1] = SetPrecision[Hz[k + 1, m - 1] - m*Hz[k, m], 50], 0];)])];
Array[ZCoeff, P2 + L2 + 1, 0];
For[k = 0, k < P2 + L2 + 1, k++, ZCoeff[k] = 0];
For[k = 0, k < P2 + L2 + 1, k++, (For[m = 0, m <= k, m++, ZCoeff[m] = ZCoeff[m] + SetPrecision[Hz[k, m]*HCoeff[k], 50]])];
Array[DZCoeff, P2 + L2, 0];
For[k = 0, k < P2 + L2, k++, DZCoeff[k] = 0];
For[k = 0, k < P2 + L2, k++, DZCoeff[k] = SetPrecision[ZCoeff[k + 1]*(k + 1), 50]];
Nn = 81;
Array[X1, Nn, 0];
Array[DfX1, Nn, 0];
Array[fX1, Nn, 0];
Array[zz, Nn, 0];
For[n = 0, n < Nn, n++, (X1[n] = 0; DfX1[n] = 0; fX1[n] = 0; zz[n] = 0;)];
For[n = 0, n < Nn, n++, (zz[n] = n*.1 - 4; X1[n] = 0; DfX1[n] = 0; fX1[n] = 0;
X1[n] = ZCoeff[0]; DfX1[n] = DZCoeff[0];
For[k = 1, k < P2 + L2 + 1, k++, (X1[n] = X1[n] + SetPrecision[ZCoeff[k]*zz[n]^k, 50];
If[k < P2 + L2, DfX1[n] = DfX1[n] + SetPrecision[DZCoeff[k]*zz[n]^k, 50], 0];)];
fX1[n] = SetPrecision[PDF[NormalDistribution[0, 1], zz[n]]/Abs[DfX1[n]], 50];)];
data1 = Table[{X1[n], fX1[n]}, {n, 0, 75}];
p1 = ListLinePlot[data1]
Here is the code snippet for plotting the density of lognormal drift + CEV diffusion SDE for comparison. For mu0=0, it will create and plot the density of driftless CEV.
Nn1 = 1000;
Array[yy, Nn1, 0];
Array[fyy, Nn1, 0];
gamma2 = 2*gamma;
X0 = x;
tau = If[mu0 != 0, (Exp[2*mu0*gamma0*t] - 1)/(2*mu0*gamma0), t];
kStar = (2*mu0)/(sigma0^2*(2 - gamma2)*(Exp[mu0*(2 - gamma2)*t] - 1));
x0 = kStar*X0^(2 - gamma2)*Exp[mu0*(2 - gamma2)*t];
For[n = 0, n < Nn1, n++, (yy[n] = (.0002 + n*.0050); YY[n] = kStar*yy[n]^(2 - gamma2);
fyy[n] = (2 - gamma2)*kStar^(1/(2 - gamma2))*(x0*YY[n]^(1 - 2*gamma2))^(1/(2*(2 - gamma2)))*BesselI[Abs[1/(2 - gamma2)], (2*(YY[n]*x0)^.5)]*Exp[-1*(x0 + YY[n])];)];
data10 = Table[{yy[n], fyy[n]}, {n, 0, 1000}];
p10 = ListLinePlot[data10, PlotRange -> All]

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

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

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

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

Updating about the recent progress regarding the work. I was able to work on the method in which I could analytically calculate weights of hermite polynomial expansion analytically using only first four hermites. And then starting from origin, we could branch a tree and then recalculate the coefficients on each branch of the tree using the analytic hermite expansion every .25 years. And we could use the locally calculated hermite weights for the particular single branch (related to a specific value of normal) of the tree. So basically after the first time step from origin, we calculate the different hermite weights according to local values of SDE variable along the tree using the same analytic hermite expansion formula and we continue to recalculate the local hermite weights along every branch after every .25 years. And the method works for several SDEs. But I have not extensively tried everything and I was sort of excited to report this new improvement in the method. I will be doing a few more tests for next day or two with different SDEs and then post the new mathematica program here for review by friends. Good thing is that it is based on only first four hermite polynomials and then takes local coefficietns/weights along the tree and we also have to rescale the variances as we add normal. I hope to post the full-fledged program in a few days.

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

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

Here is the quick raw program. Sometimes, you could see random lines in the graph coming out from origin and they are due to diffusion paths getting into negative territory and continuing and when I take derivatives numerically, these places can give unstable results. Can easily be corrected by making negative paths zero and then not continuing them. This will be perfectly resolved in the next version so please pardon these little issues. But I will work on all these aspects and release an industrial strength program in a few days. This is just an initial raw version.
You can overlay this program's graph with the mathematica code snippet for CEV diffusions with lognormal type drift (that I previously posted)  if you will like to see the results with that SDE type. I believe we are done with a lot of difficult part and hard work and I will be explaining the program and its working carefully in next few days.
Here is the code. Please note that Tt is number of time intervals in the program. Here t=.25 nand Tt=12 so Terminal time is 3.0
Clear[x, t, x0, Y, P1, P2, l1, l2, Z, mu, sigma, mu0, sigma0, alpha, beta, gamma, a0, k, k1, m, L1, L2, YAns, H, Hc, Hz, X1, fX1, DfX1, Nn, n, HCoeff, ZCoeff, DZCoeff, zz, KIndex, MIndex, MaxIndex, Y1, Y2, CoeffDX1, CoeffDX2, CoeffdtIntegral, CoeffdzIntegral, qIndex1, qIndex21, qIndex22, qIndex23, qIndex0, qqIndex1, qqIndex2, n1, m1, l3, P, k0, m0];
P1 = 6;
P2 = 6;
P3 = 6;
L1 = 0;
L2 = 0;
P = If[P1 > P2, P1, P2];
Array[Y, {P + L1 + 1, P + L1 + 1, P + L2 + 1}, {0, 0, 0}];
For[k = 0, k <= P + L1, k++, (For[m = 0, m <= P + L1, m++, (For[n = 0, n <= P + L2, n++, (Y[k, m, n] = 0);]);]);];
Array[KIndex, {P + 1}, {0}];
For[k = 0, k <= P + 1, k++, If[k > 0, (k0 = k - 1; KIndex[k] = (k0*(k0 + 1)*(2*k0 + 1) + 9*k0*(k0 + 1) + 12*k0)/12 + 1), KIndex[k] = 0]];
Array[MIndex, {P1 + 1}, {0}];
For[m = 0, m <= P + 1, m++, If[m > 0, (m0 = m - 1; MIndex[m] = m0*(m0 + 1)/2 + m0 + 1), MIndex[m] = 0]];
MaxIndex = 1 + 2*MIndex[Round[P/2 + 1]]*KIndex[Round[P/2 + 1]];
Array[Y1, MaxIndex, 0];
Array[Y2, MaxIndex, 0];
For[m = 0, m < MaxIndex, m++, (Y1[m] = 1; Y2[m] = 0;)];
Array[CoeffdtIntegral, {P + 1, P + 1}, {0, 0}];
Array[CoeffdzIntegral, {P + 1, P + 1}, {0, 0}];
For[k = 0, k < P, k++, (For[m = 0, m <= k, m++, (CoeffdtIntegral[k - m, m] = 1/(m + 1)*(1 - Sqrt[k - m]/Sqrt[2*m + (k - m) + 2]);
CoeffdzIntegral[k - m, m] = 1/Sqrt[(2*m + (k - m) + 1)*(k - m + 1)];)])];
a0 = 1;
mu0 = .250;
sigma0 = .55;
alpha = 1;
gamma = .65;
beta = 1;
Y[L1, L1, L2] = a0;
For[k = 0, k < P, k++, (For[m = 0, m <= k, m++, (l3 = L2 + k - m;
For[n = 0, n <= m, n++, (l1 = L1 + m - n; l2 = L1 + n; qIndex1 =
MIndex[m] + n + 1; QIndex1 = (qIndex1 - 1)*KIndex[P - k + 1];
qIndex21 = MIndex[m + 1] + n + 1;
QIndex21 = (qIndex21 - 1)*KIndex[P - k];
qIndex22 = MIndex[m + 1] + (n + 1) + 1;
QIndex22 = (qIndex22 - 1)*KIndex[P - k];
qIndex23 = MIndex[m] + n + 1;
QIndex23 = (qIndex23 - 1)*KIndex[P - k];
CoeffDX1 = (alpha + (l1 - L1)*beta + (l2 - L1)*2*
gamma + (l3 - L2)*gamma - (l1 - L1) -
2*(l2 - L1) - (l3 - L2));
CoeffDX2 = CoeffDX1 - 1; MuDX1 = CoeffDX1*mu0;
Sigma2DX2 = .5*CoeffDX1*CoeffDX2*sigma0^2;
For[k1 = 0, k1 < P - k, k1++, (For[m1 = 0, m1 <= k1, m1++, (For[n1 = 0, n1 <= m1, n1++, (qqIndex2 = KIndex[k1] + MIndex[m1] + n1 + 1;
qqIndex1 = KIndex[k1 + 1] + MIndex[m1 + 1] + n1 + 1;
Y2[qqIndex2 + QIndex21] = Y2[qqIndex2 + QIndex21] + Y1[qqIndex1 + QIndex1]*CoeffdtIntegral[k1 - m1, m1]*MuDX1; qqIndex1 = KIndex[k1 + 1] + MIndex[m1 + 1] + (n1 + 1) + 1;
Y2[qqIndex2 + QIndex22] = Y2[qqIndex2 + QIndex22] + Y1[qqIndex1 + QIndex1]*CoeffdtIntegral[k1 - m1, m1]*Sigma2DX2;
qqIndex1 = KIndex[k1 + 1] + MIndex[m1] + n1 + 1;
Y2[qqIndex2 + QIndex23] = Y2[qqIndex2 + QIndex23] + Y1[qqIndex1 + QIndex1]*CoeffdzIntegral[k1 - m1, m1]*SigmaDX1;)];)])];
Y[l1 + 1, l2, l3] = Y2[1 + QIndex21];
Y[l1, l2 + 1, l3] = Y2[1 + QIndex22];
Y[l1, l2, l3 + 1] = Y2[1 + QIndex23];)];)];
For[qIndex0 = 0, qIndex0 < MaxIndex,
qIndex0++, (Y1[qIndex0] = Y2[qIndex0]; Y2[qIndex0] = 0;)];)];
Nn = 81;
P2 = 6;
Array[HCoeff, {P2 + L2 + 1, Nn}, {0, 0}];
Array[Hz, {P2 + L2 + 1, P2 + L2 + 1}, {0, 0}];
For[k = 0, k < P2 + L2 + 1, k++, For[m = 0, m < P2 + L2 + 1, m++, Hz[k, m] = 0;]];
Hz[0, 0] = 1.0;
For[k = 0, k < P2 + L2, k++, (For[m = 0, m <= k, m++, (Hz[k + 1, m + 1] = Hz[k, m];
If[m >= 1, Hz[k + 1, m - 1] = Hz[k + 1, m - 1] - m*Hz[k, m], 0];)])];
x0 = 1;
t = .250;
Tt = 12;
For[k = 0, k < P2 + L2 + 1, k++, (For [nn = 0, nn < Nn, nn++, HCoeff[k, nn] = 0;])];
Array[x, Nn, 0];
For[nn = 0, nn < Nn, nn++, x[nn] = x0];
Array[ZCoeff, {P2 + L2 + 1, Nn}, 0];
Array[X1, Nn, 0];
Array[DfX1, Nn, 0];
Array[fX1, Nn, 0];
Array[zz, Nn, 0];
For[nn = 0, nn < Nn, nn++, (X12[nn] = 0; DfX12[nn] = 0;)];
Array[DZCoeff, {P2 + L2, Nn}, {0, 0}];
For[tt = 0, tt < Tt, tt++, (For[nn = 0, nn < Nn, nn++, (For[k = 0, k <= P3, k++, HCoeff[k, nn] = 0];
For[k = 0, k <= P3, k++, (For[m = 0, m <= If[k < P1, k, P1], m++, (l0 = L1 + m; l3 = L2 + k - m;
For[n = 0, n <= If[m < P2, m, P2], n++, (l1 = L1 + m - n; l2 = L1 + n;
HCoeff[l3, nn] = HCoeff[l3, nn] + Y[l1, l2, l3]*x[nn]^(alpha + (l1 - L1)*beta + (l2 - L1)*2*gamma + (l3 - L2)*gamma - (l1 - L1) - 2*(l2 - L1) - (l3 - L2))*t^(l1 + l2 - L1 + .5*l3);)];)])])];
For[k = 0, k < P2 + L2 + 1, k++, For[nn = 0, nn < Nn, nn++, ZCoeff[k, nn] = 0]];
For[nn = 0, nn < Nn, nn++, (For[k = 0, k < P2 + L2 + 1, k++,
For[m = 0, m <= k, m++, ZCoeff[m, nn] = ZCoeff[m, nn] + Hz[k, m]*HCoeff[k, nn]]];)];
For[nn = 0, nn < Nn, nn++, (X12[nn] = ZCoeff[0, nn]; zz[nn] = nn*.1 - 4;
For[k = 1, k < P2 + L2 + 1, k++,X12[nn] = X12[nn] + ZCoeff[k, nn]*(1/Tt^(.5)*zz[nn])^k;]; x[nn] = X12[nn];)];)];
For[nn = 1, nn < Nn - 1, nn++, zz[nn] = Tt^(.5)*(nn*.1 - 4);];
For[nn = 1, nn < Nn - 1, nn++, (DfX12[nn] = (X12[nn + 1] - X12[nn - 1])/(zz[n + 1] - zz[n - 1]);)];
For[nn = 0, nn < Nn, nn++, (zz[nn] = nn*.1 - 4;
fX12[nn] = PDF[NormalDistribution[0, Sqrt[Tt]], ((Tt)^(.5)*zz[nn])]/Abs[DfX12[nn]];)];
data12 = Table[{X12[nn], fX12[nn]}, {nn, 0, 80}];
p12 = ListLinePlot[data12]
Here is the output of the SDE overlaid with graph from analytic formula for the density. The SDE is dX=mu0* X dt+ sigma0 * X^gamma *dz(t)
with mu0=.25, sigma0=.55, gamma=.65. and T=3. There are Tt=12 time levels of t=.25 years duration each where the analytic coefficients are recomputed for each branch of the SDE tree based on a particular value of normal ( which ranges from -4 to +4 with .1 increments for each branch of the tree). Ok here is the overlaid graph comparison for the output with program above posted.

I am sure more work will need to be done but these results are very encouraging.

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

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

Friends, this is slightly updated version of the program that corrects for instabilities in the graph that arise from Ito-Taylor tree branches going into negative. Most of the graphs will be perfectly good now. I have also made a few minor changes. Here is the code.
Clear[x, t, x0, Y, P1, P2, l1, l2, Z, mu, sigma, mu0, sigma0, alpha, beta, gamma, a0, k, k1, m, L1, L2, YAns, H, Hc, Hz, X1, fX1, DfX1, Nn, n, HCoeff, ZCoeff, DZCoeff, zz, KIndex, MIndex, MaxIndex, Y1, Y2, CoeffDX1, CoeffDX2, CoeffdtIntegral, CoeffdzIntegral, qIndex1, qIndex21, qIndex22, qIndex23, qIndex0, qqIndex1, qqIndex2, n1, m1, l3, P, k0, m0];
P1 = 4;
P2 = 4;
P3 = 4;
L1 = 0;
L2 = 0;
P = If[P1 > P2, P1, P2];
Array[Y, {P + L1 + 1, P + L1 + 1, P + L2 + 1}, {0, 0, 0}];
For[k = 0, k <= P + L1, k++, (For[m = 0, m <= P + L1, m++, (For[n = 0, n <= P + L2, n++, (Y[k, m, n] = 0);]);]);];
Array[KIndex, {P + 1}, {0}];
For[k = 0, k <= P + 1, k++, If[k > 0, (k0 = k - 1;
KIndex[k] = (k0*(k0 + 1)*(2*k0 + 1) + 9*k0*(k0 + 1) + 12*k0)/12 + 1), KIndex[k] = 0]];
Array[MIndex, {P1 + 1}, {0}];
For[m = 0, m <= P + 1, m++, If[m > 0, (m0 = m - 1; MIndex[m] = m0*(m0 + 1)/2 + m0 + 1),
MIndex[m] = 0]];
MaxIndex = 1 + 2*MIndex[Round[P/2 + 1]]*KIndex[Round[P/2 + 1]];
Array[Y1, MaxIndex, 0];
Array[Y2, MaxIndex, 0];
For[m = 0, m < MaxIndex, m++, (Y1[m] = 1; Y2[m] = 0;)];
Array[CoeffdtIntegral, {P + 1, P + 1}, {0, 0}];
Array[CoeffdzIntegral, {P + 1, P + 1}, {0, 0}];
For[k = 0, k < P, k++, (For[m = 0, m <= k, m++, (CoeffdtIntegral[k - m, m] = 1/(m + 1)*(1 - Sqrt[k - m]/Sqrt[2*m + (k - m) + 2]);
CoeffdzIntegral[k - m, m] = 1/Sqrt[(2*m + (k - m) + 1)*(k - m + 1)];)])];
a0 = 1;
mu0 = .2;
sigma0 = 1.0;
alpha = 1;
gamma = .65;
beta = 1;
Y[L1, L1, L2] = a0;
For[k = 0, k < P, k++, (For[m = 0, m <= k, m++, (l3 = L2 + k - m;
For[n = 0, n <= m, n++, (l1 = L1 + m - n; l2 = L1 + n; qIndex1 =
MIndex[m] + n + 1; QIndex1 = (qIndex1 - 1)*KIndex[P - k + 1];
qIndex21 = MIndex[m + 1] + n + 1;
QIndex21 = (qIndex21 - 1)*KIndex[P - k];
qIndex22 = MIndex[m + 1] + (n + 1) + 1;
QIndex22 = (qIndex22 - 1)*KIndex[P - k];
qIndex23 = MIndex[m] + n + 1;
QIndex23 = (qIndex23 - 1)*KIndex[P - k];
CoeffDX1 = (alpha + (l1 - L1)*beta + (l2 - L1)*2*gamma + (l3 - L2)*gamma - (l1 - L1) - 2*(l2 - L1) - (l3 - L2));
CoeffDX2 = CoeffDX1 - 1; MuDX1 = CoeffDX1*mu0;
Sigma2DX2 = .5*CoeffDX1*CoeffDX2*sigma0^2;
For[k1 = 0, k1 < P - k, k1++, (For[m1 = 0, m1 <= k1, m1++, (For[n1 = 0, n1 <= m1, n1++, (qqIndex2 = KIndex[k1] + MIndex[m1] + n1 + 1;
qqIndex1 = KIndex[k1 + 1] + MIndex[m1 + 1] + n1 + 1;
Y2[qqIndex2 + QIndex21] = Y2[qqIndex2 + QIndex21] + Y1[qqIndex1 + QIndex1]*CoeffdtIntegral[k1 - m1, m1]*MuDX1; qqIndex1 =
KIndex[k1 + 1] + MIndex[m1 + 1] + (n1 + 1) + 1;
Y2[qqIndex2 + QIndex22] = Y2[qqIndex2 + QIndex22] + Y1[qqIndex1 + QIndex1]*CoeffdtIntegral[k1 - m1, m1]*Sigma2DX2;
qqIndex1 = KIndex[k1 + 1] + MIndex[m1] + n1 + 1;

Y2[qqIndex2 + QIndex23] = Y2[qqIndex2 + QIndex23] + Y1[qqIndex1 + QIndex1]*CoeffdzIntegral[k1 - m1, m1]*SigmaDX1;)];)])];
Y[l1 + 1, l2, l3] = Y2[1 + QIndex21];
Y[l1, l2 + 1, l3] = Y2[1 + QIndex22];
Y[l1, l2, l3 + 1] = Y2[1 + QIndex23];)];)];
For[qIndex0 = 0, qIndex0 < MaxIndex, qIndex0++, (Y1[qIndex0] = Y2[qIndex0]; Y2[qIndex0] = 0;)];)];
Nn = 81;
P2 = 4;
Array[HCoeff, {P2 + L2 + 1, Nn}, {0, 0}];
Array[Hz, {P2 + L2 + 1, P2 + L2 + 1}, {0, 0}];
For[k = 0, k < P2 + L2 + 1, k++,
For[m = 0, m < P2 + L2 + 1, m++, Hz[k, m] = 0;]];
Hz[0, 0] = 1.0;
For[k = 0, k < P2 + L2, k++, (For[m = 0, m <= k, m++, (Hz[k + 1, m + 1] = Hz[k, m];
If[m >= 1, Hz[k + 1, m - 1] = Hz[k + 1, m - 1] - m*Hz[k, m], 0];)])];
x0 = 1;
t = .1250;
Tt = 8;
For[k = 0, k < P2 + L2 + 1, k++, (For [nn = 0, nn < Nn, nn++, HCoeff[k, nn] = 0;])];
Array[x, Nn, 0];
For[nn = 0, nn < Nn, nn++, x[nn] = x0];
Array[ZCoeff, {P2 + L2 + 1, Nn}, 0];
Array[X1, Nn, 0];
Array[DfX1, Nn, 0];
Array[fX1, Nn, 0];
Array[zz, Nn, 0];
For[nn = 0, nn < Nn, nn++, (X12[nn] = 0; DfX12[nn] = 0;)];
Array[DZCoeff, {P2 + L2, Nn}, {0, 0}];
For[tt = 0, tt < Tt, tt++, (For[nn = 0, nn < Nn, nn++, (For[k = 0, k <= P3, k++, HCoeff[k, nn] = 0];
For[k = 0, k <= P3, k++, (For[m = 0, m <= If[k < P1, k, P1], m++, (l0 = L1 + m; l3 = L2 + k - m;
For[n = 0, n <= If[m < P2, m, P2], n++, (l1 = L1 + m - n; l2 = L1 + n;
HCoeff[l3, nn] = HCoeff[l3, nn] + Y[l1, l2, l3]*x[nn]^(alpha + (l1 - L1)*beta + (l2 - L1)*2*gamma + (l3 - L2)*gamma - (l1 - L1) - 2*(l2 - L1) - (l3 - L2))*t^(l1 + l2 - L1 + .5*l3);)];)])])];
For[k = 0, k < P2 + L2 + 1, k++, For[nn = 0, nn < Nn, nn++, ZCoeff[k, nn] = 0]];
For[nn = 0, nn < Nn, nn++, (For[k = 0, k < P2 + L2 + 1, k++,
For[m = 0, m <= k, m++, ZCoeff[m, nn] = ZCoeff[m, nn] + Hz[k, m]*HCoeff[k, nn]]];)];
For[nn = 0, nn < Nn, nn++, (X12[nn] = ZCoeff[0, nn]; zz[nn] = nn*.1 - 4;
For[k = 1, k < P2 + L2 + 1, k++, X12[nn] = X12[nn] + ZCoeff[k, nn]*(1/(Tt)^(.5)*zz[nn])^k;]; x[nn] = X12[nn];)];)];
CFlag = 1;
For[nn = Nn - 1, nn > 1, nn--, (If[((X12[nn] > X12[nn - 1]) && (CFlag == 1)), (X13[nn] = X12[nn]; nnStart = nn;), (X13[nn] = 0; CFlag = 0;)])];
nnStart
For[nn = 1, nn < Nn - 1, nn++, zz[nn] = (nn*.1 - 4);];
For[nn = nnStart, nn < Nn - 1, nn++, (
DfX12[nn] = (X12[nn + 1] - X12[nn - 1])/(zz[n + 1] - zz[n - 1]);)];
For[nn = nnStart, nn < Nn, nn++, (zz[nn] = nn*.1 - 4;
fX12[nn] = PDF[NormalDistribution[0, 1], zz[nn]]/Abs[DfX12[nn]];)];
data12 = Table[{X12[nn], fX12[nn]}, {nn, nnStart, 60}];
p12 = ListLinePlot[data12]
Here is the output graph of the Ito-Taylor tree program overlaid with graph from true solution using analytic formula.