Serving the Quantitative Finance Community

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

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

Friends, I have been doing work on calculation of forward expectations of functionals of SDE on a transition probabilities grid. I have been trying to put everything on a reasonable mathematical foundation. In forward simulation important problem is how to faithfully retrieve the expectations of events that occurred on the previous time grids. Usually, though not always, variance of underlying SDE distribution increases as we advance into next(or later) time grids. But expectations of SDE related functions evaluated on previous time grids were calculated with respect to SDE distribution at previous times when underlying SDE variance was (usually) slightly lesser. So in order to faithfully retrieve the densities of functionals of SDE in future time grids we have to add some appropriate variance. This added variance makes sure that we can retrieve the functionals densities with expectations evaluated at earlier times with smaller underlying variance faithfully with respect to each of the later time SDE grids when underlying SDE distribution variance has increased. But a side effect of this addition of variance is that covariances across time also increase in proportion to added variance. Expectations of Path integrals are highly sensitive to covariances of the functionals across different time grids so we have to make an accounting for how much appropriate relevant covariances have increased. This can be done in several ways since we know how much extra variance was added to expectations of functionals on each time step to convert the functionals expectation originally calculated with respect to lower variance underlying SDE grid into expectations with respect to SDE grids at later times with higher variance.
I am working on various mathematical issues so that we can find expectations of functionals of every SDE even when we do not analytically know the SDE covariances across time.
It may be still a few days before I could post a perfectly working program since I want to change a lot of things in my previous program as my understanding of the problem has become better. I will be posting comments regularly for next few days before I post a full-fledged worked out program.

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

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

There is some different progress so I thought I will mention it here for friends. Though I had a completely worked out plan how to convert the expectations of functions from an earlier SDE grid to a later time grid, I just thought that I should try measure theory for this purpose and if things could work, it would be a far more elegant solution to finding the densities of functionals of SDE. Towards that goal, I wrote down exponential likelihood function of general (possibly mean-reverting) SDEs with drift and state dependent variances (By doing Ito change of variable on log of SDE variable) and I was able to faithfully transport some very simple payoffs forward in time. However, it still remains to be seen if slightly more difficult payoffs could be transported in time and I still have to work to see that. As only seeing is believing, I still cannot be sure if the measure transport method would universally work but I am quite optimistic. I will get more detailed results in another day and will let friends know how it goes as my research progresses.

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

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

Friends, I would not be writing this post if it were not for torture that started four days ago. I was working on my projects at full speed when mind control agents targeted me with drugged water and started giving me attacks of very strong anxiety and I was totally unable to work for past four days. These attacks of very strong anxiety lasted everyday after that and therefore I decided to write more about mind control and not about mathematical finance that all of us including myself love so much. I have not so huge interest about whining about mind control but if the torture continues, I have lined up many posts about mind control that will follow.
I am copying an old post of mine and want to ask friends if mind control agencies were serving United States or trying to please their godfather. I still recall when mind control agent who made more than 800 million dollars and when he would narrate me incidents in which he would describe to his jewish backer professors at my previous school the stories of how targeted talented Pakistani muslims would beg to end mind control torture on them and then mind control agent would tell me how excited, amused and delighted some jewish professors at my old school would be when they would hear how talented muslim targets were desperately begging to end mind control that would never end. I am not lying and therefore the agent was able to make more than 800 million dollars.
Again I do not want to write such stuff but if my mind control continues, there is a lot more on the way.
I am copying an old post of mine and want to ask friends if mind control agencies were serving United States or trying to please their godfather.
.
.
Thank you Paul for your post. I will look at the website you mentioned. But I want to really thank you for letting me write bold things on your forum that most forum owners would not allow. Though I might have different ideas than you on a lot of issues, I am a staunch supporter of freedom of speech just like you are and I am sure both of us agree on one slogan,"Je suis Charlie". I consider myself lucky to have been writing on your forum since 2003.

I want to write this post to ask American tax payers to please cut funding of crooks in US army. I was reading that there has not been a successful audit of spending by US army in past decade. You would probably ask me why you should do anything to bridle spending by crooks in army. Here is why.
Crooks in US defense have spent several billion dollars in attempts to retard me over past twenty years. I went to Hong Kong in 2000 to find a job in finance and derivatives and I stayed there for three months. Crooks in US army continues to drug food and water at large number of places in Hong Kong and Kowloon and continued to play all sort of tricks with me. I still recall when I went to office of South China Morning Post to tell them how badly they had drugged the city of Hong Kong, they had already set up someone who received me and told me after hearing my story that he did not believe me.
They gave me electric shocks twenty days after I returned from Hong Kong.
Later I was asked by my Japanese employer to come to Tokyo (2005) and work in their office there(I was earlier working from Lahore). They drugged large number of places in Tokyo city and mad e my life absolute hell and I resigned only after one month of stay in Tokyo. I did not want to return to Pakistan due to huge persecution by my family and the army so I decided to go somewhere where US military could not go and got Iranian visit visa from Tokyo and went to Tehran.
Crooks in US defense gave huge amount of money to Iranian crooks(crooks everywhere know each other well) and they used huge microwave and gases torture on me. I still recall on my first night in an Iranian hotel, I left the hotel after midnight and randomly ran in streets of Tehran to somehow avoid pain due to directed microwave torture. Just after two or three days in Tehran, I was trying to find refuge in embassies like Vatican in Tehran asking them to help me on human grounds. I left Tehran after one week and flew back to Pakistan. I still recall that I arrived at Tehran airport a lot earlier than flight departure. I had previously protested somewhere against unnecessary removal of shoes of Muslim passengers and it was known to mind control agencies. They asked Iranians security on airport who approached me inside the departure lounge and took me to a room and asked me to remove my shoes and had my body search.
Since during all this, mind control agencies could not control me well and several of my neurotransmitters were coming back, they had extremely huge amount of gases in washrooms in Dubai airport where I would stop for transit while going to Karachi. No wonder ultra conservative crooks in US army align so well with Saudi(who dissolve their criticizers in acid) and Dubai.
Later I got British HSMP work visa and went to London in 2010. Crooks in US defense drugged several large cities of UK on a large scale(it was impossible to get good beverages or bottled water. They even drugged beer, wine and pork at a lot of places.) You can ask people in UK secret services(many of whom would love to tell you everything in detail) about supernatural stories how they drugged food in cities of UK in 2010. I still recall that I wrote several page account of my persecution and I would print hundred of copies of it and distribute them in London. When I was distributing those papers outside of some London Newspapers, administration of those newspapers went to extreme length to point out that certain part of sidewalk was their property and they would call police if I stepped there. I was forcefully removed from the compound in London(after I had passed all security checks) when I went there to ask law firms for help and the law firms on main street would refuse to meet me and say that I had to either email or send the letter by post and they would not accept any papers from me(somebody would approach me at the door and tell me about it and ask me to leave.).
They continued to drug food and beverages in entire city of Lahore numerous times in past ten years.
Net worth of several mind control agents who continued to retard me in Pakistan is north of several hundred million dollars each depending upon their seniority(with most senior agent worth 800 million USD).
And this drugging of entire cities by crooks of US defense never ends. Now they are drugging the twin cities of Rawalpindi/Islamabad.
I am just telling about myself. There are tens of thousands of mind control victims all over the world.
Please cut their funding before it is too late. They are not retarding just Muslims(I have no religion but of course I am Muslim by birth) only, they retard a large number of totally innocent and talented Americans. Crooks love to bite the same hand that feeds them(It is typical of crooks and large unbridled armies in human history)

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

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

When I look at similarities and parallels, I find it interesting to notice that United States is becoming, just like reactionary act of mind control, a country with more and more reactionary policies. For example I find that there is more and more concern in United States about increasing might of China and a growing need to contain China. I want to tell Americans that this is very misguided. If Americans were truly wise, they would never try to "contain" anyone, and instead of containing anyone they would just do what needs to be done to make them a "country of excellence" something they did so well over the past one century. When USSR was getting ahead in sciences and technology, America responded so perfectly by emphasizing science and technology in its universities and Americans to this day reap the fruit of great decisions for advancing science and mathematics in its academic institutions. Now when China is advancing so fast, Americans need to make sure that they truly invest in education of their children and make sure that every intelligent child becomes a star in its own right. Instead of so many other trillion dollar packages, Joe Biden would be far wiser to spend 500 billion dollars to make sure that every child in US whether he is born in inner part of some large city or some remote state would get the same great quality education, he would no longer need to worry about rising Chinese might as American fortunes would continue to surpass any Chinese miracles. I really want to say to  American friends, "Do not be reactionary. Just be original and genuine. Just do even better what you did so well over past one century."
It is interesting that the very people who could not get great education and had to become crooks in mind control agencies instead of being doctors, lawyers and successful businessmen and then these crooks do all sort of evil things in the society thinking that only way to get ahead in life is to join Godfather's paradise. Only great education in backward states would decrease the mass of crap crooks doing evil things in US society.

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

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

Friends, I have been writing on this thread since 2016 with occasional large breaks. I am going to take another break of 3-4 weeks now. I had been working on a market trading algorithm for quite some time and I took to improving it after my post of 25th September. One reason my persecution continues is that I have not been able to make enough money to live independently from my family. So I was thinking if I could write a great market trading algorithm and sell it, I could try to get out of this problem of continuing mind control. I have been putting all my research other than market trading research on internet so its value for me in terms of money sharply decreases so it was somewhat difficult to leverage it to make enough money. Good thing is that my market trading research has matured enough and I have been able to write great trading programs. Currently I am trying some methods to further improve my results. In two months' data I was able to find returns ranging from 40% of notional to larger than 100% of notional which means that returns per month would range from 20% to 50% of notional and I am not using any leverage. My algorithm is a market-making algorithm and enters the trade only using limit orders but it is not a high frequency algorithm at all. Average time between  trades ranges from one minute to 5 minutes and algorithm waits for opportune moments to start and end a trade. But you would like to be filled quickly once the algorithm wants you to do that. I still have to calculate sensitivity to order execution slippage but it should not be very large. Algorithm would work on highly liquid stocks but could be modified to deal with moderately liquid stocks. One reason for high returns could be that I do not have to scale the algorithm to very large amount of capital where investing like this would become difficult and would surely reduce the returns on large capital though I still think that if larger amounts of capital are judiciously invested, returns would still be far-better than many mainstream approaches.
Obviously, I am not going to post the algorithm since if I do that its effectiveness might end in a few weeks.
So I will be busy improving the algorithm and marketing it. I want to sell the algorithm to just 2-3 clients for money in six figures. I am very afraid that mind control agencies would do everything in their power to thwart me here.
If I could successfully sell the algorithm, I will hire a team of 5-6 mathematicians and computer programmers and train them in math finance and mathematics of trading. I want to tell friends that I actively want to continue doing research in mathematical finance and would continue posting very interesting things and research in the future with a more professional setup. I want to  keep posting my research and programs on Wilmott for several years into the future.
My post would be incomplete if I do not thank so many American people who protested to mind control agencies against my continuing persecution when it was becoming very difficult for me to continue my research  and mind control activity decreased due to protests from those good people and I was able to restart my research. And I also want to thank so many readers of this forum across the world who wished me well and wanted me to succeed and do good research. And I want to say to Americans that even if I succeeded in life, I will try to be a good human being and they will find a friend of American people in me.
One last thing, please ask mind control agencies who have access to my computer to not distribute my program in their favorites. Wrongdoing must not be encouraged like that.

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

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

I want to ask Americans a simple question. Isn't this true that if I were a jew, people all over the world would be full of praise for me. On the contrary since I am born to a muslim family, everyday I am subjected to inhuman mind control torture and several people in American government and defense continue to machinate how to fail me. Is this great American democracy all of us are never tired of praising?

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

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

Friends, please pardon my previous post. I was extremely down due to biweekly antipsychotic injection. I get half a Fluanxol 40mg injection after every fifteen days and for first seven days after the injection, I remain very down and usually unable to do anything meaningful. I become better seven days after injection and become relatively creative. I am 48 years old and it is very difficult for me to stand these regular antipsychotic injections but mind control agencies connived with my father and family and have forced injections and drugs on me for more than past twenty years. I still remember when I was 28/29 years of age, at one stage, I was given many of these injections per week but I still had energy and vitality but now I am too old to continue to take these injections regularly that are faced on me and it is impossible to  remain active and many times I fall into despair.
Anyway good thing is that seven days have passed since last injection and I have some ideas about increasing the order (in time) of solution of partial differential equations. I also have some ideas about converting my old derivation of derivatives of Fokker-planck equation into a formal solution (unlike informal working solution I had presented earlier.) I still cannot say what will come out of these ideas but I will be working on them for next several days and if there is any progress then I will post my new programs. I hope to share some interesting stuff with friends in next few days.

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

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

Sorry friends, I prematurely hit the submit button and so could not complete the previous post. Here is the actual relevant post. Please disregard the previous incomplete post.
Friends, I wanted to share a rough sketch of the proof of the method I presented in above programs. I hope friends would like it.
I am not writing Fokker-planck equation in Bessel coordinates and I am sure most friends are familiar with it.
In our setting we have an SDE variable, $B$, in Bessel coordinates on nth grid point $Z_n$ of underlying Z-grid. This is represented as $B(Z_n)$. As a function of time, $B(Z_n)$ is moving so as that CDF of the density at point $B(Z_n)$ remains constant. This is written in equation form as below

$\frac{\partial}{\partial t} \big[\int_{-\infty}^{B(Z_n)} P(B) dB \big]=0$

applying the time derivative to terms in brackets, we get

$\int_{-\infty}^{B(Z_n)} \frac{\partial P(B)}{\partial t} dB + \frac{\partial B}{\partial t} P(B(Z_n)) =0$

I will make it more formal but here, since we add drift separately and because we are in Bessel coordinates, I am just replacing $\frac{\partial P(B)}{\partial t}$ by $.5 {\sigma}^2 \frac{\partial^2 P(B)}{\partial B^2}$ in the above equation.
Making the substitution in above equation, we get

$\int_{-\infty}^{B(Z_n)} \big[.5 {\sigma}^2 \frac{\partial^2 P(B)}{\partial B^2} \big] dB + \frac{\partial B}{\partial t} P(B(Z_n)) =0$

Applying the integral on first term, we get

$\big[.5 {\sigma}^2 \frac{\partial P(B(Z_n))}{\partial B} \big] + \frac{\partial B}{\partial t} P(B(Z_n)) =0$

Now we want to make a change from $P(B)$ to $P(Z)$. In order to do that we have two equations
$P(B)=P(Z)\frac{\partial Z}{\partial B}$ and
$\frac{\partial P(B)}{\partial B}=\frac{\partial P(Z)}{\partial Z} {(\frac{\partial Z}{\partial B})}^2 + P(Z) \frac{\partial^2 Z}{\partial B^2}$
Now substituting above two equations in previous equation, we get
$.5 {\sigma}^2 \big[\frac{\partial P(Z_n)}{\partial Z} {(\frac{\partial Z}{\partial B})}^2 + P(Z_n) \frac{\partial^2 Z}{\partial B^2} \big] + \frac{\partial B}{\partial t} P(Z_n) \frac{\partial Z}{\partial B}=0$
substituting $\frac{\partial^2 Z}{\partial B^2} =-{(\frac{\partial Z}{\partial B})}^3 \frac{\partial^2 B}{\partial Z^2}$
we get
$.5 {\sigma}^2 \big[\frac{\partial P(Z_n)}{\partial Z} {(\frac{\partial Z}{\partial B})}^2 - P(Z_n) {(\frac{\partial Z}{\partial B})}^3 \frac{\partial^2 B}{\partial Z^2} \big] + \frac{\partial B}{\partial t} P(Z_n) \frac{\partial Z}{\partial B}=0$
cancelling $\frac{\partial Z}{\partial B}$ on both sides, we get
$.5 {\sigma}^2 \big[\frac{\partial P(Z_n)}{\partial Z} {(\frac{\partial Z}{\partial B})} - P(Z_n) {(\frac{\partial Z}{\partial B})}^2 \frac{\partial^2 B}{\partial Z^2} \big] +\frac{\partial B}{\partial t} P(Z_n)=0$
The first term in above is expansion that we have treated in the previous matlab program in a different analytic way by adding squared volatilities . Matching the coefficients of second and third term, we get

$\frac{\partial B}{\partial t} = .5 {\sigma}^2 \big[{(\frac{\partial Z}{\partial B})}^2 \frac{\partial^2 B}{\partial Z^2} \big]$

the last equation is the analytic solution that I have used in my previous matlab programs to get the density of SDEs in Bessel coordinates right.
.
.
.
Here I want to explain how the above analytics could be applied to Fokker-planck equation in original coordinates (and it works seamlessly).
First I write the fokker-planck equation in original coordinates as
$\frac{\partial P(x,t)}{\partial t} = -\frac{\partial \big[\mu(x)P(x,t) \big]}{\partial x}\, +.5 \, \frac{\partial^2 \big[{\sigma(x)}^2P(x,t) \big]}{\partial x^2}$

We have our partial differential equation on a grid and we want to determine the movement of an arbitrary grid point (boundary) along time so that mass within(total mass up till the grid point is conserved or associated CDF remains constant) that grid point (boundary) remains constant as a function of time. we denote this arbitrary grid point as $x_b$.
We write the conservation of probability mass up till (or associated constant CDF) the grid point (boundary) as a function of time in equation form as

$\frac{\partial}{\partial t} \big[\int_{-\infty}^{x_b} P(x,t) dx \big]=0$

applying the time derivative to terms in brackets, we get

$\int_{-\infty}^{x_b} \frac{\partial P(x,t)}{\partial t} dx + \frac{\partial x_b}{\partial t} P(x_b,t) =0$

But we know that term inside the integral is given by FP equation and is LHS of FP equation. We replace the time derivative inside the integral  by RHS of FP equation as

$\int_{-\infty}^{x_b} \big[ -\frac{\partial [\mu(x)P(x,t)]}{\partial x}\, + \, .5\frac{\partial^2 \big[{\sigma(x)}^2P(x,t) \big]}{\partial x^2}\big] dx + \frac{\partial x_b}{\partial t} P(x_b,t) =0$

Applying the integral to complete differentials inside square brackets, we get

$-\mu(x_b)P(x_b,t)\, + \, .5 \, \frac{\partial \big[{\sigma(x_b)}^2P(x_b,t) \big]}{\partial x} + \frac{\partial x_b}{\partial t} P(x_b,t) =0$
$=-\mu(x_b)P(x_b,t)\, + .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}P(x_b,t) + .5\,{\sigma(x_b)}^2 \, \frac{\partial P(x_b,t)}{\partial x} +\, \frac{\partial x_b}{\partial t} P(x_b,t) =0$

We have solved the equation in original coordinates and we could do an ODE in $\frac{\partial x_b}{\partial t}$ to find the evolution of every grid point so that probability mass till that grid point is conserved. But it will be a lot easier if we converted the above arbitrary density to underlying gaussian density as we have been doing and for that we write the following substitution equations from change of probability derivatives as($Z_b$ is associated with each point $x_b$ using common CDF as probability is being conserved.)

$P(x)=P(Z)\frac{\partial Z}{\partial x}$ and

$\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}$

Substituting the above two equations in previous equation, we get

$-\mu(x_b) \, P(Z_b)\, \frac{\partial Z}{\partial x} \,+ .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}P(Z_b) \frac{\partial Z}{\partial x} + .5\,{\sigma(x_b)}^2 \, \big[ \frac{\partial P(Z_b)}{\partial Z} {(\frac{\partial Z}{\partial x})}^2 + P(Z_b) \frac{\partial^2 Z}{\partial x^2} \big] +\, \frac{\partial x_b}{\partial t} P(Z_b) \frac{\partial Z}{\partial x}=0$

We make the following substitutions in above equation

$\frac{\partial^2 Z}{\partial x^2} =-{(\frac{\partial Z}{\partial x})}^3 \frac{\partial^2 x}{\partial Z^2}$

and

$\frac{\partial P(Z_b)}{\partial Z}=-Z_b \, P(Z_b)$

$-\mu(x_b) \, P(Z_b)\, \frac{\partial Z}{\partial x} \,+ .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}P(Z_b) \frac{\partial Z}{\partial x} + .5\,{\sigma(x_b)}^2 \, (-Z_b) \, P(Z_b) {(\frac{\partial Z}{\partial x})}^2 - .5\,{\sigma(x_b)}^2 P(Z_b) {(\frac{\partial Z}{\partial x})}^3 \frac{\partial^2 x}{\partial Z^2} +\, \frac{\partial x_b}{\partial t} P(Z_b) \frac{\partial Z}{\partial x}=0$

Cancelling  $P(Z_b) \frac{\partial Z}{\partial x}$ throughout the equation and rearranging, we get the first order very simple ODE for our grid point that conserves the probability mass up till that grid point.

$\frac{\partial x_b}{\partial t} =\mu(x_b) - .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}+ .5\,{\sigma(x_b)}^2 \, Z_b \, \frac{\partial Z}{\partial x} + \, .5\,{\sigma(x_b)}^2 \, {(\frac{\partial Z}{\partial x})}^2 \frac{\partial^2 x}{\partial Z^2}$

The above is the first order ODE (that has to be solved in time) for the grid point with associated CDF or conserved mass and this point seamlessly moves in time as a function of above ODE when we solve it. I was able to solve the above ODE just as such (without any squared addition of volatilities separately) and got seamless results for the evolution of CDF(or probability mass up till  that point) conservation point in original coordinates.
When you evolve the SDE in original coordinates, you have to take a few steps with an alternative method (until second derivative is non-negligible) and then evolve with above method. I will be posting a program of simulation of densities of SDEs with above method in a few hours.

The above method is sister equation to continuity equation of mathematical physics. You can apply it to other first order equations or higher order equations. I am sure it can also be applied to Navier Stokes equation with some hack when some boundary is in motion and momentum is conserved. Continuity equation works when boundaries are given while this method works when boundary has to be solved.
.
.
.
I was thinking of ways how to increase the order in time of solution of the FPE. One idea that struck me and I want to work on is to take a further second derivative of equation
$\frac{\partial x_b}{\partial t} =\mu(x_b) - .5 \, \frac{\partial \big[{\sigma(x_b)}^2 \big]}{\partial x}+ .5\,{\sigma(x_b)}^2 \, Z_b \, \frac{\partial Z}{\partial x} + \, .5\,{\sigma(x_b)}^2 \, {(\frac{\partial Z}{\partial x})}^2 \frac{\partial^2 x}{\partial Z^2}$
with respect to time. On RHS we will have to take a further derivative with respect to x and multiply it (using chain rule of partial derivatives) with  $\frac{\partial x_b}{\partial t}$ which would have already  been calculated using the above equation.  Once we would have calculated $\frac{\partial^2 x_b}{\partial t^2}$ we can easily use it in a Taylor-like fashion to increase the order in time of the solution of FPE. Currently it is just a suggestion as I have not done it yet and cannot say for sure before hand  if there would be any unforeseen problems. If third derivative is problematic, we could possibly restrict our solution to terms containing just first two derivatives with respect to X but all of this would have to be carefully seen.

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

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

Friends, I noticed that our 1D SDE solution and 2D SDE solution both have instabilities close to zero. I worked on 1D SDE solution to improve its behavior close to zero. Sometimes SDE does not necessarily go to zero but its dynamics take extreme stress close to zero and it becomes unstable. I noticed when the solution is calculated to -5 SDs, there is extreme stress in parameters close to -5SD and there is always a chance of their becoming unstable. It was rather far more stable to simulate the SDE analytically to only -4.0 SDs and at the end (or when desired) extrapolate the solution to -4.5 SDs. when negative Z boundary is taken to only -4SDs, there is considerably less stress on SDE dynamics and many times it remains almost always stable even at very high volatilities. I made a few other modifications and so was able to simulate CIR SDE that goes into zero with moderate volatilities. If the volatility close to zero is too high, CIR SDE is very difficult to simulate. I have created some graphs to give friends some idea about zero behavior and accuracy of SDEs. In these graphs, CIR SDEs have volatility close to maximum possible meaning if I increased volatility somewhat more, my solution would become unstable while you could still increase volatility for SDEs with volatility exponents greater than .5. I will post my program tomorrow. I want to work on improving the zero behavior for SDEs in original coordinates and finally after that I will improve and make changes to 2D SDE simulation program. Here are the graphs of densities of some SDEs with improved zero behavior.
.
.
.

Please note that in the complete above graph, the tail goes out to 17-18. It is very close to lognormal with a specially elongated tail.

Please note that in the complete above graph, the tail goes out to 8. Its tail is much smaller than the tail with exponent .95 but still it is quite pronounced.

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

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

Friends, I have caught malaria and have very high fever. I will post the complete program after I recover. Here is the tentative program and I wanted to change a lot of things in it to make it easy to understand but now I am posting as such. I will post a reworked version once I recover.
.
.
function [] = FPERevisitedTransProb08ABwNew04DownloadedAB1()

%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=mu1 x(t)^beta1 dt + mu2 x(t)^beta2 dt +sigma x(t)^gamma dz(t)

%I have not directly simulated the SDE but simulated the transformed
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.

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

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

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

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

alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;%This is for expansion of integrals for calculation of drift
%and volatility coefficients
yy(1:Nn)=x0;
w(1:Nn)=x0^(1-gamma)/(1-gamma);

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

ZProbT(1)=normcdf(.5*ZT(1)+.5*ZT(2),0,1)-normcdf(.5*ZT(1)+.5*ZT(2)-dNn,0,1);
ZProbT(NnT)=normcdf(.5*ZT(NnT)+.5*ZT(NnT-1)+dNn,0,1)-normcdf(.5*ZT(NnT)+.5*ZT(NnT-1),0,1);
ZProbT(2:NnT-1)=normcdf(.5*ZT(2:NnT-1)+.5*ZT(3:NnT),0,1)-normcdf(.5*ZT(2:NnT-1)+.5*ZT(1:NnT-2),0,1);
%Above calculate probability mass in each probability subdivision.

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

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

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

wnStart=1;%
tic

Zt1(wnStart:Nn)=0.0;
Zt2(wnStart:Nn)=0.0;
for tt=1:Tt

%[wMu0dt,dwMu0dtdw,c1] = CalculateDriftAndVolA4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
%[wMu0dt,c1,c2] = CalculateDriftAndVolA404AB(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
[wMu0dt,c1,c2,dwMu0dtdw] = CalculateDriftAndVolA404AB(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);

[wMid] = InterpolateOrderN8(8,0,Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidl),Z(NnMidh),Z(NnMidh+1),Z(NnMidh+2),Z(NnMidh+3),w(NnMidl-3),w(NnMidl-2),w(NnMidl-1),w(NnMidl),w(NnMidh),w(NnMidh+1),w(NnMidh+2),w(NnMidh+3));

Zt1(wnStart:Nn)=w(wnStart:Nn)-wMid;

[dwdZ,d2wdZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt1,Z);
[dwdZ,d2wdZ2,d3wdZ3] = First3Derivatives4thOrderEqSpaced(wnStart,Nn,dNn,Zt1,Z);

tt

C0(wnStart:Nn)=Zt1(wnStart:Nn)-dwdZ(wnStart:Nn).*Z(wnStart:Nn);
%c1(wnStart:Nn)=sigma0*sqrt(dt);
Zt2(wnStart:Nn)=C0(wnStart:Nn)+sign(dwdZ(wnStart:Nn)+c1(wnStart:Nn)).* ...
abs(sqrt(sign(dwdZ(wnStart:Nn)).*(dwdZ(wnStart:Nn)).^2+ ...
sign(c1(wnStart:Nn)).*c1(wnStart:Nn).^2)).*Z(wnStart:Nn);
[dwdZ,d2wdZ2,d3wdZ3] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt2,Z);
[dw2dZ,d2w2dZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt2,Z);
[dw2dZ,d2w2dZ2A,d3wdZ3] = First3Derivatives4thOrderEqSpaced(wnStart,Nn,dNn,Zt2,Z);
if(tt>10)
%    dwdt(wnStart:Nn)=(wMu0dt(wnStart:Nn)+Zt2(wnStart:Nn)-Zt1(wnStart:Nn)+ ...
%        .5*sigma0^2*(dw2dZ(wnStart:Nn).^(-2)).*d2w2dZ2A(wnStart:Nn).*dt)/dt;
else
%    dwdt(wnStart:Nn)=(wMu0dt(wnStart:Nn)+Zt2(wnStart:Nn)-Zt1(wnStart:Nn))/dt;%+ ...
end

wnStart

if(tt>10)
dZdw(wnStart:Nn)=1./dw2dZ(wnStart:Nn);
d2Zdw2(wnStart:Nn)=-d2w2dZ2A(wnStart:Nn).*dZdw(wnStart:Nn).^3;
d3Zdw3(wnStart:Nn)=-d3wdZ3(wnStart:Nn).*dZdw(wnStart:Nn).^4+3*d2w2dZ2A(wnStart:Nn).^2.*dZdw(wnStart:Nn).^5;
%w(wnStart:Nn)=w(wnStart:Nn)+ ...
%    wMu0dt(wnStart:Nn)+ ...
%     0*dwMu0dtdw(wnStart:Nn)./(-Z(wnStart:Nn).*dZdw(wnStart:Nn).^2+d2Zdw2(wnStart:Nn)) - ...
%     +.5*sigma0^2*dt*( (Z(wnStart:Nn).^2-1).*dZdw(wnStart:Nn).^3- ...
%    0* 3*Z(wnStart:Nn).*dZdw(wnStart:Nn).*d2Zdw2(wnStart:Nn))./(-Z(wnStart:Nn).*dZdw(wnStart:Nn).^2+d2Zdw2(wnStart:Nn));

w(wnStart:Nn)=w(wnStart:Nn)+ ...
wMu0dt(wnStart:Nn)+1*(Zt2(wnStart:Nn)-Zt1(wnStart:Nn))+ ...
1*.5*sigma0^2*(dw2dZ(wnStart:Nn).^(-2)).*d2w2dZ2A(wnStart:Nn).*dt;
else
w(wnStart:Nn)=w(wnStart:Nn)+ ...
wMu0dt(wnStart:Nn)+Zt2(wnStart:Nn)-Zt1(wnStart:Nn);
end
%     %yy(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));

%  [dwdZ,d2wdZ2A,d3wdZ3] = First3Derivatives4thOrderEqSpaced(wnStart,Nn,dNn,w,Z);
%FirstCoordinate
%     for nn=NnMidl-11:-1:wnStart+1
%         if((w(nn)-w(nn-1))>(w(nn+1)-w(nn)))
%
%             w(nn) = InterpolateOrderN8(6,Z(nn),Z(nn+1),Z(nn+2),Z(nn+3),Z(nn+4),Z(nn+5),Z(nn+6),Z(nn+7),Z(nn+8),w(nn+1),w(nn+2),w(nn+3),w(nn+4),w(nn+5),w(nn+6),w(nn+7),w(nn+8));
%
%         end
%     end
%   if(w(wnStart+1)-w(wnStart)>w(wnStart+2)-w(wnStart+1))
%       nn=wnStart;
%       w(wnStart) = InterpolateOrderN8(6,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),Z(wnStart+8),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7),w(wnStart+8));

%    w(nn)=w(nn+1)-(w(nn+2)-w(nn+1));
%   end

[wE] = InterpolateOrderN6(6,Z(Nn)+dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));

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

%   w1(wnStart:Nn-1)=w(wnStart:Nn-1);
%   w1(Nn)=w(Nn);
%   w2(wnStart:Nn-1)=w(wnStart+1:Nn);
%   w2(Nn)=wE;
%   w(w1(wnStart:Nn)>w2(:))=0;%Be careful;might not universally hold;
%    %Change 3:7/25/2020: I have improved zero correction in above.

%for nn=1:Nn
%    if(w2(nn)>w1(nn))

w(w<0)=0.0;
StartChangeFlag=0;
nnChange=0;
wnStartPrev=wnStart;
for nn=wnStart:Nn
if(w(nn)<=0)
wnStart=nn+1;
StartChangeFlag=1;
nnChange=nnChange+1;
end
end
%nn=wnStart;
%while(w(nn)<=0)
%     wnStart=nn+1;
%     nn=nn+1;
%end

InterpolatePrevFlag=1;
wnChange=wnStart-wnStartPrev;
Zs=0
ws=0;
for mm=1:nnChange

Zs(mm)=Z(wnStart)-mm*dNn;
ws(mm)=0.0
end
if(nnChange>=1)
ws = InterpolateOrderN8(8,Zs,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
ContinueFlag=1;
end
for mm=1:nnChange

if((ws(mm)>0)&&(ws(mm)<w(wnStart))&&(ContinueFlag==1))
wnStart=wnStart-1;
w(wnStart)=ws(mm);
ContinueFlag=1;
else
ContinueFlag=0;
end
end
%         if(InterpolatePrevFlag==1)
%             InterpolatePrevFlag=0;
%             %wStart1 = InterpolateOrderN8(8,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
%             NnStart=nnChange
%             wStart1 = InterpolateOrderN8(8,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
%             if(wStart1>0)
%                 wnStart=wnStart-1;
%                 w(wnStart)=wStart1;
%                 InterpolatePrevFlag=1;
%             end
%             if(wnStart==wnStartPrev)
%             %    InterpolatePrevFlag=0;
%             end
%         end
%      end

%      if(nnChange>=1)
%          wStart1 = InterpolateOrderN8(8,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
%          if(wStart1>0)
%              wnStart=wnStart-1
%              w(wnStart)=wStart1;
%              InterpolatePrevFlag=1;
%          end
%      end
%       if((nnChange>=2)&&( InterpolatePrevFlag==1))
%          wStart1 = InterpolateOrderN8(8,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
%          InterpolatePrevFlag=0;
%          if(wStart1>0)
%              wnStart=wnStart-1
%              w(wnStart)=wStart1;
%              InterpolatePrevFlag=1;
%          end
%       end
%      if((nnChange>=3)&&( InterpolatePrevFlag==1))
%          wStart1 = InterpolateOrderN8(8,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
%          InterpolatePrevFlag=0;
%          if(wStart1>0)
%              wnStart=wnStart-1
%              w(wnStart)=wStart1;
%              InterpolatePrevFlag=1;
%          end
%      end
%
w
wnStart
tt
end

wT(wnStart+NnD:Nn+NnD)=w(wnStart:Nn);
Zs=0
ws=0;
for mm=1:wnStart+NnD-1
Zs(mm)=Z(wnStart)-mm*dNn;
ws(mm)=0.0
end
wnStartT=wnStart+NnD;
%NnT=Nn+2*NnD;
ws = InterpolateOrderN8(8,Zs,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
ContinueFlag=1;
wnStart0=wnStart+NnD-1;
for mm=1:wnStart0
if((ws(mm)>0)&&(ws(mm)<wT(wnStartT))&&(ContinueFlag==1))
wnStartT=wnStartT-1;
wT(wnStartT)=ws(mm);
ContinueFlag=1;
else
ContinueFlag=0;
end
end

%      Zs=0;
%      ws=0;
%      for mm=1:NnD
%         Zs(mm)=Z(Nn)+mm*dNn;
%         ws(mm)=0.0
%      end
%
%      ws = InterpolateOrderN8(2,Zs,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
%
%
%      %for mm=1:NnD
%
%          wT(Nn+NnD+1:Nn+2*NnD)=ws(1:NnD);
%      %end

%      for mm=1:wnStart0-1
%          if((ws(mm)>0)&&(ws(mm)<w(wnStart))&&(ContinueFlag==1))
%              wnStart=wnStart-1;
%              w(wnStart)=ws(mm);
%              ContinueFlag=1;
%          else
%              ContinueFlag=0;
%          end
%      end

wnStartT

%str=input('Look at numbers')
yyT(wnStartT:NnT)=((1-gamma).*wT(wnStartT:NnT)).^(1/(1-gamma));
DfyyT(wnStartT:NnT)=0;
for nn=wnStartT+1:NnT-1

DfyyT(nn) = (yyT(nn + 1) - yyT(nn - 1))/(ZT(nn + 1) - ZT(nn - 1));
%Change of variable derivative for densities
end

pyyT(1:Nn)=0;
for nn = wnStartT:NnT-1

pyyT(nn) = (normpdf(ZT(nn),0, 1))/abs(DfyyT(nn));
end

if(tt>=23)
%plot(yy(wnStart:Nn),pyy(wnStart:Nn),'r');
%yy
%str=input('Look at the output graph');
end

toc

ItoHermiteMean=sum(yyT(wnStartT+1:NnT-1).*ZProbT(wnStartT+1:NnT-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

theta1=1;
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.
disp('Original process average from our simulation');
ItoHermiteMean=sum(yyT(wnStartT+1:NnT-1).*ZProbT(wnStartT+1:NnT-1)) %Original process average from coordinates

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

MaxCutOff=30;
NoOfBins=round(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(yyT(wnStartT+1:NnT-1),pyyT(wnStartT+1:NnT-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');

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

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

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

end


.
.
function [wMu0dt,c1,c2,dwMu0dtdw] = CalculateDriftAndVolA404AB(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt)

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

Fp2=Fp1/(1-gamma);

wMu0dt(wnStart:Nn)=(YqCoeff0(1,1,2,1).*yy(wnStart:Nn).^Fp1(1,1,2,1)+ ...
YqCoeff0(1,2,1,1).*yy(wnStart:Nn).^Fp1(1,2,1,1)+ ...
YqCoeff0(2,1,1,1).*yy(wnStart:Nn).^Fp1(2,1,1,1))*dt + ...
(YqCoeff0(1,1,3,1).*yy(wnStart:Nn).^Fp1(1,1,3,1)+ ...
YqCoeff0(1,2,2,1).*yy(wnStart:Nn).^Fp1(1,2,2,1)+ ...
YqCoeff0(2,1,2,1).*yy(wnStart:Nn).^Fp1(2,1,2,1)+ ...
YqCoeff0(1,3,1,1).*yy(wnStart:Nn).^Fp1(1,3,1,1)+ ...
YqCoeff0(2,2,1,1).*yy(wnStart:Nn).^Fp1(2,2,1,1)+ ...
YqCoeff0(3,1,1,1).*yy(wnStart:Nn).^Fp1(3,1,1,1))*dt^2 + ...
(YqCoeff0(1,1,4,1).*yy(wnStart:Nn).^Fp1(1,1,4,1)+ ...
YqCoeff0(1,2,3,1).*yy(wnStart:Nn).^Fp1(1,2,3,1)+ ...
YqCoeff0(2,1,3,1).*yy(wnStart:Nn).^Fp1(2,1,3,1)+ ...
YqCoeff0(1,3,2,1).*yy(wnStart:Nn).^Fp1(1,3,2,1)+ ...
YqCoeff0(2,2,2,1).*yy(wnStart:Nn).^Fp1(2,2,2,1)+ ...
YqCoeff0(3,1,2,1).*yy(wnStart:Nn).^Fp1(3,1,2,1)+ ...
YqCoeff0(1,4,1,1).*yy(wnStart:Nn).^Fp1(1,4,1,1)+ ...
YqCoeff0(2,3,1,1).*yy(wnStart:Nn).^Fp1(2,3,1,1)+ ...
YqCoeff0(3,2,1,1).*yy(wnStart:Nn).^Fp1(3,2,1,1)+ ...
YqCoeff0(4,1,1,1).*yy(wnStart:Nn).^Fp1(4,1,1,1))*dt^3+ ...
(YqCoeff0(1,1,5,1).*yy(wnStart:Nn).^Fp1(1,1,5,1)+ ...
YqCoeff0(1,2,4,1).*yy(wnStart:Nn).^Fp1(1,2,4,1)+ ...
YqCoeff0(2,1,4,1).*yy(wnStart:Nn).^Fp1(2,1,4,1)+ ...
YqCoeff0(1,3,3,1).*yy(wnStart:Nn).^Fp1(1,3,3,1)+ ...
YqCoeff0(2,2,3,1).*yy(wnStart:Nn).^Fp1(2,2,3,1)+ ...
YqCoeff0(3,1,3,1).*yy(wnStart:Nn).^Fp1(3,1,3,1)+ ...
YqCoeff0(1,4,2,1).*yy(wnStart:Nn).^Fp1(1,4,2,1)+ ...
YqCoeff0(2,3,2,1).*yy(wnStart:Nn).^Fp1(2,3,2,1)+ ...
YqCoeff0(3,2,2,1).*yy(wnStart:Nn).^Fp1(3,2,2,1)+ ...
YqCoeff0(4,1,2,1).*yy(wnStart:Nn).^Fp1(4,1,2,1)+ ...
YqCoeff0(1,5,1,1).*yy(wnStart:Nn).^Fp1(1,5,1,1)+ ...
YqCoeff0(2,4,1,1).*yy(wnStart:Nn).^Fp1(2,4,1,1)+ ...
YqCoeff0(3,3,1,1).*yy(wnStart:Nn).^Fp1(3,3,1,1)+  ...
YqCoeff0(4,2,1,1).*yy(wnStart:Nn).^Fp1(4,2,1,1)+ ...
YqCoeff0(5,1,1,1).*yy(wnStart:Nn).^Fp1(5,1,1,1))*dt^4;

dwMu0dtdw(wnStart:Nn)=(YqCoeff0(1,1,2,1).*Fp1(1,1,2,1).*((1-gamma)*w(wnStart:Nn)).^(Fp2(1,1,2,1)-1)+ ...
YqCoeff0(1,2,1,1).*Fp1(1,2,1,1).*((1-gamma)*w(wnStart:Nn)).^(Fp2(1,2,1,1)-1)+ ...
YqCoeff0(2,1,1,1).*Fp1(2,1,1,1).*((1-gamma)*w(wnStart:Nn)).^(Fp2(2,1,1,1)-1))*dt ;

%wMu0dtOdt(wnStart:Nn)=(YqCoeff0(1,1,2,1).*yy(wnStart:Nn).^Fp1(1,1,2,1)*0+ ...
%    YqCoeff0(1,2,1,1).*yy(wnStart:Nn).^Fp1(1,2,1,1)+ ...
%    YqCoeff0(2,1,1,1).*yy(wnStart:Nn).^Fp1(2,1,1,1))*dt ;

c1(wnStart:Nn)=((YqCoeff0(1,1,1,2).*yy(wnStart:Nn).^Fp1(1,1,1,2).*sqrt(dt))+ ...
(YqCoeff0(1,1,2,2).*yy(wnStart:Nn).^Fp1(1,1,2,2)+YqCoeff0(1,2,1,2).*yy(wnStart:Nn).^Fp1(1,2,1,2)+ ...
YqCoeff0(2,1,1,2).*yy(wnStart:Nn).^Fp1(2,1,1,2)).*dt^1.5+ ...
(YqCoeff0(1,1,3,2).*yy(wnStart:Nn).^Fp1(1,1,3,2)+YqCoeff0(1,2,2,2).*yy(wnStart:Nn).^Fp1(1,2,2,2)+ ...
YqCoeff0(2,1,2,2).*yy(wnStart:Nn).^Fp1(2,1,2,2)+YqCoeff0(1,3,1,2).*yy(wnStart:Nn).^Fp1(1,3,1,2)+ ...
YqCoeff0(2,2,1,2).*yy(wnStart:Nn).^Fp1(2,2,1,2)+YqCoeff0(3,1,1,2).*yy(wnStart:Nn).^Fp1(3,1,1,2)).*dt^2.5+ ...
(YqCoeff0(1,1,4,2).*yy(wnStart:Nn).^Fp1(1,1,4,2)+YqCoeff0(1,2,3,2).*yy(wnStart:Nn).^Fp1(1,2,3,2)+ ...
YqCoeff0(2,1,3,2).*yy(wnStart:Nn).^Fp1(2,1,3,2)+YqCoeff0(1,3,2,2).*yy(wnStart:Nn).^Fp1(1,3,2,2)+ ...
YqCoeff0(2,2,2,2).*yy(wnStart:Nn).^Fp1(2,2,2,2)+ YqCoeff0(3,1,2,2).*yy(wnStart:Nn).^Fp1(3,1,2,2)+ ...
YqCoeff0(1,4,1,2).*yy(wnStart:Nn).^Fp1(1,4,1,2)+YqCoeff0(2,3,1,2).*yy(wnStart:Nn).^Fp1(2,3,1,2)+ ...
YqCoeff0(3,2,1,2).*yy(wnStart:Nn).^Fp1(3,2,1,2)+YqCoeff0(4,1,1,2).*yy(wnStart:Nn).^Fp1(4,1,1,2)).*dt^3.5);

c2(wnStart:Nn)=(YqCoeff0(1,1,2,3).*yy(wnStart:Nn).^Fp1(1,1,2,3)+YqCoeff0(1,2,1,3).*yy(wnStart:Nn).^Fp1(1,2,1,3)+ ...
YqCoeff0(2,1,1,3).*yy(wnStart:Nn).^Fp1(2,1,1,3)).*dt^2;%+ ...
(YqCoeff0(1,1,3,3).*yy(wnStart:Nn).^Fp1(1,1,3,3)+YqCoeff0(1,2,2,3).*yy(wnStart:Nn).^Fp1(1,2,2,3)+ ...
YqCoeff0(2,1,2,3).*yy(wnStart:Nn).^Fp1(2,1,2,3) + YqCoeff0(1,3,1,3).*yy(wnStart:Nn).^Fp1(1,3,1,3)+ ...
YqCoeff0(2,2,1,3).*yy(wnStart:Nn).^Fp1(2,2,1,3)+YqCoeff0(3,1,1,3).*yy(wnStart:Nn).^Fp1(3,1,1,3)).*dt^3;

% (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+ ...
%     (YCoeff0(1,1,3,2).*YY(1:paths).^Fp(1,1,3,2)+YCoeff0(1,2,2,2).*YY(1:paths).^Fp(1,2,2,2)+ ...
%     YCoeff0(2,1,2,2).*YY(1:paths).^Fp(2,1,2,2)+YCoeff0(1,3,1,2).*YY(1:paths).^Fp(1,3,1,2)+ ...
%     YCoeff0(2,2,1,2).*YY(1:paths).^Fp(2,2,1,2)+YCoeff0(3,1,1,2).*YY(1:paths).^Fp(3,1,1,2)).*dtM^2.5+ ...
%     (YCoeff0(1,1,4,2).*YY(1:paths).^Fp(1,1,4,2)+YCoeff0(1,2,3,2).*YY(1:paths).^Fp(1,2,3,2)+ ...
%     YCoeff0(2,1,3,2).*YY(1:paths).^Fp(2,1,3,2)+YCoeff0(1,3,2,2).*YY(1:paths).^Fp(1,3,2,2)+ ...
%     YCoeff0(2,2,2,2).*YY(1:paths).^Fp(2,2,2,2)+ YCoeff0(3,1,2,2).*YY(1:paths).^Fp(3,1,2,2)+ ...
%     YCoeff0(1,4,1,2).*YY(1:paths).^Fp(1,4,1,2)+YCoeff0(2,3,1,2).*YY(1:paths).^Fp(2,3,1,2)+ ...
%     YCoeff0(3,2,1,2).*YY(1:paths).^Fp(3,2,1,2)+YCoeff0(4,1,1,2).*YY(1:paths).^Fp(4,1,1,2)).*dtM^3.5) .*HermiteP1(2,1:paths)

end
.
.
.
function [dwdZ,d2wdZ2,d3wdZ3] = First3Derivatives4thOrderEqSpaced(wnStart,Nn,dNn,w,Z)

% [wS] = InterpolateOrderN6(4,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5));
% [wE] = InterpolateOrderN6(4,Z(Nn)+dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));

% [wS] = InterpolateOrderN8(7,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
% [wE] = InterpolateOrderN8(7,Z(Nn)+dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),Z(Nn-6),Z(Nn-7),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5),w(Nn-6),w(Nn-7));

[wS] = InterpolateOrderN6(6,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5));
[wE] = InterpolateOrderN6(6,Z(Nn)+dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));

[wS2] = InterpolateOrderN6(6,Z(wnStart)-2*dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5));
[wE2] = InterpolateOrderN6(6,Z(Nn)+2*dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));

%  Zi=[Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),Z(wnStart+8),Z(wnStart+9)];
%       wi=[w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7),w(wnStart+8),w(wnStart+9)];
% %      Zi=[Z(Nn-6),Z(Nn-5),Z(Nn-4),Z(Nn-3),Z(Nn-2),Z(Nn-1),Z(Nn)];
% %      wi=[w(Nn-6),w(Nn-5),w(Nn-4),w(Nn-3),w(Nn-2),w(Nn-1),w(Nn)];
%
%       Zi=[Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6)];
%       wi=[w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6)];
%       Zq=Z(wnStart)-dNn;
%       wS=spline(Zi,wi,Zq);
%
%       Zq=Z(wnStart)-2*dNn;
%       wS2=spline(Zi,wi,Zq);
%
%
%
% %      Zi=[Z(Nn-9),Z(Nn-8),Z(Nn-7),Z(Nn-6),Z(Nn-5),Z(Nn-4),Z(Nn-3),Z(Nn-2),Z(Nn-1),Z(Nn)];
% %      wi=[w(Nn-9),w(Nn-8),w(Nn-7),w(Nn-6),w(Nn-5),w(Nn-4),w(Nn-3),w(Nn-2),w(Nn-1),w(Nn)];
%       Zi=[Z(Nn-6),Z(Nn-5),Z(Nn-4),Z(Nn-3),Z(Nn-2),Z(Nn-1),Z(Nn)];
%       wi=[w(Nn-6),w(Nn-5),w(Nn-4),w(Nn-3),w(Nn-2),w(Nn-1),w(Nn)];
%       Zq=Z(Nn)+dNn;
%       wE=spline(Zi,wi,Zq);
%
%
%       Zq=Z(Nn)+2*dNn;
%       wE2=spline(Zi,wi,Zq);

ZS=Z(wnStart)-dNn;
ZE=Z(Nn)+dNn;

% d2wdZ2(wnStart)=(Z(wnStart)*wS-2*w(wnStart)*Z(wnStart)+w(wnStart+1).*Z(wnStart))/(dNn.^2);
% d2wdZ2(wnStart+1:Nn-1)=(w(wnStart+1:Nn-1).*Z(wnStart+1:Nn-1)-2*w(wnStart+1:Nn-1).*Z(wnStart+1:Nn-1)+1*w(wnStart+2:Nn).*Z(wnStart+1:Nn-1))/(dNn.^2);
% d2wdZ2(Nn)=(wE.*Z(Nn)-2*w(Nn).*Z(Nn)+w(Nn-1).*Z(Nn))/(dNn.^2);

dwdZ(wnStart+2:Nn-2)=(1*w(wnStart:Nn-4)-8*w(wnStart+1:Nn-3)+0*w(wnStart+2:Nn-2)+8*w(wnStart+3:Nn-1)-1*w(wnStart+4:Nn))/(12*dNn);
dwdZ(wnStart+1)=(-3*w(wnStart)-10*w(wnStart+1)+18*w(wnStart+2)-6*w(wnStart+3)+1*w(wnStart+4))/(12*dNn);
dwdZ(wnStart)=(-25/12*w(wnStart)+4*w(wnStart+1)-3*w(wnStart+2)+4/3*w(wnStart+3)-1/4*w(wnStart+4))/dNn;%/(12*dNn);
dwdZ(Nn-1)=(-1*w(Nn-4)+6*w(Nn-3)-18*w(Nn-2)+10*w(Nn-1)+3*w(Nn))/(12*dNn);
dwdZ(Nn)=(1/4*w(Nn-4)-4/3*w(Nn-3)+3*w(Nn-2)-4*w(Nn-1)+25/12*w(Nn))/dNn;%/(12*dNn);

%−25/12	4	−3	4/3	−1/4
%(-3*f[i-1]-10*f[i+0]+18*f[i+1]-6*f[i+2]+1*f[i+3])/(12*1.0*h**1)
% (-1*f[i-3]+6*f[i-2]-18*f[i-1]+10*f[i+0]+3*f[i+1])/(12*1.0*h**1)
d2wdZ2(wnStart+2:Nn-2)=(-1*w(wnStart:Nn-4)+16*w(wnStart+1:Nn-3)-30*w(wnStart+2:Nn-2)+16*w(wnStart+3:Nn-1)-1*w(wnStart+4:Nn))/(12*dNn^2);
d2wdZ2(wnStart+1)=(10*w(wnStart)-15*w(wnStart+1)-4*w(wnStart+2)+14*w(wnStart+3)-6*w(wnStart+4)+1*w(wnStart+5))/(12*dNn^2);
d2wdZ2(wnStart)=(15/4*w(wnStart)-77/6*w(wnStart+1)+107/6*w(wnStart+2)-13*w(wnStart+3)+61/12*w(wnStart+4)-5/6*w(wnStart+5))/(dNn^2);
d2wdZ2(Nn-1)=(10*w(Nn)-15*w(Nn-1)-4*w(Nn-2)+14*w(Nn-3)-6*w(Nn-4)+1*w(Nn-5))/(12*dNn^2);
%d2wdZ2(Nn)=(-1*w(Nn-2)+16*w(Nn-1)-30*w(Nn)+16*wE-1*wE2)/(12*dNn^2);
d2wdZ2(Nn)=(15/4*w(Nn)-77/6*w(Nn-1)+107/6*w(Nn-2)-13*w(Nn-3)+61/12*w(Nn-4)-5/6*w(Nn-5))/(dNn^2);

% (1*f[i-4]-6*f[i-3]+14*f[i-2]-4*f[i-1]-15*f[i+0]+10*f[i+1])/(12*1.0*h**2)
% (10*f[i-1]-15*f[i+0]-4*f[i+1]+14*f[i+2]-6*f[i+3]+1*f[i+4])/(12*1.0*h**2)
% 15/4	−77/6	107/6	−13	61/12	−5/6

%
%
%  d2wdZ2(wnStart)=(wS-2*w(wnStart)+w(wnStart+1))/(dNn.^2);
%  d2wdZ2(wnStart+1:Nn-1)=(w(wnStart:Nn-2)-2*w(wnStart+1:Nn-1)+1*w(wnStart+2:Nn))/(dNn.^2);
%  d2wdZ2(Nn)=(wE-2*w(Nn)+w(Nn-1))/(dNn.^2);
%
%
%  dwdZ(wnStart)=(-1*wS+1*w(wnStart+1))/(2*dNn);
%  dwdZ(wnStart+1:Nn-1)=(-1*w(wnStart:Nn-2)+1*w(wnStart+2:Nn))/(2*dNn);
%  dwdZ(Nn)=(wE-w(Nn-1))/(2*dNn);
%
d3wdZ3(wnStart+2:Nn-2)=(-1*w(wnStart:Nn-4)+2*w(wnStart+1:Nn-3)+0*w(wnStart+2:Nn-2)-2*w(wnStart+3:Nn-1)+w(wnStart+4:Nn))/(2*dNn^3);
d3wdZ3(wnStart+1)=(-1*wS+2*w(wnStart)+0*w(wnStart+1)-2*w(wnStart+2)+w(wnStart+3))/(2*dNn^3);
d3wdZ3(wnStart)=(-1*wS2+2*wS+0*w(wnStart)-2*w(wnStart+1)+w(wnStart+2))/(2*dNn^3);
d3wdZ3(Nn-1)=(-1*w(Nn-3)+2*w(Nn-2)+0*w(Nn-1)-2*w(Nn)+wE)/(2*dNn^3);
d3wdZ3(Nn)=(-1*w(Nn-2)+2*w(Nn-1)+0*w(Nn)-2*wE+wE2)/(2*dNn^3);

end

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

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

Made some very minor changes to the program.
.
function [] = FPERevisitedTransProb08ABwNew04DownloadedAB1()

%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=mu1 x(t)^beta1 dt + mu2 x(t)^beta2 dt +sigma x(t)^gamma dz(t)

%I have not directly simulated the SDE but simulated the transformed
%Besse1l process version of the SDE and then changed coordinates to retreive
%the SDE in original coo
%rdinates.
%The present program will analytically evolve only the Bessel Process version of the
%SDE in transformed coordinates.

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

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

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

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

alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.
alpha1=1-gamma;%This is for expansion of integrals for calculation of drift
%and volatility coefficients
yy(1:Nn)=x0;
w(1:Nn)=x0^(1-gamma)/(1-gamma);

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

ZProbT(1)=normcdf(.5*ZT(1)+.5*ZT(2),0,1)-normcdf(.5*ZT(1)+.5*ZT(2)-dNn,0,1);
ZProbT(NnT)=normcdf(.5*ZT(NnT)+.5*ZT(NnT-1)+dNn,0,1)-normcdf(.5*ZT(NnT)+.5*ZT(NnT-1),0,1);
ZProbT(2:NnT-1)=normcdf(.5*ZT(2:NnT-1)+.5*ZT(3:NnT),0,1)-normcdf(.5*ZT(2:NnT-1)+.5*ZT(1:NnT-2),0,1);
%Above calculate probability mass in each probability subdivision.

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

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

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

wnStart=1;%
tic

Zt1(wnStart:Nn)=0.0;
Zt2(wnStart:Nn)=0.0;
for tt=1:Tt

%[wMu0dt,dwMu0dtdw,c1] = CalculateDriftAndVolA4(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
%[wMu0dt,c1,c2] = CalculateDriftAndVolA404AB(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);
[wMu0dt,c1,c2,dwMu0dtdw] = CalculateDriftAndVolA404AB(w,wnStart,Nn,YqCoeff0,Fp1,gamma,dt);

[wMid] = InterpolateOrderN8(8,0,Z(NnMidl-3),Z(NnMidl-2),Z(NnMidl-1),Z(NnMidl),Z(NnMidh),Z(NnMidh+1),Z(NnMidh+2),Z(NnMidh+3),w(NnMidl-3),w(NnMidl-2),w(NnMidl-1),w(NnMidl),w(NnMidh),w(NnMidh+1),w(NnMidh+2),w(NnMidh+3));

Zt1(wnStart:Nn)=w(wnStart:Nn)-wMid;

%  [dwdZ,d2wdZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt1,Z);
[dwdZ,d2wdZ2,d3wdZ3] = First3Derivatives4thOrderEqSpaced(wnStart,Nn,dNn,Zt1,Z);

tt

C0(wnStart:Nn)=Zt1(wnStart:Nn)-dwdZ(wnStart:Nn).*Z(wnStart:Nn);
%c1(wnStart:Nn)=sigma0*sqrt(dt);
Zt2(wnStart:Nn)=C0(wnStart:Nn)+sign(dwdZ(wnStart:Nn)+c1(wnStart:Nn)).* ...
abs(sqrt(sign(dwdZ(wnStart:Nn)).*(dwdZ(wnStart:Nn)).^2+ ...
sign(c1(wnStart:Nn)).*c1(wnStart:Nn).^2)).*Z(wnStart:Nn);
%   [dwdZ,d2wdZ2,d3wdZ3] = First3Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt2,Z);
%   [dw2dZ,d2w2dZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt2,Z);
[dw2dZ,d2w2dZ2A,d3wdZ3] = First3Derivatives4thOrderEqSpaced(wnStart,Nn,dNn,Zt2,Z);
if(tt>10)
%    dwdt(wnStart:Nn)=(wMu0dt(wnStart:Nn)+Zt2(wnStart:Nn)-Zt1(wnStart:Nn)+ ...
%        .5*sigma0^2*(dw2dZ(wnStart:Nn).^(-2)).*d2w2dZ2A(wnStart:Nn).*dt)/dt;
else
%    dwdt(wnStart:Nn)=(wMu0dt(wnStart:Nn)+Zt2(wnStart:Nn)-Zt1(wnStart:Nn))/dt;%+ ...
end

wnStart

if(tt>10)
dZdw(wnStart:Nn)=1./dw2dZ(wnStart:Nn);
d2Zdw2(wnStart:Nn)=-d2w2dZ2A(wnStart:Nn).*dZdw(wnStart:Nn).^3;
d3Zdw3(wnStart:Nn)=-d3wdZ3(wnStart:Nn).*dZdw(wnStart:Nn).^4+3*d2w2dZ2A(wnStart:Nn).^2.*dZdw(wnStart:Nn).^5;
%w(wnStart:Nn)=w(wnStart:Nn)+ ...
%    wMu0dt(wnStart:Nn)+ ...
%     0*dwMu0dtdw(wnStart:Nn)./(-Z(wnStart:Nn).*dZdw(wnStart:Nn).^2+d2Zdw2(wnStart:Nn)) - ...
%     +.5*sigma0^2*dt*( (Z(wnStart:Nn).^2-1).*dZdw(wnStart:Nn).^3- ...
%    0* 3*Z(wnStart:Nn).*dZdw(wnStart:Nn).*d2Zdw2(wnStart:Nn))./(-Z(wnStart:Nn).*dZdw(wnStart:Nn).^2+d2Zdw2(wnStart:Nn));

w(wnStart:Nn)=w(wnStart:Nn)+ ...
wMu0dt(wnStart:Nn)+1*(Zt2(wnStart:Nn)-Zt1(wnStart:Nn))+ ...
1*.5*sigma0^2*(dw2dZ(wnStart:Nn).^(-2)).*d2w2dZ2A(wnStart:Nn).*dt;
else
w(wnStart:Nn)=w(wnStart:Nn)+ ...
wMu0dt(wnStart:Nn)+Zt2(wnStart:Nn)-Zt1(wnStart:Nn);
end
%     %yy(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));

%  [dwdZ,d2wdZ2A,d3wdZ3] = First3Derivatives4thOrderEqSpaced(wnStart,Nn,dNn,w,Z);
%FirstCoordinate
%     for nn=NnMidl-11:-1:wnStart+1
%         if((w(nn)-w(nn-1))>(w(nn+1)-w(nn)))
%
%             w(nn) = InterpolateOrderN8(6,Z(nn),Z(nn+1),Z(nn+2),Z(nn+3),Z(nn+4),Z(nn+5),Z(nn+6),Z(nn+7),Z(nn+8),w(nn+1),w(nn+2),w(nn+3),w(nn+4),w(nn+5),w(nn+6),w(nn+7),w(nn+8));
%
%         end
%     end
%   if(w(wnStart+1)-w(wnStart)>w(wnStart+2)-w(wnStart+1))
%       nn=wnStart;
%       w(wnStart) = InterpolateOrderN8(6,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),Z(wnStart+8),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7),w(wnStart+8));

%    w(nn)=w(nn+1)-(w(nn+2)-w(nn+1));
%   end

[wE] = InterpolateOrderN6(6,Z(Nn)+dNn,Z(Nn),Z(Nn-1),Z(Nn-2),Z(Nn-3),Z(Nn-4),Z(Nn-5),w(Nn),w(Nn-1),w(Nn-2),w(Nn-3),w(Nn-4),w(Nn-5));

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

%   w1(wnStart:Nn-1)=w(wnStart:Nn-1);
%   w1(Nn)=w(Nn);
%   w2(wnStart:Nn-1)=w(wnStart+1:Nn);
%   w2(Nn)=wE;
%   w(w1(wnStart:Nn)>w2(:))=0;%Be careful;might not universally hold;
%    %Change 3:7/25/2020: I have improved zero correction in above.

%for nn=1:Nn
%    if(w2(nn)>w1(nn))

w(w<0)=0.0;
StartChangeFlag=0;
nnChange=0;
wnStartPrev=wnStart;
for nn=wnStart:Nn
if(w(nn)<=0)
wnStart=nn+1;
StartChangeFlag=1;
nnChange=nnChange+1;
end
end
%nn=wnStart;
%while(w(nn)<=0)
%     wnStart=nn+1;
%     nn=nn+1;
%end

InterpolatePrevFlag=1;
wnChange=wnStart-wnStartPrev;
Zs=0
ws=0;
for mm=1:nnChange

Zs(mm)=Z(wnStart)-mm*dNn;
ws(mm)=0.0
end
if(nnChange>=1)
ws = InterpolateOrderN8(8,Zs,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
ContinueFlag=1;
end
for mm=1:nnChange

if((ws(mm)>0)&&(ws(mm)<w(wnStart))&&(ContinueFlag==1))
wnStart=wnStart-1;
w(wnStart)=ws(mm);
ContinueFlag=1;
else
ContinueFlag=0;
end
end
%         if(InterpolatePrevFlag==1)
%             InterpolatePrevFlag=0;
%             %wStart1 = InterpolateOrderN8(8,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
%             NnStart=nnChange
%             wStart1 = InterpolateOrderN8(8,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
%             if(wStart1>0)
%                 wnStart=wnStart-1;
%                 w(wnStart)=wStart1;
%                 InterpolatePrevFlag=1;
%             end
%             if(wnStart==wnStartPrev)
%             %    InterpolatePrevFlag=0;
%             end
%         end
%      end

%      if(nnChange>=1)
%          wStart1 = InterpolateOrderN8(8,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
%          if(wStart1>0)
%              wnStart=wnStart-1
%              w(wnStart)=wStart1;
%              InterpolatePrevFlag=1;
%          end
%      end
%       if((nnChange>=2)&&( InterpolatePrevFlag==1))
%          wStart1 = InterpolateOrderN8(8,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
%          InterpolatePrevFlag=0;
%          if(wStart1>0)
%              wnStart=wnStart-1
%              w(wnStart)=wStart1;
%              InterpolatePrevFlag=1;
%          end
%       end
%      if((nnChange>=3)&&( InterpolatePrevFlag==1))
%          wStart1 = InterpolateOrderN8(8,Z(wnStart)-dNn,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
%          InterpolatePrevFlag=0;
%          if(wStart1>0)
%              wnStart=wnStart-1
%              w(wnStart)=wStart1;
%              InterpolatePrevFlag=1;
%          end
%      end
%
w
wnStart
tt
end

wT(wnStart+NnD:Nn+NnD)=w(wnStart:Nn);
Zs=0
ws=0;
for mm=1:wnStart+NnD-1
Zs(mm)=Z(wnStart)-mm*dNn;
ws(mm)=0.0
end
wnStartT=wnStart+NnD;
%NnT=Nn+2*NnD;
ws = InterpolateOrderN8(8,Zs,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
ContinueFlag=1;
wnStart0=wnStart+NnD-1;
for mm=1:wnStart0
if((ws(mm)>0)&&(ws(mm)<wT(wnStartT))&&(ContinueFlag==1))
wnStartT=wnStartT-1;
wT(wnStartT)=ws(mm);
ContinueFlag=1;
else
ContinueFlag=0;
end
end

%      Zs=0;
%      ws=0;
%      for mm=1:NnD
%         Zs(mm)=Z(Nn)+mm*dNn;
%         ws(mm)=0.0
%      end
%
%      ws = InterpolateOrderN8(2,Zs,Z(wnStart),Z(wnStart+1),Z(wnStart+2),Z(wnStart+3),Z(wnStart+4),Z(wnStart+5),Z(wnStart+6),Z(wnStart+7),w(wnStart),w(wnStart+1),w(wnStart+2),w(wnStart+3),w(wnStart+4),w(wnStart+5),w(wnStart+6),w(wnStart+7));
%
%
%      %for mm=1:NnD
%
%          wT(Nn+NnD+1:Nn+2*NnD)=ws(1:NnD);
%      %end

%      for mm=1:wnStart0-1
%          if((ws(mm)>0)&&(ws(mm)<w(wnStart))&&(ContinueFlag==1))
%              wnStart=wnStart-1;
%              w(wnStart)=ws(mm);
%              ContinueFlag=1;
%          else
%              ContinueFlag=0;
%          end
%      end

wnStartT

%str=input('Look at numbers')
yyT(wnStartT:NnT)=((1-gamma).*wT(wnStartT:NnT)).^(1/(1-gamma));
DfyyT(wnStartT:NnT)=0;
for nn=wnStartT+1:NnT-1

DfyyT(nn) = (yyT(nn + 1) - yyT(nn - 1))/(ZT(nn + 1) - ZT(nn - 1));
%Change of variable derivative for densities
end

pyyT(1:Nn)=0;
for nn = wnStartT:NnT-1

pyyT(nn) = (normpdf(ZT(nn),0, 1))/abs(DfyyT(nn));
end

if(tt>=23)
%plot(yy(wnStart:Nn),pyy(wnStart:Nn),'r');
%yy
%str=input('Look at the output graph');
end

toc

ItoHermiteMean=sum(yyT(wnStartT+1:NnT-1).*ZProbT(wnStartT+1:NnT-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

theta1=1;
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.
disp('Original process average from our simulation');
ItoHermiteMean=sum(yyT(wnStartT+1:NnT-1).*ZProbT(wnStartT+1:NnT-1)) %Original process average from coordinates

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

MaxCutOff=30;
NoOfBins=round(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(yyT(wnStartT+1:NnT-1),pyyT(wnStartT+1:NnT-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');
%plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g',Z(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'b');

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

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

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

end



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

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

Friends I have malaria and very high fever, they still continue to charge my body and release sickening gases regularly in my room.  Mind control crooks have absolutely no humanity. Please protest against such inhuman torture on innocent civilians on hands of mind control crooks.

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

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

Friends, I have recovered from disease but have extreme weakness. I hope to restart my research work in 2-3 days. I want to thank friends who wished for my illness to end.
It was very difficult to fight the disease and also fight the mind control torture. I made a list of mind control activities of  agents and will share that with friends later.