Serving the Quantitative Finance Community

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

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

Friends here are a few notes about my derivation in previous post. I am sure most friends already are aware of this but still writing.
1. you can easily substitute Fokker-planck equation in Bessel and original coordinates to get a nice solution to moving boundary of the conserved probability mass.
2. What we derived in previous post is analogous to being sister equation of continuity equation of mathematical physics but just the one with a moving boundary. And I am very sure it can be used at a large number of places where some quantity(not just probability mass) is being conserved and a boundary is moving to get the solution to moving boundary.
Friends, I love sharing my research here but I want to request all the good people to please give me credit for my original research even if I do not write a formal paper on it though I might as well write a formal paper in next few weeks or months in Wilmott journal.
I hope to be coming back soon with new posts with solution to two dimensional system of stochastic volatility type SDEs.
I hope that my not writing a paper or keeping a low key/profile would never be a reason for friends to not give me credit for my original research. I have always stated that once my research is in public on internet, everyone is free to do their further research and extensions but they have to give me credit for my original research. I think friends would agree to my fair request.
I was after this problem since 2014.

Amin
Topic Author
Posts: 2906
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.

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

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

With a heavy heart, I want to tell all good Americans that my persecution continues unabated. There are marks and light swelling on the body tissue on my back due to high frequency waves constantly being targeted to go inside from my back. This is a constant torture and keeps target people irritated and is very unsettling for them but lowly animals of American army love such tactics since it becomes extremely difficult for the target to concentrate on anything as the target is kept irritated.
Yesterday, they drugged water bottles at a bakery where I went to take food and things like this this continue every day.
I have been begging for help for years and years but to no avail and the torture by lowly crooks never ends at all. These lowly racist crooks who are incapable of doing any good in this world want to stop every black, muslim or foreigner from doing anything productive in this world. Is there any justice in this world or the these racist crooks are the last words of justice? Please help me.

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

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

I made this post a few months ago about inviting prominent investigative journalists to end the ill of mind control from human societies. Every western newspaper runs large stories about persecution of people in Xinjiang which is a great thing to do but why does these newspapers become tight-lipped when it comes to systematic torture by American army on thousands of talented people of other color, religion, or nationality. Please speak vocally to end this torture now before it is too late.

Here is my previous post:

It seems that crooks in American army and mind control  are unable to pay a heed to civilized calls by good American people and others to end their animal cruel practices and remain bent on mischief. If crooks in American army and mind control remain adamant on continuing evil practices, the only solution seems that some brave investigative journalists expose them by running a comprehensive story on mind control practices by American army. It seems that doing a civilized dialogue with hardened crooks in mind control gives them a strong sense of weakness of good people and makes the crooks even more adamant to continue their evil practices.  Reminds me of  Abu-Gharib. Very similarly, Crooks in defense had no conception that they were doing any wrong thing when there was a culture of openly urinating on human captives. It was an open thing in army and most crooks thought it was indeed a very right thing to do( and I am sure some of those crooks still believe that it was a right thing they did even after being reprimanded by the broader American society). It was not until brave people at CNN did a daring story against animal practices by crooks that evil practices stopped. Though American army crooks retard intelligent people of all color and creed, these ultra right wing army crooks love to retard blacks and muslims with great relish. Blacks are only 8-10% of US population but make a very large proportion of victims in united states since many ultra right wing crooks in US army would rather die than let those intelligent blacks succeed in American society in a big way.
I want to request CNN, New York Times, Washington Post and large reputed European media outlets to run a comprehensive story on mind control torture, retarding and victimization of intelligent people of US and other nations by crooks in US army on behest of some powerful people in United States and due to rightwing extremist biases of these crooks. If you would like to do investigative research about animal practices of US army, one great resource would be mainland European embassies in Muslim countries who keep a detailed account of mind control persecution of intelligent muslims in these countries by crooks in US army. Many of the staff in mainland European embassies are very good human beings who abhor such practices and would love to cooperate with good journalists in exposing the evil animal practices of US army crooks. Only the accounts of people at mainland European embassies in muslim countries thorough animal practices of American army's mind control wing would be enough to drop a great bombshell in the media and general public all across the world. I want to warn all the good people who try to have a civilized dialogue with crooks in American army to end their evil practices that their being nice with crooks is a very misguided approach. Since good people do not have the power to forcefully end the evil practices of US army crooks, only way to end the evil practices of US army crooks is by exposing them openly in US public and all across the world.

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

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

Friends, here is the program with ODE evolution of densities of SDEs in original coordinates. I have done first 15 steps simulation in bessel coordinates since second derivative would be too small for stable evolution in original coordinates. Zero boundary is still not covered and friends would have to work on that a bit carefully. For densities that go into zero, program is still not appropriate in this basic form.
Again in this program, I have used the ODE for evolution of constant CDF points on the grid just as the ODE is derived in previous few posts.

Here is the program.
.
function [] = FPERevisitedTransProb08ABwNew04B()

%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/2;   % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=128*2*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=dt*4*2*4;
TtM=Tt/4/2/4;

dNn=.2/1;   % Normal density subdivisions width. would change with number of subdivisions
Nn=50;  % No of normal density subdivisions
NnMidl=25;%One half density Subdivision left from mid of normal density(low)
NnMidh=26;%One half density subdivision right from the mid of normal density(high)
NnMid=4.0;

x0=1.0;   % starting value of SDE
beta1=0.0;
beta2=1.0;   % Second drift term power.
gamma=.95;%50;   % volatility power.
kappa=4.0;%.950;   %mean reversion parameter.
theta=.250;%mean reversion target
sigma0=1.50;%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);
x(1:Nn)=x0;

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

Z
str=input('Look at Z');
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1)-normcdf(.5*Z(1)+.5*Z(2)-dNn,0,1);
ZProb(Nn)=normcdf(.5*Z(Nn)+.5*Z(Nn-1)+dNn,0,1)-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);
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.

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

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

x(isnan(x)==1)=0.00;

[xMu0dt,c1] = CalculateDriftAndVolA4Original(x,wnStart,Nn,YCoeff0,Fp,gamma,dt);

%Loop below tackles first ten time steps in Bessel coordinates
if(tt<=10)
w(wnStart:Nn)=x(wnStart:Nn).^(1-gamma)/(1-gamma);
[wMu0dt,dwMu0dtdw,c1] = CalculateDriftAndVolA4(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);

C0(wnStart:Nn)=Zt1(wnStart:Nn)-dwdZ(wnStart:Nn).*Z(wnStart:Nn);

Zt2(wnStart:Nn)=C0(wnStart:Nn)+abs(sqrt((dwdZ(wnStart:Nn)).^2+sigma0^2*dt)).*Z(wnStart:Nn);

%x2(wnStart:Nn)=((1-gamma)*(Zt2(wnStart:Nn)-0*Zt1(wnStart:Nn)+wMid)).^(1.0/(1-gamma));

end
[dxdZ,d2xdZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,x,Z);

%     [xMid] = 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),x(NnMidl-3),x(NnMidl-2),x(NnMidl-1),x(NnMidl),x(NnMidh),x(NnMidh+1),x(NnMidh+2),x(NnMidh+3));
%
%     Zt1(wnStart:Nn)=x(wnStart:Nn)-xMid;
%
%     [dxdZ,d2xdZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,x,Z);
%
%     C0(wnStart:Nn)=Zt1(wnStart:Nn)-dxdZ(wnStart:Nn).*Z(wnStart:Nn);%-.5*d2xdZ2A(wnStart:Nn).*(Z(wnStart:Nn).^2-1);
%
%     %Zt2(wnStart:Nn)=C0(wnStart:Nn)+abs(sqrt((dxdZ(wnStart:Nn)).^2+sigma0^2*x(wnStart:Nn).^(2*gamma).*dt)).*Z(wnStart:Nn);%+.5.*(sqrt((d2xdZ2A(wnStart:Nn).^2)+(sigma0^2.*(gamma).*x(wnStart:Nn).^(2*gamma-1)*dt).^2)).*(Z(wnStart:Nn).^2-1);
%
%
%     %Zt2(wnStart:Nn)=C0(wnStart:Nn)+abs(sqrt((dxdZ(wnStart:Nn).*x(wnStart:Nn).^gamma).^2+c1(wnStart:Nn).^2)).*Z(wnStart:Nn);
%     [dxdZ,d2xdZ2A] = First2Derivatives2ndOrderEqSpacedA(wnStart,Nn,dNn,Zt2,Z);
%     x2=Zt2(wnStart:Nn)+xMid;

%    tt
if(tt>10)

x(wnStart:Nn)= +x(wnStart:Nn)-sigma0.^2*gamma*x(wnStart:Nn).^(2*gamma-1)*dt ...
+xMu0dt(wnStart:Nn) ...
+.5*sigma0^2*x(wnStart:Nn).^(2*gamma).*(dxdZ(wnStart:Nn)).^(-2).*(d2xdZ2A(wnStart:Nn)).*dt ...
+.5*sigma0^2*x(wnStart:Nn).^(2*gamma).*Z(wnStart:Nn).*(dxdZ(wnStart:Nn)).^(-1)*dt;

else
%Bessel coordinates evolution associated with first ten steps.
w(wnStart:Nn)= wMid+Zt2(wnStart:Nn)+wMu0dt(wnStart:Nn);
x(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
%x(wnStart:Nn)= +x2(wnStart:Nn)-sigma0.^2*gamma*x(wnStart:Nn).^(2*gamma-1)*dt+  ...
%    xMu0dt(wnStart:Nn);%+Zt2(wnStart:Nn)-Zt1(wnStart:Nn);%+ ...
end
%yy(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
w(wnStart:Nn)=x(wnStart:Nn).^(1-gamma)/(1-gamma);

%      [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));
%
%       w1(wnStart:Nn-1)=w(wnStart:Nn-1);
%       w1(Nn)=x(Nn);
%       w2(wnStart:Nn-1)=w(wnStart+1:Nn);
%       w2(Nn)=wE;
%       w(w1(:)>w2(:))=0;%Be careful;might not universally hold;
%       %  Change 3:7/25/2020: I have improved zero correction in above.
%       w(w<0)=0.0;
for nn=wnStart:Nn
if(x(nn)<=0)
wnStart=nn+1;
end
end

end
%yy(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
yy(wnStart:Nn)=x(wnStart:Nn);
Dfyy(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1

Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end

pyy(1:Nn)=0;
for nn = wnStart:Nn-1

pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
end

toc

ItoHermiteMean=sum(yy(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

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(yy(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-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(yy(wnStart+1:Nn-1),pyy(wnStart+1:Nn-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
.
.
.
Here is the output of the program
ItoHermiteMean =
0.249565254868068
true Mean only applicable to standard SV mean reverting type models otherwise disregard
TrueMean =
0.250251596970927
Original process average from monte carlo
MCMean =
0.249542499753951
Original process average from our simulation
ItoHermiteMean =
0.249565254868068
true Mean only applicble to standard SV mean reverting type models otherwise disregard
TrueMean =
0.250251596970927
IndexMax =
651

and here is the output graph(I have changed the scale)

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

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

There was an immediate relief in my back after after I wrote my previous two posts complaining about it and requesting investigative journalists to consider writing a story about mind control torture. I remained far better all through yesterday and even later today but when I slept last night, they continued to target my back with EM frequency weapons and there was a lot of pain in my back at several places when I woke up from sleep this morning.
Today when I went out for breakfast, I also bought some water bottles and diet coke and both water bottles and diet coke had mind control chemicals. They had already realized possible places I could go to get breakfast and had already contaminated water and other drinks before I reached there.
So yes, there was some let up on torture on my back but it would restart in full swing in a few days again when the whole thing goes out of limelight otherwise they continued to add mind control drugs to food and water just as usual.

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

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

I remained quite better all through yesterday except for a few things.
1. It was not otherwise an extremely serious problem but when I was doing new research yesterday, they continued to target my back with lasers.
2. I slept in air-conditioning in a different room(my room has terrible air-conditioning and gas that comes out of AC in my room is just thoroughly killing.) which has far better air-conditioning but they still use some light gases in air-conditioner and when I woke up, I was feeling frail and without any energy due to gas I had inhaled. I would have loved to live without air-conditioning but it is way too humid and warm at the moment.
3. When I went out in my car, they charged my head and face with gas that settles on exposed body. It makes me totally dull and takes all energy out of my body.

But still it was a far better day since they decreased laser on my back remarkably especially when I am just moving around and doing nothing very significant.

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

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

Friends, I have shifted further posts with my ongoing problems with mind control agencies to off-topic forum. Please continue to read there.
https://forum.wilmott.com/viewtopic.php?f=15&t=94796&p=867743#p867743
I hope friends would be sympathetic to me and try to get me off mind control. Thank you.

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

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

As I earlier mentioned, in a few days, I will be posting another worked out matlab program for solution of arbitrary stochastic volatility type system of stochastic differential equations. Please put the thread in your favorites as a lot more interesting stuff is on the way.

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

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

My post in off-topic where I candidly discuss mind control, religion, holocaust, Jews and anti-semitism. https://forum.wilmott.com/viewtopic.php?f=15&t=94796&p=867809#p867809

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

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

Friends, I am posting a monte carlo program for Stochastic volatility type SDEs. I have given two drift terms in both asset and SV  SDEs. You can specify your correlation rho as well. Here are the SDEs.

dX(t)=(a X(t)^alpha1 + b X(t)^alpha2) dt + sqrt(1-rho^2) sigma1 V(t)^gammaV X(t)^gammaX dZ1(t) + rho sigma1 V(t)^gammaV X(t)^gammaX dZ2(t)

dV(t)=(mu1 V(t)^beta1 + mu2 V(t)^beta2)*dt + sigma0 V(t)^gamma dZ2(t)

I will be posting analytic solution to the same PDE associated with these two SDEs framework in another two days. I have worked out the solution details of analytic solution.

Here is the monte carlo program. I have not made it computationally efficient and emphasis is on clarity of how I derived the solution so I have written every term as such. This is a quick post and I will make more detailed posts with 2D graphs of overlaid solution from monte carlo and analytic PDE solutions.

function [] = MonteCarlo2D()

% dX(t)=(a X(t)^alpha1 + b X(t)^alpha2) dt + sqrt(1-rho^2) sigma1 V(t)^gammaV X(t)^gammaX dZ1(t) + rho sigma1 V(t)^gammaV X(t)^gammaX dZ2(t)
% dV(t)=(mu1 V(t)^beta1 + mu2 V(t)^beta2)*dt + sigma0 V(t)^gamma dZ2(t)

%%%%%%%%%%%%%
X0=1.0;
alpha1=0;
alpha2=1;
kappaX=2.0;
thetaX=.5;
a=kappaX*thetaX;
b=-kappaX;
a=0;
b=0;
rho=-.7;
sigma1=.85;
gammaV=.5;
gammaX=.95;
%%%%%%%%%%%%%
V0=1.0;
thetaV=.25;
kappaV=4.0;
mu1=kappaV*thetaV;
mu2=-kappaV;
beta1=0;
beta2=1;
gamma=.95;
sigma0=1.25;
%%%%%%%%%%%%

dt=.03125;
Tt=32;

T=Tt*dt;
%paths=100000;

%%%%%%%%%%%%%%%%%%%%%%%%5

%rng(29079137, 'twister')
rng(15898837, 'twister')
paths=1000000;
V(1:paths)=V0;  %Original process monte carlo.
X(1:paths)=X0;
Random1(1:paths)=0;
Random2(1:paths)=0;

for tt=1:Tt
Random1=randn(size(Random1));
Random2=randn(size(Random2));

X(1:paths)=X(1:paths)+ ...
(a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dt + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) * sqrt(dt) + ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*sqrt(dt) + ...
(a*alpha1* X(1:paths).^(alpha1-1)+b*alpha2* X(1:paths).^(alpha2-1)).* ...
(((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dt^2/2)+ ...
(sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) *(1-1/sqrt(3)).*dt^1.5+ ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*(1-1/sqrt(3)).*dt^1.5))+ ...
.5*(a*alpha1*(alpha1-1)* X(1:paths).^(alpha1-2)+b*alpha2*(alpha2-1).* X(1:paths).^(alpha2-2)).* ...
( sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX)) *dt^2/2 + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random1(1:paths) * 1/sqrt(3).* dt^1.5 + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*(Random1(1:paths).^2-1) * dt/2 + ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths)*dt/2)+ ...
.5*sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
(sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random1(1:paths).*1/sqrt(3).*dt^1.5)+ ...
sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random1(1:paths) * 1/sqrt(3).* dt^1.5 + ...
sigma0*V(1:paths).^gamma.*Random1(1:paths).*Random2(1:paths)*dt/2)+ ...
.5*sqrt(1-rho^2)* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
(sigma0^2*V(1:paths).^(2*gamma).*Random1(1:paths)*1/sqrt(3)*dt^1.5)+ ...
sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random1(1:paths)*1/sqrt(3)*dt^1.5+ ...
rho* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random2(1:paths) * 1/sqrt(3).* dt^1.5 + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths) * dt/2 + ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*(Random2(1:paths).^2-1)*dt/2)+ ...
.5*rho* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
(sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random2(1:paths).*1/sqrt(3).*dt^1.5)+ ...
rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths) * 1/sqrt(3).* dt^1.5 + ...
sigma0*V(1:paths).^gamma.*(Random2(1:paths).^2-1)*dt/2)+ ...
.5*rho* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
sigma0^2.*V(1:paths).^(2*gamma).*Random2(1:paths) * 1/sqrt(3).* dt^1.5+ ...
rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random2(1:paths)*1/sqrt(3)*dt^1.5;

V(1:paths)=V(1:paths)+ ...
(mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dt + ...
sigma0*V(1:paths).^gamma .*Random2(1:paths)*sqrt(dt) + ...
(mu1.*beta1*V(1:paths).^(beta1-1) + mu2.*beta2.*V(1:paths).^(beta2-1)).* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dt^2/2 + ...
sigma0*V(1:paths).^gamma .*Random2(1:paths)*(1-1/sqrt(3))*dt^1.5) + ...
.5*(mu1.*beta1.*(beta1-1).*V(1:paths).^(beta1-2) + mu2.*beta2.*(beta2-1).*V(1:paths).^(beta2-2)).* ...
sigma0^2.*V(1:paths).^(2*gamma).*dt^2/2 + ...
sigma0*gamma*V(1:paths).^(gamma-1) .* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths).*1/sqrt(3)*dt^1.5 + ...
sigma0.*V(1:paths).^gamma .*(Random2(1:paths).^2-1)*dt/2) + ...
.5*sigma0*gamma*(gamma-1).*V(1:paths).^(gamma-2) .* ...
sigma0^2.*V(1:paths).^(2*gamma) .*Random2(1:paths).*1/sqrt(3)*dt^1.5;

end

thetaV+exp(-kappaV*T)*(V0-thetaV)
sum(V(1:paths))/paths
thetaX+exp(-kappaX*T)*(X0-thetaX)
sum(X(1:paths))/paths

end

Here is the output of this quick program.
>> MonteCarlo2D
ans =
0.263736729166551
ans =
0.263871276320159
ans =
0.567667641618306
ans =
0.999643417736976

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

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

Friends my above post has to be read in context of another post I made on the same topic exactly three years earlier. Here is the link to my previous post: Muslims, mind control, Muslim-Jewish relationship and American defense https://forum.wilmott.com/viewtopic.php?f=15&t=101453
When I wrote that previous post I tried to be as discrete as possible. Please note that my persecution started in year 1997 with mind control. So even when I was writing the above post after 21 years of mind control torture, persecution and insults, I never wanted to give any mild hints of anything anti-jewish. I was appealing to sympathies and humanity of good jews and other people thinking they will intervene to end my persecution. If I had any intentions to arouse resentment against jews, my wording of above post would have been extremely different. But my taking a mild stance only emboldened the evil jewish crooks and my persecution only increased and never decreased. I can take only enough abuse and torture for the rest of my life and therefore I had to write the post that I made two days ago. I am sure some people will call me anti-semite even after knowing the context of my persecution but I can tell you that if I am an anti-semite then any human being other than Christ will become an anti-semite when you will put him in my shoes.
I will request good and responsible people among Jews to take some concrete steps to stop the evil jewish crooks who have gone on a mind control rampage across the United States and rest of the world.
For the context of my persecution, please also read: My Letter To HRW About My Human Rights Abuse https://forum.wilmott.com/viewtopic.php?f=15&t=102298

In the thread below I tried in vain to request prime minister of my country to end mind control here but he is helpless before crook generals in Pakistan army who take hundreds of millions of dollars in bribes from CIA. https://forum.wilmott.com/viewtopic.php?f=15&t=94796&start=1125#p858119

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

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

Friends, when I posted the last program, it was very late at night, it was almost morning therefore I could not add any graphs for the monte carlo densities. I have added the graphs to monte carlo densities now and posting the program again. Again it is a program with exactly the same main body and some initial parameters are different and graphs are added. Here is the old program with graphs now
.
function [] = MonteCarlo2D()

% dX(t)=(a X(t)^alpha1 + b X(t)^alpha2) dt + sqrt(1-rho^2) sigma1 V(t)^gammaV X(t)^gammaX dZ1(t) + rho sigma1 V(t)^gammaV X(t)^gammaX dZ2(t)
% dV(t)=(mu1 V(t)^beta1 + mu2 V(t)^beta2)*dt + sigma0 V(t)^gamma dZ2(t)

%%%%%%%%%%%%%
X0=1.0;
alpha1=0;
alpha2=1;
%kappaX=2.0;
%thetaX=.5;
kappaX=0.0;
thetaX=0.0;
a=kappaX*thetaX;
b=-kappaX;
%a=0;
%b=0;
rho=-.99;
sigma1=.85;
gammaV=.5;
gammaX=.65;
%%%%%%%%%%%%%
V0=1.0;
thetaV=.25;
kappaV=4.0;
mu1=kappaV*thetaV;
mu2=-kappaV;
beta1=0;
beta2=1;
gamma=.95;
sigma0=1.25;
%%%%%%%%%%%%

dt=.03125;
Tt=32;

T=Tt*dt;
%paths=100000;

%%%%%%%%%%%%%%%%%%%%%%%%5

%rng(29079137, 'twister')
rng(15898837, 'twister')
paths=200000;
V(1:paths)=V0;  %Original process monte carlo.
X(1:paths)=X0;
Random1(1:paths)=0;
Random2(1:paths)=0;

for tt=1:Tt
Random1=randn(size(Random1));
Random2=randn(size(Random2));

X(1:paths)=X(1:paths)+ ...
(a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dt + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) * sqrt(dt) + ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*sqrt(dt) + ...
(a*alpha1* X(1:paths).^(alpha1-1)+b*alpha2* X(1:paths).^(alpha2-1)).* ...
(((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dt^2/2)+ ...
(sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) *(1-1/sqrt(3)).*dt^1.5+ ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*(1-1/sqrt(3)).*dt^1.5))+ ...
.5*(a*alpha1*(alpha1-1)* X(1:paths).^(alpha1-2)+b*alpha2*(alpha2-1).* X(1:paths).^(alpha2-2)).* ...
( sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX)) *dt^2/2 + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random1(1:paths) * 1/sqrt(3).* dt^1.5 + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*(Random1(1:paths).^2-1) * dt/2 + ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths)*dt/2)+ ...
.5*sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
(sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random1(1:paths).*1/sqrt(3).*dt^1.5)+ ...
sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random1(1:paths) * 1/sqrt(3).* dt^1.5 + ...
sigma0*V(1:paths).^gamma.*Random1(1:paths).*Random2(1:paths)*dt/2)+ ...
.5*sqrt(1-rho^2)* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
(sigma0^2*V(1:paths).^(2*gamma).*Random1(1:paths)*1/sqrt(3)*dt^1.5)+ ...
sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random1(1:paths)*1/sqrt(3)*dt^1.5+ ...
rho* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random2(1:paths) * 1/sqrt(3).* dt^1.5 + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths) * dt/2 + ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*(Random2(1:paths).^2-1)*dt/2)+ ...
.5*rho* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
(sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random2(1:paths).*1/sqrt(3).*dt^1.5)+ ...
rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths) * 1/sqrt(3).* dt^1.5 + ...
sigma0*V(1:paths).^gamma.*(Random2(1:paths).^2-1)*dt/2)+ ...
.5*rho* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
sigma0^2.*V(1:paths).^(2*gamma).*Random2(1:paths) * 1/sqrt(3).* dt^1.5+ ...
rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random2(1:paths)*1/sqrt(3)*dt^1.5;

V(1:paths)=V(1:paths)+ ...
(mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dt + ...
sigma0*V(1:paths).^gamma .*Random2(1:paths)*sqrt(dt) + ...
(mu1.*beta1*V(1:paths).^(beta1-1) + mu2.*beta2.*V(1:paths).^(beta2-1)).* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dt^2/2 + ...
sigma0*V(1:paths).^gamma .*Random2(1:paths)*(1-1/sqrt(3))*dt^1.5) + ...
.5*(mu1.*beta1.*(beta1-1).*V(1:paths).^(beta1-2) + mu2.*beta2.*(beta2-1).*V(1:paths).^(beta2-2)).* ...
sigma0^2.*V(1:paths).^(2*gamma).*dt^2/2 + ...
sigma0*gamma*V(1:paths).^(gamma-1) .* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths).*1/sqrt(3)*dt^1.5 + ...
sigma0.*V(1:paths).^gamma .*(Random2(1:paths).^2-1)*dt/2) + ...
.5*sigma0*gamma*(gamma-1).*V(1:paths).^(gamma-2) .* ...
sigma0^2.*V(1:paths).^(2*gamma) .*Random2(1:paths).*1/sqrt(3)*dt^1.5;

end

SVolMeanAnalytic=thetaV+exp(-kappaV*T)*(V0-thetaV)
SVolMean=sum(V(1:paths))/paths
AssetMeanAnalytic=thetaX+exp(-kappaX*T)*(X0-thetaX)
AssetMean=sum(X(1:paths))/paths

MaxCutOff=30;
%NoOfBins=round(600*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too
NoOfBins=2000;
[XDensity,IndexOutX,IndexMaxX] = MakeDensityFromSimulation_Infiniti_NEW(X,paths,NoOfBins,MaxCutOff );
plot(IndexOutX(1:IndexMaxX),XDensity(1:IndexMaxX),'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,thetaX=%.3f,kappaX=%.2f,gammaX=%.3f,gammaV=%.3f,sigma1=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', X0,thetaX,kappaX,gammaX,gammaV,sigma1,T,dt,AssetMean,AssetMeanAnalytic));%,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('Green line is density of Asset SDE.');

MaxCutOff=30;
%NoOfBins=round(600*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa));%Decrease the number of bins if the graph is too
NoOfBins=2000;
[VDensity,IndexOutV,IndexMaxV] = MakeDensityFromSimulation_Infiniti_NEW(V,paths,NoOfBins,MaxCutOff );
plot(IndexOutV(1:IndexMaxV),VDensity(1:IndexMaxX),'r');
%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('V0 = %.4f,thetaV=%.3f,kappaV=%.2f,gamma=%.3f,sigma0=%.2f,T=%.2f,dt=%.5f,SVolMean=%.4f,SVolMeanAnalytic=%.4f', V0,thetaV,kappaV,gamma,sigma0,T,dt,SVolMean,SVolMeanAnalytic));%,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 Asset SDE.');

end
.
.
Here is the final output of the program.
SVolMeanAnalytic =
0.263736729166551
SVolMean =
0.263748723062681
AssetMeanAnalytic =
1
AssetMean =
1.000894173714266 + 0.000003747246978i
IndexMax =
2001
Green line is density of Asset SDE.
IndexMax =
2001
red line is density of Asset SDE.

and following is the asset graph you will see.

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

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

Friends, today I received biweekly antipsychotic injection that is forced on me to keep me in control. I have begged so many times on this forum that I have been given antipsychotic drugs and injections for more than twenty years and now I am (48 years old) and it is too difficult  to continue to take anti-psychotic drugs and still do everything nicely. I am 30 kg over weight due to antipsychotic injections and mind control. At this stage in my life, it is very difficult for me to continue to take mind control and other abuse so I really request good Jews to intervene and force the crooks to end my mind control and mind control of thousands of other victims like me.
I have fought a very long battle against mind control agencies and I believe my true contribution to our societies (American and other)  is not in the field of mathematics, rather it is in exposing the actors behind mind control and exposing the cruel nature of mind control practices and I am sure my efforts will have quite a reasonable effect in stopping mind control of new victims who would otherwise be boldly(and surely) retarded if there were no awareness in American society about mind control. I want to request all good Americans to do whatever they can to end the menace of mind control from all human societies.
I had previously written about mind control by barack obama and requested the new president Biden to end mind control practices in United States and all over the world. Please read my posts in light of the information I have revealed in past few days.

Obama and Mind Control     https://forum.wilmott.com/viewtopic.php?f=15&p=864314#p864253

and finally appeal to new president Joe Biden to end mind control
https://forum.wilmott.com/viewtopic.php?f=15&p=864314#p864314

And please read here how my tweets to American senators are stopped by American defense and mind control crooks.
https://forum.wilmott.com/viewtopic.php?f=15&p=864314#p864276

and finally Blacks and mind control
Blacks and Mind Control

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

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

Sorry Double
Last edited by Amin on August 15th, 2021, 9:08 am, edited 1 time in total.