Serving the Quantitative Finance Community

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

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

I have started working on backing out volatility from FPE equation that could be input to transition probabilities grid. Meanwhile, I worked a little bit more with optimization of coefficients for various terms of analytic solution of FPE equation. These new coefficients work excellent with mean-correction. I have not perfectly optimized for them but they capture the shape of the density in an excellent way.
Here is the value of coefficients that I have copied from my program

C22(wnStart:Nn)=-.2*(1-gamma)+.0875*sigma0^2./((1-gamma)*w(wnStart:Nn));
C33(wnStart:Nn)=0.0;
C44(wnStart:Nn)=0.1050;

I will do a definitive optimization with a new algorithm after completing the part to back out volatility from FPE. But with mean-correction, you can still use the program for everything.
Please note that C22 takes bessel drift term with a slight modification.

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

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

This is a post I made in off-topic and I am posting it here since I consider it an important message for good Americans and good people of other countries and therefore I copied it here..

https://forum.wilmott.com/viewtopic.php?f=15&t=94796&p=863284#p863284

I also make these claim again that most of these corrupt people have absolutely no sincere intention of doing many things that are in true American interest. With so many corrupt Pakistani generals on their payroll, they still do not force them to stop negative activities in Afghanistan as long as Pakistani generals continue to follow their (based on their hardliner hatred) demands in damaging the interests of our own country, Pakistan. So it seems that it is more important for hardliner generals in American military to promote their agendas of personal hatred than promoting genuine US interests or establishing cherished American/human values and practices in our backward countries. Instead of stopping negative activities by our army that they easily can if they were earnest, these corrupt generals in US army continue to cite these negative activities by our army to increase their funding and destroy reputation of our good-natured public so that they can continue on the other front of pursuing their hateful agendas against the civilian public in our countries.
Can Americans believe that Shakeel Afridi who helped US military in reaching the greatest terrorist on earth still rots in jail despite Trump made tall claims of getting him back in US in a matter of hours. But later convinced by other hardliner conservatives, he totally gave up on the idea since a lot of people in US do not want Shakeel Afridi to be free. In a sophisticated US society, people like Shakeel Afridi would otherwise fail to earn much respect but once good people in US know that this was the guy who helped get the most dangerous terrorist, they would be full of praise and respect for him and this could by association help improve the image of many other good and peaceful Muslim Americans like him and many hardliner conservatives do not want anything like that to ever happen and therefore want Shakeel Afridi to instead rot in jail. I really want to ask good Americans to suggest to the new president Joe Biden to ask the corrupt crooks in US military to finally do one good thing by asking the Pakistani generals on their payroll to free Shakeel Afridi which is a genuine American demand that these corrupt crooks could have easily achieved several long years ago if they were following genuine American interests instead of the agendas based on their personal hatred.
And Diplomatic staff in many countries like mine are directly chosen by people in American army who have no other interest other than finding diplomats who would best help promote their agendas of hardliner conservative hatred by damaging the rightful interest of Pakistani people. I really want to request good Americans when diplomats are shuffled after new president takes oath, please find a nice gentleman/woman for our country who would not only promote genuine friendship between two countries, but also try to advance the good human/American values of democracy, freedom of speech, secularism and all other good things that modern humanity cherishes. I know good Americans are probably the nicest people on earth and I hope they will always be able to control bad elements in their society(that by the way exist in every society). Bad people are only able to prevail in many good and educated societies because good people do not know the complete truth about actions of bad people.
Tomorrow, my mother would be operated for angiography in AFIC by sister's family friend, I hope mind control agents would take it easy on their persecution and drugging of water while I would have to be with my mother.

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

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

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

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

$\frac{\partial p(X)}{\partial t}=\frac{\partial [p(Z) \frac{\partial Z}{\partial X}]}{\partial t}= \frac{\partial p(Z)}{\partial Z} {(\frac{\partial Z}{\partial X})}^2 \frac{\partial X}{\partial t}+ p(Z) \frac{\partial^2 Z}{\partial X^2} \frac{\partial X}{\partial t}$
$=Zp(Z) {(\frac{\partial Z}{\partial X})}^2 \frac{\partial X}{\partial t}+ p(Z) \frac{\partial^2 Z}{\partial X^2} \frac{\partial X}{\partial t}$   Equation(9)

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

Now we write the density of X(t) in terms of density of standard normal Z by substituting from above equations (6-9) in Fokker-Planck Equation (10)

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

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

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

We multiply both sides of equations with dt and rearrange the above equation(12) to write

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

We multiply both sides of the above equation 13 with ${(\frac{\partial X}{\partial Z})}^3$ to get

$dX(Z (\frac{\partial X}{\partial Z}) + {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^2 Z}{\partial X^2} )$
$= -\mu(X) (Z (\frac{\partial X}{\partial Z})+ {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^2 Z}{\partial X^2}) dt$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} (Z (\frac{\partial X}{\partial Z})+ {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^2 Z}{\partial X^2}) dt$
$- .5 {\sigma(X)}^2 dt$
$+ .5 {\sigma(X)}^2 (Z^2 -3 Z {(\frac{\partial X}{\partial Z})}^2 \frac{\partial^2 Z}{\partial X^2} +{(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3}) dt$
$- \frac{\partial \mu(X)}{\partial X} {(\frac{\partial X}{\partial Z})}^2 dt$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) {(\frac{\partial X}{\partial Z})}^2 dt$  Equation(14)

We use the relationship  ${(\frac{\partial X}{\partial Z})}^3 \frac{\partial^2 Z}{\partial X^2} =- \frac{\partial^2 X}{\partial Z^2}$ in above equation (14) to get
the following equation

$dX(Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2} )$
$=- \mu(X) (Z (\frac{\partial X}{\partial Z})- \frac{\partial^2 X}{\partial Z^2}) dt$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} (Z (\frac{\partial X}{\partial Z})- \frac{\partial^2 X}{\partial Z^2}) dt$
$- .5 {\sigma(X)}^2 dt$
$+ .5 {\sigma(X)}^2 (Z^2 +3 Z {(\frac{\partial Z}{\partial X})} \frac{\partial^2 X}{\partial Z^2} +{(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3}) dt$
$- \frac{\partial \mu(X)}{\partial X} {(\frac{\partial X}{\partial Z})}^2 dt$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) {(\frac{\partial X}{\partial Z})}^2 dt$  Equation(15)

The first four terms of the equation can be solved to give a differential of following squared form

$dX(Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2} )$
$= -\mu(X) (Z (\frac{\partial X}{\partial Z})- \frac{\partial^2 X}{\partial Z^2}) dt$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} (Z (\frac{\partial X}{\partial Z})- \frac{\partial^2 X}{\partial Z^2}) dt$
$- .5 {\sigma(X)}^2 dt$
$=.5 d\Big[\int_{t_0}^{t_1} dw - \int_{t_0}^{t_1} \mu(X) ds - \int_{t_0}^{t_1} 2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} ds$
$+ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]^2$  Equation(16)

I want to explain that the above term (equation 16)
$.5 d\Big[\int_{t_0}^{t_1} dX - \int_{t_0}^{t_1} \mu(X) ds - \int_{t_0}^{t_1} 2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} ds$
$+ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]^2$
equals
$+.5 \Big[( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]^2$ at time $t_0$ since the first three integral terms vanish at start time.
and equation 16 equals at time $t_1$
$.5 \Big[X(t_1)-X(t_0) - \int_{t_0}^{t_1} \mu(X) ds - \int_{t_0}^{t_1} 2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} ds$
$+ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]^2$

So we can write the above equation 15 after substituting equation 16 in it as
$.5 d\Big[\int_{t_0}^{t_1} dX - \int_{t_0}^{t_1} \mu(X) ds - \int_{t_0}^{t_1} 2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} ds$
$+ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]^2$  Equation(16)
$=+ .5 {\sigma(X)}^2 (Z^2 +3 Z {(\frac{\partial Z}{\partial X})} \frac{\partial^2 X}{\partial Z^2} +{(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3}) dt$
$- \frac{\partial \mu(X)}{\partial X} {(\frac{\partial X}{\partial Z})}^2 dt$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) {(\frac{\partial X}{\partial Z})}^2 dt$  Equation(17)

We can write the above equation after noting that some terms are getting constant multiplier coefficients not immediately known. So we re-write equation (17) as

$.5 \Big[X(t_1)-X(t_0) - \int_{t_0}^{t_1} \mu(X) ds - \int_{t_0}^{t_1} 2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} ds$
$+ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]^2$
$= +{\Big[ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]}^2$
$+ .5 {\sigma(X)}^2 (Z^2 +3C_2 Z {(\frac{\partial Z}{\partial X})} \frac{\partial^2 X}{\partial Z^2} +C_4 {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3}) dt$
$- C_3\frac{\partial \mu(X)}{\partial X} {(\frac{\partial X}{\partial Z})}^2 dt$
$+C_5({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) {(\frac{\partial X}{\partial Z})}^2 dt$  Equation(18)

After taking square root on both sides and re-arranging, we get the evolution equation for X(t) as

$X(t_1)=X(t_0) + \int_{t_0}^{t_1} \mu(X) ds + \int_{t_0}^{t_1} 2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} ds$
$- Z (\frac{\partial X}{\partial Z}) + \frac{\partial^2 X}{\partial Z^2}$
$+\sqrt{\Big[ +{\big[ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\big]}^2 + {\sigma(X)}^2 \big[ Z^2 +3C_2 Z {(\frac{\partial Z}{\partial X})} \frac{\partial^2 X}{\partial Z^2} +C_4 {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3} \big] dt - C_3\frac{\partial \mu(X)}{\partial X} {(\frac{\partial X}{\partial Z})}^2 dt+C_5({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) {(\frac{\partial X}{\partial Z})}^2 dt \Big]}$
Equation(19)

Please note that all the terms on RHS appear under square-root. The last matlab program I posted has two terms outside the square-root and are being linearly added and that is an error. I will post a new program with this solution in a day and it works far better closer to zero.
Here I mention that Lamperti/Bessel form SDE is given as
$dX(t)=\mu(X) dt + \sigma dz(t)$  Equation(20)
and it has a solution given as

$X(t_1)=X(t_0) + \int_{t_0}^{t_1} \mu(X) ds - Z (\frac{\partial X}{\partial Z}) + \frac{\partial^2 X}{\partial Z^2}$
$+\sqrt{\Big[ +{\big[ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\big]}^2+ {\sigma(X)}^2 \big[Z^2 +3C_2 Z {(\frac{\partial Z}{\partial X})} \frac{\partial^2 X}{\partial Z^2} +C_4 {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3}\big] dt- C_3\frac{\partial \mu(X)}{\partial X} {(\frac{\partial X}{\partial Z})}^2) dt \Big]}$ Equation (21)

Please pardon any error with signs. I quickly wrote the post but I will carefully check it again and fix any errors in a day or two.

Friends, I redid some of the analysis of analytic solution of Fokker-Planck equation. I would further test it on computer and see how the existing algorithm needs to be improved. There might be some minor errors but I am presenting it in the hope that it would be useful. I will come with improvements in FPE solution algorithm code in a few days and post it here. Here is the analytic analysis I have redone.

I am starting with equation(13) of copied post.
$dX(Z{(\frac{\partial Z}{\partial X})}^2 + \frac{\partial^2 Z}{\partial X^2} )$
$= -\mu(X) (Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2}) dt$
$+2 \sigma(X) \frac{\partial \sigma(X)}{\partial X} (Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2}) dt$
$+ .5 {\sigma(X)}^2 ((Z^2-1) {(\frac{\partial Z}{\partial X})}^3 -3 Z \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} + \frac{\partial^3 Z}{\partial X^3}) dt$
$- \frac{\partial \mu(X)}{\partial X} |\frac{\partial Z}{\partial X}| dt$
$+({(\frac{\partial \sigma(X)}{\partial X})}^2+\sigma(X) \frac{\partial^2 \sigma(X)}{\partial X^2}) |\frac{\partial Z}{\partial X}| dt$  Equation(13)

For simplicity, I just consider the FPE of SDEs in Bessel coordinates with constant volatility. I am sure once this works for SDEs in Bessel coordinates, we would be able to extrapolate it to original coordinates. So we write the Above equation(13) of copied post for SDEs in Bessel coordinates as

$dX(Z{(\frac{\partial Z}{\partial X})}^2 + \frac{\partial^2 Z}{\partial X^2} )$
$= -\mu(X) (Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2}) dt$
$+ .5 {\sigma(X)}^2 ((Z^2-1) {(\frac{\partial Z}{\partial X})}^3 -3 Z \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} + \frac{\partial^3 Z}{\partial X^3}) dt$
$- \frac{\partial \mu(X)}{\partial X} |\frac{\partial Z}{\partial X}| dt$    Equation(13 A)

I re-write the equation as (I have changed the sign on second line below. It was probably wrong but I will fix this later in previous equations)

$dX(Z{(\frac{\partial Z}{\partial X})}^2 + \frac{\partial^2 Z}{\partial X^2} )$
$-\mu(X) (Z {(\frac{\partial Z}{\partial X})}^2+ \frac{\partial^2 Z}{\partial X^2}) dt$
$=+ .5 {\sigma(X)}^2 ((Z^2-1) {(\frac{\partial Z}{\partial X})}^3 -3 Z \frac{\partial Z}{\partial X} \frac{\partial^2 Z}{\partial X^2} + \frac{\partial^3 Z}{\partial X^3}) dt$
$- \frac{\partial \mu(X)}{\partial X} |\frac{\partial Z}{\partial X}| dt$    Equation(13 B)

I multiply the whole equation with ${(\frac{\partial X}{\partial Z})}^3$ and also re-arrange to get

$dX(Z (\frac{\partial X}{\partial Z}) + {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^2 Z}{\partial X^2} )$
$-\mu(X) (Z (\frac{\partial X}{\partial Z})+ {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^2 Z}{\partial X^2}) dt$
$- .5 {\sigma(X)}^2 dt$
$=+ .5 {\sigma(X)}^2 (Z^2 -3 Z {(\frac{\partial X}{\partial Z})}^2 \frac{\partial^2 Z}{\partial X^2} +{(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3}) dt$
$- \frac{\partial \mu(X)}{\partial X} {(\frac{\partial X}{\partial Z})}^2 dt$   Equation(14A)

We want to convert the first three lines of the above Equation(14 A) into a complete differential of square with following terms and we want to calculate the squared differential below so that we would know what terms have to be added/subtracted to first three lines of Equation(14 A) in order to make it a complete squared differential
$.5 d\Big[\int_{t_0}^{t} dX(s) - \int_{t_0}^{t} \mu(X) ds + ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]^2$  Equation(16A)

We apply the square inside the differential and get the following

$.5 d\Big[ {[\int_{t_0}^{t} dX]}^2 -2 [\int_{t_0}^{t} dX] [ \int_{t_0}^{t} \mu(X) ds] +2 [\int_{t_0}^{t} dX] ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})$
$+ {[\int_{t_0}^{t} \mu(X) ds]}^2 -2 [\int_{t_0}^{t} \mu(X) ds] ( Z (\frac{\partial X}{\partial Z})- \frac{\partial^2 X}{\partial Z^2})$
$+ {( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2}^2)}^2 \Big]$  Equation(16A)

We apply the differential on the above equation to get
$[\int_{t_0}^{t} dX(s)] dX(s)+ .5 {\sigma(x)}^2 ds- dX(s) [ \int_{t_0}^{t} \mu(X) ds]- [\int_{t_0}^{t} dX(s)] \mu(X) ds +dX(s) ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})$
$+ [\int_{t_0}^{t} \mu(X) ds] \mu(X) ds - \mu(X) ds ( Z (\frac{\partial X}{\partial Z})- \frac{\partial^2 X}{\partial Z^2})$  Equation(17A)

We only have second,  fifth and sixth term present in the equation coming from FPE equation. We will have to add rest of the terms on both sides and then complete the squared differential on LHS.
Considering first term, $][\int_{t_0}^{t} dX(s)] dX(s)$ we have
$[\int_{t_0}^{t} dX(s)] = \mu(X) (t-t_0) + \sigma(X) (Z(t)-Z(t_0))$
and
$dX(s)] = \mu(X) ds + \sigma(X) dZ(s)$
so we have the first terms in Equation(17A) as
$[\int_{t_0}^{t} dX(s)] dX(s)= {\mu(X)}^2 (t-t_0) dt + {\sigma(X)}^2 (Z(t)-Z(t_0))dZ(t) +\mu(X) \sigma(X) (t-t_0)dZ(t) +\mu(X) \sigma(X) (Z(t)-Z(t_0))d(t)$
and the third term in Equation(17A) as
$dX(s) [ \int_{t_0}^{t} \mu(X) ds]= {\mu(X)}^2 (t-t_0) dt + \mu(X) \sigma(X) (t-t_0) dZ(t)$\
and fourth term in Equation(17A) as
$[\int_{t_0}^{t} dX(s)] \mu(X) ds={\mu(X)}^2 (t-t_0) dt + +\mu(X) \sigma(X) (Z(t)-Z(t_0))d(t)$
and sixth term in Equation(17A) as
$[\int_{t_0}^{t} \mu(X) ds] \mu(X) ds={\mu(X)}^2 (t-t_0) dt$

So we could have from Equation() and Equation()
$.5 d\Big[\int_{t_0}^{t} dX(s) - \int_{t_0}^{t} \mu(X) ds + ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]^2$
$= +{\Big[ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]}^2$
$+ .5 {\sigma(X)}^2 (Z^2 +3C_2 Z {(\frac{\partial Z}{\partial X})} \frac{\partial^2 X}{\partial Z^2} +C_4 {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3}) dt$
$- C_3\frac{\partial \mu(X)}{\partial X} {(\frac{\partial X}{\partial Z})}^2 dt$
$+{\mu(X)}^2 (t-t_0) dt + {\sigma(X)}^2 (Z(t)-Z(t_0))dZ(t) +\mu(X) \sigma(X) (t-t_0)dZ(t) +\mu(X) \sigma(X) (Z(t)-Z(t_0))d(t)$
$-{\mu(X)}^2 (t-t_0) dt - \mu(X) \sigma(X) (t-t_0) dZ(t)$
$-{\mu(X)}^2 (t-t_0) dt - \mu(X) \sigma(X) (Z(t)-Z(t_0))d(t)$
$+{\mu(X)}^2 (t-t_0) dt$

After simplification, we have
$.5 d\Big[\int_{t_0}^{t} dX(s) - \int_{t_0}^{t} \mu(X) ds + ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]^2$
$= +{\Big[ ( Z (\frac{\partial X}{\partial Z}) - \frac{\partial^2 X}{\partial Z^2})\Big]}^2$
$+ .5 {\sigma(X)}^2 (Z^2 +3C_2 Z {(\frac{\partial Z}{\partial X})} \frac{\partial^2 X}{\partial Z^2} +C_4 {(\frac{\partial X}{\partial Z})}^3 \frac{\partial^3 Z}{\partial X^3}) dt$
$- C_3\frac{\partial \mu(X)}{\partial X} {(\frac{\partial X}{\partial Z})}^2 dt$
$+{\sigma(X)}^2 (.5 dH_2(t)-Z(t_0)dZ(t) )$

We can apply time integral on both sides and then take a square-root to solve the equation. The above analysis is very preliminary and only suggestive and may have some algebra errors. I will come with a detailed program in a few days.

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

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

I want to bring the attention of friends to this slightly old post that is very relevant especially when next president takes office.

viewtopic.php?f=15&t=94796&start=1155#p861398

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

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

Eleven days after receiving the antipsychotic injection(mixed with mind control drugs), I was feeling slightly better today after a very bad patch of more than ten days and had very little bouts of anxiety and uneasiness today.
My mother's heart surgery operation went well on Monday and she had to receive intensive care from hospital staff and our family and it will take her several months to recover properly.
In my research, I tried to back out volatility from the solution of FP equation so that the volatility could be used in the solution of SDE with transition probabilities. I made some progress but could not fully solve the problem. I also found a way to solve for covariance/correlation between SDE solution at two different stages(t1 and t2) in time. I will share it with friends after verifying it further. I also think that it would help us find the transition probabilities volatility analytically for diffusions with mean-reverting and more complicated drifts in the SDE and get us to solve for analytically beyond simple noises. I hope to do more experiments tomorrow and then come back to friends in a few days with the worked out programs.
I also want to request friends to approach the new president and his staff to end mind control of good, nice and talented individuals both in America and abroad. He has far more pressing things to take care at the start but I hope that plight of mind control victims also improves. I, for one, is on brutal mind control torture for past twenty two years and there are many others like me who have been on mind control for more than a decade. And I want to thank all those good and nice humans who have tried for my mind control to end in the past and I have been able to do some creative work in mathematics because my mind control decreased due to efforts from those great human beings. I will be indebted to all those good people for the rest of my life who tried to end my mind control. It is not possible for me to repay this debt other than to have regard and good wishes for those people who helped decrease my mind control.

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

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

I was able to code and verify the algorithm for calculation of Covariance between two arbitrary points on different time stages. I applied it in a model/SDE-free way to our underlying normal with unit variance increments and matched the covariances perfectly. Algorithm has nothing specific about unit variance increments normal and it can be applied/carried over to any SDE freely. In a day or two, I will be posting the coded algorithm with more details.
Covariance between two points $X_n(t_1)$ and $X_m(t_2)$ on the solution of any SDE can be calculated as
$=(X_n(t_1)-\mu_X(t_1)) P(X_n(t_1)) P_t(t_1,t_2,n,m) (X_m(t_2)-\mu_X(t_2))$
In our transition probabilities framework, we have that  $P(X_n(t_1))$ is the probability mass associated with grid point $X_n$ at time $t_1$ and $P_t(t_1,t_2,n,m)$ is the transition probability(In an algorithm, we can calculate it in string of single steps as well) between grid point $X_n$ at time $t_1$  and  grid point $X_m$ at arbitrary time $t_2$. For any SDE, these transition probabilities between two arbitrary points on different time grids can be calculated just once(on our unit variance increments normal) and then stored on disk and reused again and again. And we do not have to be in a transition probabilities framework to be able to calculate these covariances across time, we can in fact just as easily calculate these covariances in our Fokker-planck solution setting without adding any explicit transition probabilities framework.
When we talk about covariance across two time stages in usual setting, we usually do not mean covariance between two grid points on different time grids but we usually mean covariance averaged on entire time grids and we can easily do that in our framework.
Please give me another two or three days and I will be posting detailed programs to calculate these covariances across time for an SDE.

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

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

Friends, please pardon for the delay in posting the program to calculate covariances between SDE solution at two different dates. I have done this on transition probabilities grid but mean-reverting SDEs are still not working perfectly well on transition probabilities grid though the algorithm should work perfectly well once we get mean-reverting SDEs to work on this grid. For CEV type noises, it seems to work perfectly well but there is nothing very special about them. You might want to port the solution to FPE solution setting. If you consider it important, please mention here and I will port it to that setting as I will probably do at some later time anyway.
I tried to see why mean-reverting SDEs are not working in transition probabilities setting, and if the problem is in my calculation of drift? I used an alternative way and did the calculations with Girsanov theorem and multiplied the brownian noise(since in Bessel setting, all we have is simple brownian noise with drift) with Girsanov exponential derived from the SDE. And when I ran the program with Girsanov exponential, it gave me the same exact result that I was getting earlier with previous method. Therefore I realize that my transition probabilities algorithm requires more sophistication in updating the previous value of SDE when there is non-linear(mean-reverting type etc.) drift. Basic calculations of drift term are perfectly fine but algorithm for updating the previous values needs to be improved. Please give me a few days as I have some ideas about the solution that I mean to try. I have also included the calculations(and supporting code)  for Girsanov theorem update algorithm for friends who might be interested in looking at it.
Please note that you will have to caculate the standard normal transition increments and related probabilities from an older file and load them from disk as we did in previously posted program on this thread.

Here are the code files
function [] = SDETransProb08WmtGrid200b2()

%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
%Besse1 process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.
%When you run this program, please make sure that you look at graph zoomed
%at main body of the graph. Sometimes all the SDE is properly simulated but
%the tail explodes for some reason so you would know. These are minor
%issues that will be fixed in next versions of the program.
%I have put the terminal date at a very small value. Since calculations are
%done for all possible pairs of covariances on a very large number of time
%steps, increasing time to a large value is going to take quite a while.
%Start with a small date. You might want to modify this algorithm to
%calculate covariance pairs certain time spacing larger than minimum time spacing.
T=1/16;

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

dNn=.05;   % Normal density subdivisions width. would change with number of subdivisions
Nn=200;  % No of normal density subdivisions
NnMidl=100;%One half density Subdivision left from mid of normal density(low)
NnMidh=101;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;

Nn2=280;
dNn2=.05;
Nn2Midl=140;
Nn2Midh=141;
Nn2Mid=4.0;

x0=1.00;   % starting value of SDE
beta1=0.0;
beta2=1.00;   % Second drift term power.
gamma=.95;%50;   % volatility power.
kappa=0.0;%.950;   %mean reversion parameter.
theta=1.0;%mean reversion target
sigma0=.850;%Volatility value

%you can specify any general mu1 and mu2 and beta1 and beta2.
mu1=+1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
%mu1=0;
%mu2=0;

alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;%This is for expansion of integrals for calculation of drift
%and volatility coefficients

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

%Z(1:Nn)=(((1:Nn)-5.5)*dNn-NnMid);
%Z(1:Nn)=(((1:Nn)-20.5)*dNn-NnMid);
Z(1:Nn)=(((1:Nn)-20.5)*dNn-NnMid);

%Z(1)=-5.118453;
%Z(Nn)=5.118453;

Z2(1:Nn2)=(((1:Nn2)-60.5)*dNn2-Nn2Mid);

Z
Z2
str=input('Look at Zs')
%ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
%ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);
%ZProb(2:Nn-1)=normcdf(.5*Z(2:Nn-1)+.5*Z(3:Nn),0,1)-normcdf(.5*Z(2:Nn-1)+.5*Z(1:Nn-2),0,1);
%Above calculate probability mass in each probability subdivision.

ZProb(1)=normcdf(-4.95);
ZProb(Nn)=1-normcdf(4.95);
ZProb(2)=normcdf(-4.9)-normcdf(-4.95);
ZProb(Nn-1)=normcdf(4.95)-normcdf(4.9);
ZProb(3:Nn-2)=normcdf(.5*Z(3:Nn-2)+.5*Z(4:Nn-1),0,1)-normcdf(.5*Z(3:Nn-2)+.5*Z(2:Nn-3),0,1);

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

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

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

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

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

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Current program has 50 grid points. I have divided every grid subdivision
%into three further subdivisions with probabilities P0,Pb,Pa.

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

Pprev(1:Nn)=ZProb(1:Nn);%probability mass at starting node.
Pnew(1:Nn2)=0.0;
EZwnew(1:Nn2)=0.0;
EZwprev(1:Nn)=0.0;
Ewnew(1:Nn2)=0;
Ewprev(1:Nn)=0;
Ew1new(1:Nn2)=0;
Ew1prev(1:Nn)=0;

CovarZt(1:Tt,1:Tt)=0.0;
Covarwt(1:Tt,1:Tt)=0.0;
%ECovZnew0(1:Tt,1:Tt,1:Nn2)=0.0;
ECovZprev0(1:Tt,1:Tt,1:Nn)=0.0;
%ECovwnew0(1:Tt,1:Tt,1:Nn2)=0.0;
ECovwprev0(1:Tt,1:Tt,1:Nn)=0.0;

VarZ(1:Tt)=0;
Varw(1:Tt)=0;
ValZprev(1:Nn)=0;
Valwprev(1:Nn)=0;
ValZnew(1:Nn2)=0;
Valwnew(1:Nn2)=0;
ValZnow(1:Nn)=0;
Valwnow(1:Nn)=0;
Zwt0(1:Nn)=0;
wt0(1:Nn)=0;

tic

for tt=1:Tt
t2=tt*dt;
t1=(tt-1)*dt;
if(tt==1)

[wMu0dt,c1,c22] = CalculateDriftAndVolA8Trans(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
Pnew(1:Nn2)=pt0(1:Nn2);
EZwnew(1:Nn2)=pt0(1:Nn2).*Zt0(1:Nn2);
Ewnew(1:Nn2)=(x0^(1-gamma)/(1-gamma)+wMu0dt(1)+c1(1).*Zt0(1:Nn2)+c22(1).*(Zt0(1:Nn2).^2-1)).*pt0(1:Nn2);
Ew1new(1:Nn2)=(x0^(1-gamma)/(1-gamma)+wMu0dt(1)+c1(1).*Zt0(1:Nn2)+c22(1).*(Zt0(1:Nn2).^2-1)).*pt0(1:Nn2);

end
if(tt>1)
%Convert the expected values in each cell to standard Values associated with the grid cell after removing the
%integrated probability in the cell.
w(1:Nn)=Ewprev(1:Nn)./Pprev(1:Nn);
w1(1:Nn)=Ew1prev(1:Nn)./Pprev(1:Nn);
[wMu0dt,wMu1dt,dwMu0dtdw,d2wMu0dtdw2,c1] = CalculateDriftAndVolA404(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);

[w1Mu0dt,w1Mu1dt,dw1Mu0dtdw1,d2w1Mu0dtdw2,cw1] = CalculateDriftAndVolA404(w1,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);

pt(1:Nn,1:Nn2)=ptt(tt,1:Nn,1:Nn2);
Zt(1:Nn,1:Nn2)=Ztt(tt,1:Nn,1:Nn2);
%#Girsanov
%For solution of drifted SDEs with Girsanov,uncomment the next three lines.
%       wmStart=1;
%       Mm=Nn2;
%       [GirsanovExp] = CalculateGirsanovExponential(wnStart,Nn,wmStart,Mm,w,Zt,dt,mu1,mu2,beta1,beta2,gamma,sigma0);
for mm=1:Nn

%#Girsanov
%Uncomment the next line to see that expected value of Girsanov exponential
%should be unity.
%           sum(GirsanovExp(mm,1:Nn2).*pt(mm,1:Nn2))

A0=ptA(mm,tt);
B0=ptB(mm,tt);
%A0 and B0 are starting and ending boundaries for non-zero
%elements in larger target array. These are boundaries for
%transition probabilities originating from cell mm on first grid
% at time tt. Out of the 280 target grid
%cell, only a few ten are non-zero. These non-zero elements
%boundaries(A0 and B0) were calculated in a different module and saved on
%disk just when we calculated the transition probabilities.
%Later we retrieved them from disk to use in this program.

Pnew(A0:B0)=Pnew(A0:B0)+Pprev(mm).*pt(mm,A0:B0);
EZwnew(A0:B0)=EZwnew(A0:B0)+EZwprev(mm).*pt(mm,A0:B0)+ ...
Zt(mm,A0:B0).*pt(mm,A0:B0).*Pprev(mm);

%Below is the transition probabilities simulation
%with improvised term. These improvised terms make an improvement
%over the base model which is given in next paragraph of code.
Ewnew(A0:B0)=Ewnew(A0:B0)+Ewprev(mm).*pt(mm,A0:B0)+ ...
(wMu0dt(mm)+c1(mm).*Zt(mm,A0:B0)).* ...
pt(mm,A0:B0).*Pprev(mm);

%#Girsanov
%For Girsanov, uncomment the next block and comment the previous block.
%            Ewnew(A0:B0)=Ewnew(A0:B0)+Ewprev(mm).*pt(mm,A0:B0)+ ...
%                GirsanovExp(mm,A0:B0).*sigma0.*sqrt(dt).*Zt(mm,A0:B0).* ...
%                pt(mm,A0:B0).*Pprev(mm);

Ewnew(Ewnew<0)=0;
Ewnew(isnan(Ewnew)==1)=0;
%Below is the base transition probabilities simulation
%with same terms as in monte carlo simulation. cw1 term is analogue
%of c1 which is coefficient of Zt in our bessel form SDE.
Ew1new(A0:B0)=Ew1new(A0:B0)+Ew1prev(mm).*pt(mm,A0:B0)+ ...
(w1Mu0dt(mm)+cw1(mm).*Zt(mm,A0:B0)).* ...
pt(mm,A0:B0).*Pprev(mm);

Ew1new(Ew1new<0)=0;
Ew1new(isnan(Ew1new)==1)=0;

end

wa(1:Nn-1)=w(1:Nn-1);
wa(Nn)=w(Nn);
wb(1:Nn-1)=w(2:Nn);
wb(Nn)=w(Nn-1)+w(Nn-1)-w(Nn-2);
w(wa(:)>wb(:))=wb(wa(:)>wb(:));%Be careful;might not universally hold;
%if(w(Nn)>(2*w(Nn-1)-w(Nn-2)))
%    w(Nn)=(2*w(Nn-1)-w(Nn-2));
%end

end

Pprev(1)=sum(Pnew(1:41));
Pprev(2:Nn-1)=Pnew(42:40+Nn-1);
Pprev(Nn)=sum(Pnew(40+Nn:80+Nn));%All the probability mass in last 41
%subdivisions on target grid is assigned to end boundary point on
%the grid that will become smaller first grid for next simulation step.

EZwprev(1)=sum(EZwnew(1:41));
EZwprev(2:Nn-1)=EZwnew(42:40+Nn-1);
EZwprev(Nn)=sum(EZwnew(40+Nn:80+Nn));

Zwt(tt,1:Nn)=EZwprev(1:Nn)./Pprev(1:Nn);
MuZ(tt)=sum(EZwprev(1:Nn));

Ewprev(1)=sum(Ewnew(1:41));
Ewprev(2:Nn-1)=Ewnew(42:40+Nn-1);
Ewprev(Nn)=sum(Ewnew(40+Nn:80+Nn));

wt(tt,1:Nn)=Ewprev(1:Nn)./Pprev(1:Nn);
Muw(tt)=sum(Ewprev(1:Nn));

Ew1prev(1)=sum(Ew1new(1:41));
Ew1prev(2:Nn-1)=Ew1new(42:40+Nn-1);
Ew1prev(Nn)=sum(Ew1new(40+Nn:80+Nn));

Pnew(1:Nn2)=0;
EZwnew(1:Nn2)=0.0;
Ewnew(1:Nn2)=0.0;
Ew1new(1:Nn2)=0.0;

%#Covariances Calculation.
%Below, we calculate covariances between the each date with all of the
%later following dates

if(tt==1)
VarZ(tt)=sum((Zwt(tt,1:Nn)-MuZ(tt)).^2.*ZProb(1:Nn));
Varw(tt)=sum((wt(tt,1:Nn)-Muw(tt)).^2.*ZProb(1:Nn));
end
if(tt>1)
VarZ(tt)=sum((Zwt(tt,1:Nn)-MuZ(tt)).^2.*ZProb(1:Nn));
Varw(tt)=sum((wt(tt,1:Nn)-Muw(tt)).^2.*ZProb(1:Nn));

%Below loop with variable tt1 indicates the starting time for various
%covariances
%Below loop variable tt2 indicates the algorithm calculations on dates
%following the strting time.
for tt1=1:tt-1
for tt2=tt:Tt
tt1
tt2
if(tt1<tt-1)
%If the starting time for covariance calculations is smaller than one date
%ago, it means that we do not have to initialize the algorithm and we are
%in the middle of the calculations.
ValZprev(1:Nn)=ECovZprev0(tt1,tt2,1:Nn);
Valwprev(1:Nn)=ECovwprev0(tt1,tt2,1:Nn);
end
if(tt1==tt-1)
%If the starting time for covariance calculations is equal to one date
%ago, it means that we haveto initialize the algorithm and we are
%starting the calculations. Here we assign the starting values for covariance
%calculation.
Zwt0(1:Nn)=Zwt(tt-1,1:Nn);
ValZprev(1:Nn)=(Zwt0(1:Nn)-MuZ(tt-1)).*ZProb(1:Nn);
wt0(1:Nn)=wt(tt-1,1:Nn);
Valwprev(1:Nn)=(wt0(1:Nn)-Muw(tt-1)).*ZProb(1:Nn);
end
ValZnew(1:Nn2)=0;
Valwnew(1:Nn2)=0;
%Below we calculate the running string for the transition probability increment
%between first covariance date and last covariance date.
for mm=1:Nn
A0=ptA(mm,tt);
B0=ptB(mm,tt);
pt(mm,1:Nn2)=ptt(tt,mm,1:Nn2);
ValZnew(A0:B0)=ValZnew(A0:B0)+ ...
ValZprev(mm).*pt(mm,A0:B0);
Valwnew(A0:B0)=Valwnew(A0:B0)+ ...
Valwprev(mm).*pt(mm,A0:B0);
end
%ECovZnew0(tt1,tt2,1:Nn2)=ValZnew(1:Nn2);
%ECovwnew0(tt1,tt2,1:Nn2)=Valwnew(1:Nn2);
ECovZprev0(tt1,tt2,1)=sum(ValZnew(1:41));
ECovZprev0(tt1,tt2,2:Nn-1)=ValZnew(42:40+Nn-1);
ECovZprev0(tt1,tt2,Nn)=sum(ValZnew(40+Nn:80+Nn));

ECovwprev0(tt1,tt2,1)=sum(Valwnew(1:41));
ECovwprev0(tt1,tt2,2:Nn-1)=Valwnew(42:40+Nn-1);
ECovwprev0(tt1,tt2,Nn)=sum(Valwnew(40+Nn:80+Nn));
end
%Below when second covariance pair date is equal to present date, we do
%final(ending) calculations for covariance calculation with all first pairs from
%all of the earlier dates.
ValZnow(1:Nn)=Zwt(tt,1:Nn)-MuZ(tt);
Valwnow(1:Nn)=wt(tt,1:Nn)-Muw(tt);
ValZprev(1:Nn)=ECovZprev0(tt1,tt,1:Nn);
Valwprev(1:Nn)=ECovwprev0(tt1,tt,1:Nn);
CovarZt(tt1,tt)=sum(ValZprev(1:Nn).*ValZnow(1:Nn));
Covarwt(tt1,tt)=sum(Valwprev(1:Nn).*Valwnow(1:Nn));

end
end
end

%At the end, we want to see the results for our covariance calulations
for tt=1:Tt-1
%Below, we see the variance of standard increments normal on a particular date
VarZ(tt)
%Now, we want to see the covariance of all the following dates with present date
%and for normal noise this should be equal to the variance at starting
%date.
CovarZt(tt,tt+1:Tt)
tt
str=input('Look at CovarZt');
%Below, we see the variance of a general SDE(but this current program works well
%with driftless noises only) on a particular date
Varw(tt)
%Now, we want to see the covariance of all the following dates with present date
%and for CEV type driftless SDEs, this should be equal to the variance at starting
%date.
Covarwt(tt,tt+1:Tt)
tt
str=input('Look at Covarwt');
end

plot((1:Nn),ZProb(1:Nn),'g',(1:Nn),Pprev(1:Nn),'r');
%plot((1:Nn),ZProb(1:Nn),'g',(1:Nn),Pprev(1:Nn),'r');
%str=input('Look at distributions');

Pnew=Pprev;
%Convert the expected values in each cell to standard Values associated with the grid cell after removing the
%integrated probability in the cell.
w(wnStart:Nn)=Ewprev(wnStart:Nn)./Pprev(wnStart:Nn);
w1(wnStart:Nn)=Ew1prev(wnStart:Nn)./Pprev(wnStart:Nn);
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
yy1(wnStart:Nn)=((1-gamma)*w1(wnStart:Nn)).^(1/(1-gamma));

%below D's (the names of variables starting with D) are
%change of probability derivatives.

y_w(1:Nn)=0;
y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
y_w1(1:Nn)=0;
y_w1(wnStart:Nn) = ((1-gamma)*w1(wnStart:Nn)).^(1/(1-gamma));
Dfy_w(wnStart:Nn)=0;
Dfy_w1(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));
Dfy_w1(nn) = (y_w1(nn + 1) - y_w1(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
py_w(1:Nn)=0;
py_w1(1:Nn)=0;
for nn = wnStart:Nn-1
py_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
py_w1(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w1(nn));
end

toc

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

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

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

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

end

YY(YY<0)=0;
disp('Original process average from monte carlo');
MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
MCVar=sum((YY(:)-MCMean).^2)/paths
disp('Original process average from our simulation');
ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates
ItoHermiteVar=sum((y_w(wnStart+1:Nn-1)-ItoHermiteMean).^2.*ZProb(wnStart+1:Nn-1))
disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average

MaxCutOff=30;
NoOfBins=round(300*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );

%plot(y_w1(wnStart+1:Nn-1),py_w1(wnStart+1:Nn-1),'b',y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
plot(y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',y_w1(wnStart+1:Nn-1),py_w1(wnStart+1:Nn-1),'b',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');

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

legend({'Improvised terms Model Density','Monte Carlo Density'},'Location','northeast')

str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
end
.
Here is the Girsanov exponential calculation function file.
.
function [GirsanovExp] = CalculateGirsanovExponential(wnStart,Nn,wmStart,Mm,w,Zt,dt,mu1,mu2,beta1,beta2,gamma,sigma0)

%muG is simple drift expression.
muG(wnStart:Nn)=mu1.*((1-gamma).*w(wnStart:Nn)).^((beta1-gamma)/(1-gamma))+ ...
mu2.*((1-gamma).*w(wnStart:Nn)).^((beta2-gamma)/(1-gamma))- ...
.5*sigma0^2*gamma.*((1-gamma).*w(wnStart:Nn)).^(-1);

dmuGdw(wnStart:Nn)=(1-gamma).*(mu1.*((beta1-gamma)/(1-gamma)).*((1-gamma).*w(wnStart:Nn)).^((beta1-gamma)/(1-gamma)-1)+ ...
mu2.*((beta2-gamma)/(1-gamma)).*((1-gamma).*w(wnStart:Nn)).^((beta2-gamma)/(1-gamma)-1)- ...
.5*sigma0^2*gamma.*(-1).*((1-gamma).*w(wnStart:Nn)).^(-2));

d2muGdw2(wnStart:Nn)=(1-gamma)^2.*(mu1.*((beta1-gamma)/(1-gamma)).*((beta1-gamma)/(1-gamma)-1).*((1-gamma).*w(wnStart:Nn)).^((beta1-gamma)/(1-gamma)-2)+ ...
mu2.*((beta2-gamma)/(1-gamma)).*((beta2-gamma)/(1-gamma)).*((1-gamma).*w(wnStart:Nn)).^((beta2-gamma)/(1-gamma)-2)- ...
.5*sigma0^2*gamma.*(-1).*(-2).*((1-gamma).*w(wnStart:Nn)).^(-3));

%mu2G is squared drift expression.
mu2G(wnStart:Nn)=mu1^2.*((1-gamma).*w(wnStart:Nn)).^(2*(beta1-gamma)/(1-gamma))+ ...
mu2^2.*((1-gamma).*w(wnStart:Nn)).^(2*(beta2-gamma)/(1-gamma))+ ...
.25*sigma0^4*gamma^2.*((1-gamma).*w(wnStart:Nn)).^(-2)+ ...
2*mu1*mu2*((1-gamma).*w(wnStart:Nn)).^((beta1-gamma)/(1-gamma)+(beta2-gamma)/(1-gamma))- ...
2*mu1*.5*gamma*sigma0^2.*((1-gamma)*w(wnStart:Nn)).^((beta1-gamma)/(1-gamma)-1)- ...
2*mu2*.5*gamma*sigma0^2.*((1-gamma)*w(wnStart:Nn)).^((beta2-gamma)/(1-gamma)-1);

dmu2Gdw(wnStart:Nn)=(1-gamma).*(mu1^2.*(2*(beta1-gamma)/(1-gamma)).*((1-gamma).*w(wnStart:Nn)).^(2*(beta1-gamma)/(1-gamma)-1)+ ...
mu2^2.*(2*(beta2-gamma)/(1-gamma)).*((1-gamma).*w(wnStart:Nn)).^(2*(beta2-gamma)/(1-gamma)-1)+ ...
.25*sigma0^4*gamma^2.*(-2).*((1-gamma).*w(wnStart:Nn)).^(-3)+ ...
2*mu1*mu2*((beta1-gamma)/(1-gamma)+(beta2-gamma)/(1-gamma)).*((1-gamma).*w(wnStart:Nn)).^((beta1-gamma)/(1-gamma)+(beta2-gamma)/(1-gamma)-1)- ...
2*mu1*.5*gamma*sigma0^2.*((beta1-gamma)/(1-gamma)-1).*((1-gamma)*w(wnStart:Nn)).^((beta1-gamma)/(1-gamma)-2)- ...
2*mu2*.5*gamma*sigma0^2.*((beta2-gamma)/(1-gamma)-1).*((1-gamma)*w(wnStart:Nn)).^((beta2-gamma)/(1-gamma)-2));

d2mu2Gdw2(wnStart:Nn)=(1-gamma)^2.*(mu1^2.*(2*(beta1-gamma)/(1-gamma)).*(2*(beta1-gamma)/(1-gamma)-1).* ...
((1-gamma).*w(wnStart:Nn)).^(2*(beta1-gamma)/(1-gamma)-2)+ ...
mu2^2.*(2*(beta2-gamma)/(1-gamma)).*(2*(beta2-gamma)/(1-gamma)-1).*((1-gamma).*w(wnStart:Nn)).^(2*(beta2-gamma)/(1-gamma)-2)+ ...
.25*sigma0^4*gamma^2.*(-2).*(-3).*((1-gamma).*w(wnStart:Nn)).^(-4)+ ...
2*mu1*mu2*((beta1-gamma)/(1-gamma)+(beta2-gamma)/(1-gamma)).*((1-gamma).*((beta1-gamma)/(1-gamma)+(beta2-gamma)/(1-gamma)-1).* ...
w(wnStart:Nn)).^((beta1-gamma)/(1-gamma)+(beta2-gamma)/(1-gamma)-2)- ...
2*mu1*.5*gamma*sigma0^2.*((beta1-gamma)/(1-gamma)-1).*((beta1-gamma)/(1-gamma)-2).*((1-gamma)*w(wnStart:Nn)).^((beta1-gamma)/(1-gamma)-3)- ...
2*mu2*.5*gamma*sigma0^2.*((beta2-gamma)/(1-gamma)-1).*((beta2-gamma)/(1-gamma)-2).*((1-gamma)*w(wnStart:Nn)).^((beta2-gamma)/(1-gamma)-3));

GTermZ(wnStart:Nn,wmStart:Mm)=0; %This is dz part of Girsanov exponential
GTermt(wnStart:Nn,wmStart:Mm)=0;%This is dt part of Girsanov exponential
GirsanovExp(wnStart:Nn,wmStart:Mm)=0;
for nn=wnStart:Nn
GTermZ(nn,wmStart:Mm)=muG(nn)./sigma0.*Zt(nn,wmStart:Mm)*sqrt(dt)+ ...
dmuGdw(nn)./sigma0.*sigma0.*(Zt(nn,wmStart:Mm).^2-1)/(2)*dt+ ...
.5*d2muGdw2(nn)./sigma0.*sigma0^2.*1/sqrt(3).*Zt(nn,wmStart:Mm).*dt^1.5;
GTermt(nn,wmStart:Mm)=mu2G(nn)./sigma0^2.*dt+ ...
dmu2Gdw(nn)./sigma0^2.*sigma0.*(1-1/sqrt(3)).*Zt(nn,wmStart:Mm).*dt^1.5+ ...
.5*d2mu2Gdw2(nn)./sigma0^2.*sigma0^2*dt^2/2;

GirsanovExp(nn,wmStart:Mm)=1+GTermZ(nn,wmStart:Mm)-.5*GTermt(nn,wmStart:Mm)+ ...
.5*GTermZ(nn,wmStart:Mm).^2+.5*.25*GTermt(nn,wmStart:Mm).^2- ...
.5*2*GTermZ(nn,wmStart:Mm)*.5.*GTermt(nn,wmStart:Mm);

end

end


Here are the final few lines you will see on your screen when you run the first function.
MCMean =
0.999075562104895
MCVar =
0.045976768541907
Original process average from our simulation
ItoHermiteMean =
1.000003155540884
ItoHermiteVar =
0.046045848931585
true Mean only applicble to standard SV mean reverting type models otherwise disregard
TrueMean =
1
IndexMax =
922

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

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

Sorry Friends, the algorithm I posted yesterday to calculate variances had an extra unnecessary nested loop. I have removed the redundant loop and now the algorithm is way faster and you can far easily use it for calculation of covariances on longer time horizons.
Here is the new modified efficient program in place of program posted yesterday.
function [] = SDETransProb08WmtGrid200b2()

%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
%Besse1 process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.
%When you run this program, please make sure that you look at graph zoomed
%at main body of the graph. Sometimes all the SDE is properly simulated but
%the tail explodes for some reason so you would know. These are minor
%issues that will be fixed in next versions of the program.
%I have put the terminal date at a very small value. Since calculations are
%done for all possible pairs of covariances on a very large number of time
%steps, increasing time to a large value is going to take quite a while.
%Start with a small date. You might want to modify this algorithm to
%calculate covariance pairs certain time spacing larger than minimum time spacing.
T=1/4;

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

dNn=.05;   % Normal density subdivisions width. would change with number of subdivisions
Nn=200;  % No of normal density subdivisions
NnMidl=100;%One half density Subdivision left from mid of normal density(low)
NnMidh=101;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;

Nn2=280;
dNn2=.05;
Nn2Midl=140;
Nn2Midh=141;
Nn2Mid=4.0;

x0=1.00;   % starting value of SDE
beta1=0.0;
beta2=1.00;   % Second drift term power.
gamma=.95;%50;   % volatility power.
kappa=0.0;%.950;   %mean reversion parameter.
theta=1.0;%mean reversion target
sigma0=.850;%Volatility value

%you can specify any general mu1 and mu2 and beta1 and beta2.
mu1=+1*theta*kappa;   %first drift coefficient.
mu2=-1*kappa;    % Second drift coefficient.
%mu1=0;
%mu2=0;

alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;%This is for expansion of integrals for calculation of drift
%and volatility coefficients

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

%Z(1:Nn)=(((1:Nn)-5.5)*dNn-NnMid);
%Z(1:Nn)=(((1:Nn)-20.5)*dNn-NnMid);
Z(1:Nn)=(((1:Nn)-20.5)*dNn-NnMid);

%Z(1)=-5.118453;
%Z(Nn)=5.118453;

Z2(1:Nn2)=(((1:Nn2)-60.5)*dNn2-Nn2Mid);

Z
Z2
str=input('Look at Zs')
%ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
%ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);
%ZProb(2:Nn-1)=normcdf(.5*Z(2:Nn-1)+.5*Z(3:Nn),0,1)-normcdf(.5*Z(2:Nn-1)+.5*Z(1:Nn-2),0,1);
%Above calculate probability mass in each probability subdivision.

ZProb(1)=normcdf(-4.95);
ZProb(Nn)=1-normcdf(4.95);
ZProb(2)=normcdf(-4.9)-normcdf(-4.95);
ZProb(Nn-1)=normcdf(4.95)-normcdf(4.9);
ZProb(3:Nn-2)=normcdf(.5*Z(3:Nn-2)+.5*Z(4:Nn-1),0,1)-normcdf(.5*Z(3:Nn-2)+.5*Z(2:Nn-3),0,1);

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

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

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

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

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

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Current program has 50 grid points. I have divided every grid subdivision
%into three further subdivisions with probabilities P0,Pb,Pa.

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

Pprev(1:Nn)=ZProb(1:Nn);%probability mass at starting node.
Pnew(1:Nn2)=0.0;
EZwnew(1:Nn2)=0.0;
EZwprev(1:Nn)=0.0;
Ewnew(1:Nn2)=0;
Ewprev(1:Nn)=0;
Ew1new(1:Nn2)=0;
Ew1prev(1:Nn)=0;

CovarZt(1:Tt,1:Tt)=0.0;
Covarwt(1:Tt,1:Tt)=0.0;
%ECovZnew0(1:Tt,1:Tt,1:Nn2)=0.0;
ECovZprev0(1:Tt,1:Nn)=0.0;
%ECovwnew0(1:Tt,1:Tt,1:Nn2)=0.0;
ECovwprev0(1:Tt,1:Nn)=0.0;

VarZ(1:Tt)=0;
Varw(1:Tt)=0;
ValZprev(1:Nn)=0;
Valwprev(1:Nn)=0;
ValZnew(1:Nn2)=0;
Valwnew(1:Nn2)=0;
ValZnow(1:Nn)=0;
Valwnow(1:Nn)=0;
Zwt0(1:Nn)=0;
wt0(1:Nn)=0;

tic

for tt=1:Tt
t2=tt*dt;
t1=(tt-1)*dt;
if(tt==1)

[wMu0dt,c1,c22] = CalculateDriftAndVolA8Trans(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
Pnew(1:Nn2)=pt0(1:Nn2);
EZwnew(1:Nn2)=pt0(1:Nn2).*Zt0(1:Nn2);
Ewnew(1:Nn2)=(x0^(1-gamma)/(1-gamma)+wMu0dt(1)+c1(1).*Zt0(1:Nn2)+c22(1).*(Zt0(1:Nn2).^2-1)).*pt0(1:Nn2);
Ew1new(1:Nn2)=(x0^(1-gamma)/(1-gamma)+wMu0dt(1)+c1(1).*Zt0(1:Nn2)+c22(1).*(Zt0(1:Nn2).^2-1)).*pt0(1:Nn2);

end
if(tt>1)
%Convert the expected values in each cell to standard Values associated with the grid cell after removing the
%integrated probability in the cell.
w(1:Nn)=Ewprev(1:Nn)./Pprev(1:Nn);
w1(1:Nn)=Ew1prev(1:Nn)./Pprev(1:Nn);
[wMu0dt,wMu1dt,dwMu0dtdw,d2wMu0dtdw2,c1] = CalculateDriftAndVolA404(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);

[w1Mu0dt,w1Mu1dt,dw1Mu0dtdw1,d2w1Mu0dtdw2,cw1] = CalculateDriftAndVolA404(w1,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);

pt(1:Nn,1:Nn2)=ptt(tt,1:Nn,1:Nn2);
Zt(1:Nn,1:Nn2)=Ztt(tt,1:Nn,1:Nn2);
%#Girsanov
%For solution of drifted SDEs with Girsanov,uncomment the next three lines.
%       wmStart=1;
%       Mm=Nn2;
%       [GirsanovExp] = CalculateGirsanovExponential(wnStart,Nn,wmStart,Mm,w,Zt,dt,mu1,mu2,beta1,beta2,gamma,sigma0);
for mm=1:Nn

%#Girsanov
%Uncomment the next line to see that expected value of Girsanov exponential
%should be unity.
%           sum(GirsanovExp(mm,1:Nn2).*pt(mm,1:Nn2))

A0=ptA(mm,tt);
B0=ptB(mm,tt);
%A0 and B0 are starting and ending boundaries for non-zero
%elements in larger target array. These are boundaries for
%transition probabilities originating from cell mm on first grid
% at time tt. Out of the 280 target grid
%cell, only a few ten are non-zero. These non-zero elements
%boundaries(A0 and B0) were calculated in a different module and saved on
%disk just when we calculated the transition probabilities.
%Later we retrieved them from disk to use in this program.

Pnew(A0:B0)=Pnew(A0:B0)+Pprev(mm).*pt(mm,A0:B0);
EZwnew(A0:B0)=EZwnew(A0:B0)+EZwprev(mm).*pt(mm,A0:B0)+ ...
Zt(mm,A0:B0).*pt(mm,A0:B0).*Pprev(mm);

%Below is the transition probabilities simulation
%with improvised term. These improvised terms make an improvement
%over the base model which is given in next paragraph of code.
Ewnew(A0:B0)=Ewnew(A0:B0)+Ewprev(mm).*pt(mm,A0:B0)+ ...
(wMu0dt(mm)+c1(mm).*Zt(mm,A0:B0)).* ...
pt(mm,A0:B0).*Pprev(mm);

%#Girsanov
%For Girsanov, uncomment the next block and comment the previous block.
%            Ewnew(A0:B0)=Ewnew(A0:B0)+Ewprev(mm).*pt(mm,A0:B0)+ ...
%                GirsanovExp(mm,A0:B0).*sigma0.*sqrt(dt).*Zt(mm,A0:B0).* ...
%                pt(mm,A0:B0).*Pprev(mm);

Ewnew(Ewnew<0)=0;
Ewnew(isnan(Ewnew)==1)=0;
%Below is the base transition probabilities simulation
%with same terms as in monte carlo simulation. cw1 term is analogue
%of c1 which is coefficient of Zt in our bessel form SDE.
Ew1new(A0:B0)=Ew1new(A0:B0)+Ew1prev(mm).*pt(mm,A0:B0)+ ...
(w1Mu0dt(mm)+cw1(mm).*Zt(mm,A0:B0)).* ...
pt(mm,A0:B0).*Pprev(mm);

Ew1new(Ew1new<0)=0;
Ew1new(isnan(Ew1new)==1)=0;

end

wa(1:Nn-1)=w(1:Nn-1);
wa(Nn)=w(Nn);
wb(1:Nn-1)=w(2:Nn);
wb(Nn)=w(Nn-1)+w(Nn-1)-w(Nn-2);
w(wa(:)>wb(:))=wb(wa(:)>wb(:));%Be careful;might not universally hold;
%if(w(Nn)>(2*w(Nn-1)-w(Nn-2)))
%    w(Nn)=(2*w(Nn-1)-w(Nn-2));
%end

end

Pprev(1)=sum(Pnew(1:41));
Pprev(2:Nn-1)=Pnew(42:40+Nn-1);
Pprev(Nn)=sum(Pnew(40+Nn:80+Nn));%All the probability mass in last 41
%subdivisions on target grid is assigned to end boundary point on
%the grid that will become smaller first grid for next simulation step.

EZwprev(1)=sum(EZwnew(1:41));
EZwprev(2:Nn-1)=EZwnew(42:40+Nn-1);
EZwprev(Nn)=sum(EZwnew(40+Nn:80+Nn));

Zwt(tt,1:Nn)=EZwprev(1:Nn)./Pprev(1:Nn);
MuZ(tt)=sum(EZwprev(1:Nn));

Ewprev(1)=sum(Ewnew(1:41));
Ewprev(2:Nn-1)=Ewnew(42:40+Nn-1);
Ewprev(Nn)=sum(Ewnew(40+Nn:80+Nn));

wt(tt,1:Nn)=Ewprev(1:Nn)./Pprev(1:Nn);
Muw(tt)=sum(Ewprev(1:Nn));

Ew1prev(1)=sum(Ew1new(1:41));
Ew1prev(2:Nn-1)=Ew1new(42:40+Nn-1);
Ew1prev(Nn)=sum(Ew1new(40+Nn:80+Nn));

Pnew(1:Nn2)=0;
EZwnew(1:Nn2)=0.0;
Ewnew(1:Nn2)=0.0;
Ew1new(1:Nn2)=0.0;
%#Covariances Calculation.
%Below, we calculate covariances between the each date with all of the
%later following dates

if(tt==1)
VarZ(tt)=sum((Zwt(tt,1:Nn)-MuZ(tt)).^2.*ZProb(1:Nn));
Varw(tt)=sum((wt(tt,1:Nn)-Muw(tt)).^2.*ZProb(1:Nn));
end
if(tt>1)
VarZ(tt)=sum((Zwt(tt,1:Nn)-MuZ(tt)).^2.*ZProb(1:Nn));
Varw(tt)=sum((wt(tt,1:Nn)-Muw(tt)).^2.*ZProb(1:Nn));

%Below loop with variable tt1 indicates the starting time for various
%covariances
for tt1=1:tt-1
tt1
tt
if(tt1<tt-1)
%If the starting time for covariance calculations is smaller than one date
%ago, it means that we do not have to initialize the algorithm and we are
%in the middle of the calculations.
ValZprev(1:Nn)=ECovZprev0(tt1,1:Nn);
Valwprev(1:Nn)=ECovwprev0(tt1,1:Nn);
end
if(tt1==tt-1)
%If the starting time for covariance calculations is equal to one date
%ago, it means that we haveto initialize the algorithm and we are
%starting the calculations. Here we assign the starting values for covariance
%calculation.
Zwt0(1:Nn)=Zwt(tt-1,1:Nn);
ValZprev(1:Nn)=(Zwt0(1:Nn)-MuZ(tt-1)).*ZProb(1:Nn);
wt0(1:Nn)=wt(tt-1,1:Nn);
Valwprev(1:Nn)=(wt0(1:Nn)-Muw(tt-1)).*ZProb(1:Nn);
end
ValZnew(1:Nn2)=0;
Valwnew(1:Nn2)=0;
%Below we calculate the running string for the transition probability increment
%between first covariance date and last covariance date.
for mm=1:Nn
A0=ptA(mm,tt);
B0=ptB(mm,tt);
pt(mm,1:Nn2)=ptt(tt,mm,1:Nn2);
ValZnew(A0:B0)=ValZnew(A0:B0)+ ...
ValZprev(mm).*pt(mm,A0:B0);
Valwnew(A0:B0)=Valwnew(A0:B0)+ ...
Valwprev(mm).*pt(mm,A0:B0);
end
%ECovZnew0(tt1,tt2,1:Nn2)=ValZnew(1:Nn2);
%ECovwnew0(tt1,tt2,1:Nn2)=Valwnew(1:Nn2);
ECovZprev0(tt1,1)=sum(ValZnew(1:41));
ECovZprev0(tt1,2:Nn-1)=ValZnew(42:40+Nn-1);
ECovZprev0(tt1,Nn)=sum(ValZnew(40+Nn:80+Nn));

ECovwprev0(tt1,1)=sum(Valwnew(1:41));
ECovwprev0(tt1,2:Nn-1)=Valwnew(42:40+Nn-1);
ECovwprev0(tt1,Nn)=sum(Valwnew(40+Nn:80+Nn));

%Below when second covariance pair date is equal to present date, we do
%final(ending) calculations for covariance calculation with all first pairs from
%all of the earlier dates.
ValZnow(1:Nn)=Zwt(tt,1:Nn)-MuZ(tt);
Valwnow(1:Nn)=wt(tt,1:Nn)-Muw(tt);
ValZprev(1:Nn)=ECovZprev0(tt1,1:Nn);
Valwprev(1:Nn)=ECovwprev0(tt1,1:Nn);
CovarZt(tt1,tt)=sum(ValZprev(1:Nn).*ValZnow(1:Nn));
Covarwt(tt1,tt)=sum(Valwprev(1:Nn).*Valwnow(1:Nn));

end
end
% %#Covariances Calculation.
% %Below, we calculate covariances between the each date with all of the
% %later following dates
%
%     if(tt==1)
%         VarZ(tt)=sum((Zwt(tt,1:Nn)-MuZ(tt)).^2.*ZProb(1:Nn));
%         Varw(tt)=sum((wt(tt,1:Nn)-Muw(tt)).^2.*ZProb(1:Nn));
%     end
%     if(tt>1)
%         VarZ(tt)=sum((Zwt(tt,1:Nn)-MuZ(tt)).^2.*ZProb(1:Nn));
%         Varw(tt)=sum((wt(tt,1:Nn)-Muw(tt)).^2.*ZProb(1:Nn));
%
% %Below loop with variable tt1 indicates the starting time for various
% %covariances
% %Below loop variable tt2 indicates the algorithm calculations on dates
% %following the strting time.
%         for tt1=1:tt-1
%             for tt2=tt:Tt
%                 tt1
%                 tt2
%                 if(tt1<tt-1)
% %If the starting time for covariance calculations is smaller than one date
% %ago, it means that we do not have to initialize the algorithm and we are
% %in the middle of the calculations.
%                     ValZprev(1:Nn)=ECovZprev0(tt1,tt2,1:Nn);
%                     Valwprev(1:Nn)=ECovwprev0(tt1,tt2,1:Nn);
%                 end
%                 if(tt1==tt-1)
% %If the starting time for covariance calculations is equal to one date
% %ago, it means that we haveto initialize the algorithm and we are
% %starting the calculations. Here we assign the starting values for covariance
% %calculation.
%                     Zwt0(1:Nn)=Zwt(tt-1,1:Nn);
%                     ValZprev(1:Nn)=(Zwt0(1:Nn)-MuZ(tt-1)).*ZProb(1:Nn);
%                     wt0(1:Nn)=wt(tt-1,1:Nn);
%                     Valwprev(1:Nn)=(wt0(1:Nn)-Muw(tt-1)).*ZProb(1:Nn);
%                 end
%                 ValZnew(1:Nn2)=0;
%                 Valwnew(1:Nn2)=0;
% %Below we calculate the running string for the transition probability increment
% %between first covariance date and last covariance date.
%                 for mm=1:Nn
%                     A0=ptA(mm,tt);
%                     B0=ptB(mm,tt);
%                     pt(mm,1:Nn2)=ptt(tt,mm,1:Nn2);
%                     ValZnew(A0:B0)=ValZnew(A0:B0)+ ...
%                         ValZprev(mm).*pt(mm,A0:B0);
%                     Valwnew(A0:B0)=Valwnew(A0:B0)+ ...
%                         Valwprev(mm).*pt(mm,A0:B0);
%                 end
%                 %ECovZnew0(tt1,tt2,1:Nn2)=ValZnew(1:Nn2);
%                 %ECovwnew0(tt1,tt2,1:Nn2)=Valwnew(1:Nn2);
%                 ECovZprev0(tt1,tt2,1)=sum(ValZnew(1:41));
%                 ECovZprev0(tt1,tt2,2:Nn-1)=ValZnew(42:40+Nn-1);
%                 ECovZprev0(tt1,tt2,Nn)=sum(ValZnew(40+Nn:80+Nn));
%
%                 ECovwprev0(tt1,tt2,1)=sum(Valwnew(1:41));
%                 ECovwprev0(tt1,tt2,2:Nn-1)=Valwnew(42:40+Nn-1);
%                 ECovwprev0(tt1,tt2,Nn)=sum(Valwnew(40+Nn:80+Nn));
%            end
% %Below when second covariance pair date is equal to present date, we do
% %final(ending) calculations for covariance calculation with all first pairs from
% %all of the earlier dates.
%            ValZnow(1:Nn)=Zwt(tt,1:Nn)-MuZ(tt);
%            Valwnow(1:Nn)=wt(tt,1:Nn)-Muw(tt);
%            ValZprev(1:Nn)=ECovZprev0(tt1,tt,1:Nn);
%            Valwprev(1:Nn)=ECovwprev0(tt1,tt,1:Nn);
%            CovarZt(tt1,tt)=sum(ValZprev(1:Nn).*ValZnow(1:Nn));
%            Covarwt(tt1,tt)=sum(Valwprev(1:Nn).*Valwnow(1:Nn));
%
%         end
%     end
end

%At the end, we want to see the results for our covariance calulations
for tt=1:Tt-1
%Below, we see the variance of standard increments normal on a particular date
VarZ(tt)
%Now, we want to see the covariance of all the following dates with present date
%and for normal noise this should be equal to the variance at starting
%date.
CovarZt(tt,tt+1:Tt)
tt
str=input('Look at CovarZt');
%Below, we see the variance of a general SDE(but this current program works well
%with driftless noises only) on a particular date
Varw(tt)
%Now, we want to see the covariance of all the following dates with present date
%and for CEV type driftless SDEs, this should be equal to the variance at starting
%date.
Covarwt(tt,tt+1:Tt)
tt
str=input('Look at Covarwt');
end

plot((1:Nn),ZProb(1:Nn),'g',(1:Nn),Pprev(1:Nn),'r');
%plot((1:Nn),ZProb(1:Nn),'g',(1:Nn),Pprev(1:Nn),'r');
%str=input('Look at distributions');

Pnew=Pprev;
%Convert the expected values in each cell to standard Values associated with the grid cell after removing the
%integrated probability in the cell.
w(wnStart:Nn)=Ewprev(wnStart:Nn)./Pprev(wnStart:Nn);
w1(wnStart:Nn)=Ew1prev(wnStart:Nn)./Pprev(wnStart:Nn);
yy(wnStart:Nn)=((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
yy1(wnStart:Nn)=((1-gamma)*w1(wnStart:Nn)).^(1/(1-gamma));

%below D's (the names of variables starting with D) are
%change of probability derivatives.

y_w(1:Nn)=0;
y_w(wnStart:Nn) = ((1-gamma)*w(wnStart:Nn)).^(1/(1-gamma));
y_w1(1:Nn)=0;
y_w1(wnStart:Nn) = ((1-gamma)*w1(wnStart:Nn)).^(1/(1-gamma));
Dfy_w(wnStart:Nn)=0;
Dfy_w1(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));
Dfy_w1(nn) = (y_w1(nn + 1) - y_w1(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
py_w(1:Nn)=0;
py_w1(1:Nn)=0;
for nn = wnStart:Nn-1
py_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density
py_w1(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w1(nn));
end

toc

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

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

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

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

end

YY(YY<0)=0;
disp('Original process average from monte carlo');
MCMean=sum(YY(:))/paths %origianl coordinates monte carlo average.
MCVar=sum((YY(:)-MCMean).^2)/paths
disp('Original process average from our simulation');
ItoHermiteMean=sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates
ItoHermiteVar=sum((y_w(wnStart+1:Nn-1)-ItoHermiteMean).^2.*ZProb(wnStart+1:Nn-1))
disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average

MaxCutOff=30;
NoOfBins=round(300*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too
[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti_NEW(YY,paths,NoOfBins,MaxCutOff );

%plot(y_w1(wnStart+1:Nn-1),py_w1(wnStart+1:Nn-1),'b',y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
plot(y_w(wnStart+1:Nn-1),py_w(wnStart+1:Nn-1),'r',y_w1(wnStart+1:Nn-1),py_w1(wnStart+1:Nn-1),'b',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');

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

legend({'Improvised terms Model Density','Monte Carlo Density'},'Location','northeast')

str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
end
.
Here is the last output of the program.
MCMean =
0.999338539579104
MCVar =
0.195295297580708
Original process average from our simulation
ItoHermiteMean =
1.000013291041896
ItoHermiteVar =
0.195216609800823
true Mean only applicble to standard SV mean reverting type models otherwise disregard
TrueMean =
1
IndexMax =
922

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

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

I am writing this post specifically to let friends know that mind control is not decreasing at all. Today since I was being able to do better work, they specially used my family to give me large amount of mind control drugs and I lost all momentum and mental stimulation. And when I tried to post on internet, they controlled my internet on my computer and I could not post despite being able to login and repeatedly trying for fifteen minutes and then I am making this post using my ceĺl phone something I almost never use for internet. Mind control agency also actively continurs to drug my fòod in the markets and several suvh incidents happened in past few days. I want to mention here that date for giving me next antipsychotic injection is due in another ten days and my family on orders of mind control agency will be very strict and might even use force to continue to give me injections. I always hopef that I have had my last injection and nobody will force injections on me now but things srem to go other way and mind control agencies are still having a great home run. Please force the mind control agency to end my mind control and mind control of other countless victims like me. It is such s pity that mind control continues and thrives even after arrival of Biden administration. I want to request all good Americans to please not let crash hopes of a large number of victims that we had with the new Biden government.

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

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

I have not written about my activities in past one week other than the two posts I made yesterday. My routine has been very disarranged for the time. My mother had left the hospital and she is recovering very slowly due to old age. And somebody has to attend to her to move her around and I am stronger than my sister and my father who are here to attend to her and therefore I have to help her with her movements. And my mother usually has to go to washroom at 2:30 am or 3:30 am so this becomes sort of difficult and therefore my sleep and eating patterns were a bit disturbed (though I was still doing reasonably good research in some patches) and therefore I was not thinking much about writing here at Wilmott until things remarkably deteriorated. My father has already left Islamabad due to his other engagements and I would have to spend more time looking after my mother in coming days.
I want to tell friends about another interesting thing. For at least past two years almost since I bought this (laptop) computer, I have had a suspicious rogue network connection on my computer. I tried to disable and remove it(from settings) tens of times but it would reappear and then I just gave up on removing it. So my computer usually shows two network connections telling me that my wi-fi has internet access and explicitly mentions that other network connection has no internet access. But to my great chagrin, I have noticed several times when I am surfing internet that for brief periods of time this unidentified network connection would disable my wi-fi connection and all of a sudden this rogue network connection would have internet access and all my internet traffic would seem to be routed through this unidentified network connection that almost always had no internet as a rule. May be the hardliner and machinatory crooks of pentagon use these tricks to show people that some internet traffic is coming out of my computer while I have absolutely no idea about it. I would request friends to be very careful if they know some bad material I am surfing because some of it might have been intentionally planted out of my computer by machinatory crooks. These crooks have more control over my computer than I have myself. In fact there was a time that I used to keep my computers on me all the time(They still broke into my computer from internet or other wi-fi/bluetooth devices) but my family who would allow these crooks unrestricted access to my computers would call me schizophrenic so I just gave up on caring too much and would freely leave my computers at home while I would be away and these crooks got much bolder.
I have been exposing the hardliner and racist crooks of Pentagon for so many years now. The more evil these crooks are, the more noble aura they create around themselves with the help of their lies, trickery and cunning. These righteous crooks are guided by their hardliner right wing ethics and they are well aware that world does not follow their righteous right wing ethics and therefore they are steeped in the tactics based on lies, trickery , cunning  and machinations so that they could justify to rest of the world what is otherwise truly only their righteous hardliner right wing racist agendas.
I would really request good Americans to stress on the new president to bridle these righteous crooks and let so many good-natured, intelligent people become free who became victims to these crooks only due to their intelligence and righteous ethics of the crooks in pentagon.

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

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

Very sorry friends. I thought I had made the last post in Off-topic thread but not knowingly posted it in technical forum. Though I continue to post some non-technical material in technical forum sometimes, I would really prefer it to stay technical if it were not for my other woes. Sorry again and ,in the future, I would mostly just give links to Off-topic thread here in technical forum if it is not something extremely important. Please pardon this mistake.

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

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

viewtopic.php?f=15&t=94796&p=864038#p864038
Friends, I have outlined ideas about how to plant a forest in Sahara desert to cool down the earth's atmosphere. It seems lunacy but it can actually be a wildly profitable enterprise on its own even not including any benefit to earth's Climate and atmosphere. I have given simple back of the envelope calculations about costs and profits associated with planting a jungle in Sahara desert.
viewtopic.php?f=15&t=94796&p=864038#p864038
viewtopic.php?f=15&t=100778

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

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

Friends, I believe I have a good idea about how to complete the research on transition probabilities project. Basically we have a Z-grid that expands at the rate of unit normal. CEV type noises work well on the grid since CEV noise volatility remains proportional to underlying Z-grid volatility. For noises, given a certain noise/innovation volatility, the grid expands faster when existing Variance associated(proportional to Z-grid) with Z-grid is smaller -- This is why a noise grid expands very fast at the start of simulation when existing variance is low and barely budges when the time is large and the grid has expanded to a large variance i.e the effect of a certain noise volatility is inversely proportional to existing variance. However when there is drift in an SDE that significantly affects variance, the underlying Z-grid and simulated SDE are no longer proportional to each other. For example when a mean-reverting SDE is simulated, the drift decreases effective variance while underlying Z-grid continues to expand at the same unit variance rate. As we observed earlier, on an over-expanded grid, the effect added by a certain noise volatility is smaller than it would be when contraction in grid due to mean-reverting drift would also be accounted so we get an over all decreased variance in case of mean-reverting SDEs (when we do not explicitly account for effect of variance contraction on underlying Z-grid). Since it would be relatively difficult to alter the base Z-grid(though we may possibly be able to try special grids for mean-reversion SDEs), I am thinking of ways how to just add the effect of contraction in existing variance to the noise volatility without altering the base Z-grid(made of unit normals).
Again, regarding the underlying Z-grid, it might be relatively simpler to make a Z-grid with adaptive variance that includes the effect of contraction by the drift on the fly at the same time as we do the main simulation but that would be far more expensive since all linear equations would have to be solved at the same time when doing the main simulation.
Friends, my mother has almost recovered from diabetic/sugar complications but she is still in the hospital and would be returning home in a day or two. I have been able to work today after a long time since I did not go to care for her to the hospital. I would be caring for her for a lot of time when she comes home but that should not be obstacle to my research since I would continue to work on my computer while attending her from time to time.
I also want to tell friends that an antipsychotic injection is due in a few days and I have plans to resist getting this injection. I want to request friends to force mind control agencies to not force my family to give me any more injections. Please ask mind control agencies to not give me antipsychotics any more and let me live normally.

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

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

After I shared the last post in which I mentioned that variance of Z-grid has to be altered to take into account change in variance due to drift, I had some other thoughts that I also want to share with friends. I have several time tried simple simulations(without random numbers) on Z-lines ( similar to what we did in solution of Fokker-planck equation) by solving stochastic integrals in forward time(something I shared a month or two ago but had some errors in integrals of second hermite polynomial but for SDEs in  bessel form we do not need those integrals of second hermite polynomials). We basically continue to integrate those stochastic integrals of standard brownain motion along Z-lines to solve the SDEs. It works remarkably well for simple CEV type noises in bessel form but I never shared it since it did not work for mean-reverting SDEs. Now I realize that this method can also work for mean-reverting SDEs if we could account for decrease in variance of the underlying brownain motion integrals due to mean-reversion. Only difference with fokker-planck solution is that Fokker-planck solution automatically calculates any change in variance due to drift while other monte carlo-like(with random numbers replaced with stochastic integrals) would have to explicitly account for any contraction due to drift and would have to make extra calculations for that. After completing the part on transition probabilities, I hope to present the programs that solve for densities of SDEs on Z-lines(a tree with just one branching at the origin) with simple stochastic integrals.

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

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

I want to quickly give friends a toy example about how a decrease in variance of SDE due to drift effects its future evolution. Suppose we have a very simple SDE of the form  $dw(t)=\mu(w(t)) dt+ \sigma dz(t)$
if we neglect higher order terms and only look at most important terms, a naive evolution of SDE along a particular Z-line with principal order terms would be
$W(t_1+dt)=W(t_1)+ mu(w(t_1)) *dt + \sigma * (Z \sqrt{t_1+dt}-Z \sqrt{t_1})$
But suppose that due to effect of earlier non-linear drift, the effective equivalent Brownian variance has already decreased before this simulation interval slightly from $t_1$ to a smaller value $s_1$. In that case our simulation of principal terms should change to
$W(t_1+dt)=W(t_1)+ mu(w(t_1)) *dt + \sigma * (Z \sqrt{s_1+dt}-Z \sqrt{s_1})$

since $s_1$ is smaller than $t_1$,   the value of volatility part $\sigma * (Z \sqrt{s_1+dt}-Z \sqrt{s_1})$  is now larger than if we had not modified the variance of underlying Z(t) and kept it at unaltered value with the resulting smaller volatility increment(if we had not switched to decreased effective variance) equal to  $\sigma * (Z \sqrt{t_1+dt}-Z \sqrt{t_1})$
In a proper simulation of bessel form SDEs, we would have to include stochastic integrals of order 1.5 and possibly other values but the essence remains the same.
Though there might be some other minor issues, if we can keep track of effective equivalent brownian variance, we can easily evolve most difficult SDEs along Z-lines like we did in fokker-planck with a simple monte-carlo like discretization with analytic stochastic integrals. We would just have to keep track of decrease or increase of local variance due to drift and project it on underlying brownian motion or underlying standard increments based normal. The only major different with random numbers based monte carlo is that full-fledged monte carlo increments do not have to explicitly take into account decrease of underlying Z-variance/brownian variance while in our new analytic form(without random numbers) we have to explicitly account for this change in variance due to drift.
This is why unless we alter the effective underlying Z-variance on every step of simulation, our resulting simulation of mean-reverting SDEs has a smaller variance(due to smaller than true volatility effect at each step)