Serving the Quantitative Finance Community
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, I have been able to solve the SDEs of CEV noises and drifting SDEs (that take one term in their drift) very well by adding variances using a very simple algorithm that employs hermite polynomials. The method seems to work perfectly. I will make a detailed post explaining the logic of the method in a few hours. I will also post the matlab program. I will come back in a few hours.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Here are some graphs of solution of drifting SDEs over a large range of parameters. I have only matched variance by squared addition of coefficients of hermite polynomials. Please note that in some graphs density goes into zero but then comes back into positive territory again and shows up in the graph again. This can easily be obviated by removing points on grid (further into negative) after density has reached zero.
I still could not get mean-reverting SDEs to work with adding variances.
All parameters of drifting SDEs are on the title. True mean on title written as TMean is not valid since it was only known for mean reverting SDEs.
.
.
.
Here are the graphs.
I still could not get mean-reverting SDEs to work with adding variances.
All parameters of drifting SDEs are on the title. True mean on title written as TMean is not valid since it was only known for mean reverting SDEs.
.
.
.
Here are the graphs.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Here is the program used to make graphs above. I will write a thoroughly detailed explanation tomorrow. I did not make graphs of CEV noises since they were relatively less interesting as compared to drifting SDEs even though the fit to monte carlo density was extremely excellent for CEV noises. But you can run this program for all sort of parameters.
.
.
.
.
When you run this program, you will see the last graph in the previous post. I did rescale the graph.
.
.
.
Code: Select all
function [] = FPESeriesCoeffsBesselWithVariance04()
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=mu1 x(t)^beta1 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/2/2/2; % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=32*2*10;%16*2*4;%*4*4*1;%*4;%128*2*2*2; % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
T=Tt*dt;
Ts=Tt;
ds(1:Ts)=dt;
%
% %Below I have done calculations for smaller step size at start.
%
% %ds(1:64)=dt/16;
% %ds(65:128)=dt/4;
% if(Tt<=4)
% Ts=(Tt*8);%+(64-16);
% ds(1:Ts)=dt/8;
% end
% if(Tt>4)&&(Tt<=20)
% Ts=(32)+((Tt-4)*2);
% ds(1:32)=dt/8;
% ds(33:Ts)=dt/2;
% end
% if(Tt>20)
% Ts=(32)+((20-4)*2)+(Tt-20);
% ds(1:32)=dt/8;
% ds(33:64)=dt/2;
% ds(65:Ts)=dt;
% end
% T
% sum(ds(1:Ts))
% Ts
%str=input('Check if time is right');
OrderA=4; %
OrderM=4; %
if(T>=1)
dtM=1/16;%.0625/1;
TtM=T/dtM;
else
dtM=T/16;%.0625/1;
TtM=T/dtM;
end
dtM=dt*8;
TtM=Tt/8;
dNn=.2/2; % Normal density subdivisions width. would change with number of subdivisions
Nn=47*2; % No of normal density subdivisions
NnMid=4.0;
x0=2.0; % starting value of SDE
beta1=0.65;
beta2=1.0; % Second drift term power.
%beta1=0.950;
%beta2=0.0; % Second drift term power.
gamma=.95;%50; % volatility power.
kappa=1.00;%.950; %mean reversion parameter.
theta=1.0;%mean reversion target
sigma0=.250;%Volatility value
%you can specify any general mu1 and mu2 and beta1 and beta2.
mu1=.10; %first drift coefficient.
mu2=0;%-1*kappa; % Second drift coefficient.
%mu1=theta*kappa;
%mu2=-1*kappa;
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;
w0=x0^(1-gamma)/(1-gamma);
Z(1:Nn)=(((1:Nn)-7.5/1)*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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ExpnOrder=5;
SeriesOrder=5;
wnStart=1;
tic
for tt=1:Ts
tt
if(tt==1)
[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
c0=wMu0dt+w0;
c(1)=sigma0*sqrt(ds(tt));
c(2:10)=0.0;
w0=c0;
elseif(tt>=2)
[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
%Below we get Z-series representation of drift term.
[b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
%Below we convert Z_series of drift term to hermite representation
[bH0,bH] = ConvertZCoeffsToHCoeffs(b0,b,SeriesOrder);
[wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
%Below we get Z-series representation of coefficient of first hermite polynomial in diffusion term.
[g10,g1] = CalculateDriftbCoeffs08A(wVol0dt,dwVol0dtdw,c,SeriesOrder);
g10=g10+sigma0*sqrt(dt);
%The Z_series(called g1) below is coefficient of first hermite polynomial
%g10+g1(1)*Z+g1(2)*Z^2+g1(3)*Z^3+g1(4)*Z^4+g1(5)*Z^5
[g1H0,g1H] = ConvertZCoeffsToHCoeffs(g10,g1,SeriesOrder);
%g1Sqr is variance of Z_series g1 that makes coefficient of first
%hermite polynomial
g1Sqr=(sign(g1H0).*g1H0^2+sign(g1H(1)).*g1H(1).^2+sign(g1H(2)).*g1H(2).^2*2+ ...
sign(g1H(3)).*g1H(3).^2*6+sign(g1H(4)).*g1H(4).^2*24+sign(g1H(5)).*g1H(5).^2*120);
g11=sign(g1Sqr).*sqrt(abs(g1Sqr));
[wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
%Below we get Z-series representation of coefficient of second hermite polynomial in diffusion term.
[g20,g2] = CalculateDriftbCoeffs08A(wVol2dt,dwVol2dtdw,c,SeriesOrder);
%The Z_series(called g2) below is coefficient of second hermite polynomial
%g20+g2(1)*Z+g2(2)*Z^2+g2(3)*Z^3+g2(4)*Z^4+g2(5)*Z^5
[g2H0,g2H] = ConvertZCoeffsToHCoeffs(g20,g2,SeriesOrder);
%g2Sqr is variance of Z_series g2 that makes coefficient of second
%hermite polynomial
g2Sqr=(sign(g2H0).*g2H0^2+sign(g2H(1)).*g2H(1).^2+sign(g2H(2)).*g2H(2).^2*2+ ...
sign(g2H(3)).*g2H(3).^2*6+sign(g2H(4)).*g2H(4).^2*24+sign(g2H(5)).*g2H(5).^2*120);
g22=sign(g2Sqr).*sqrt(abs(g2Sqr));
%Below are derived hermite polynomial representation of SDE variable.
[cH0,cH] = ConvertZCoeffsToHCoeffs(c0,c,SeriesOrder);
%Drift hermite coefficients are added to SDE variable hermite
%coefficients variables linearly. This is done before adding variances
cH0=cH0+bH0;
cH=cH+bH;
%Since our diffusion term has only first and second hermites, we add
%them in a squared fashion to hermite coefficients of SDE variable.
%Other than first and second hermite, coefficients of other hermite
%polynomials representing the SDE variable remain the same.
wH0=cH0;
wH(1)=sign(cH(1)+g11).*sqrt(abs(sign(cH(1)).*cH(1).^2+g1Sqr));
wH(2)=sign(cH(2)+g22).*sqrt(abs(sign(cH(2)).*cH(2).^2+g2Sqr));
wH(3:SeriesOrder)=cH(3:SeriesOrder);
[c0,c] = ConvertHCoeffsToZCoeffs(wH0,wH,SeriesOrder);%
%Above we convert from hermite coefficients representation of SDE to
%simple Z-series representation.
w0=c0;
end
end
wnStart=1;
%yy0=((1-gamma).*w0).^(1/(1-gamma));
w(1:Nn)=c0;
for nn=1:ExpnOrder
w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^nn;
end
%str=input('Look at numbers')
yy0(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
yy=yy0;
Dfyy(wnStart:Nn)=0;
%Dfyydt(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
% Dfyydt(nn) = (yydt(nn + 1) - yydt(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
%Dfyy(Nn)=Dfyy(Nn-1);
%Dfyy(1)=Dfyy(2);
%Dfyydt(Nn)=Dfyydt(Nn-1);
%Dfyydt(1)=Dfyydt(2);
pyy(1:Nn)=0;
%pyydt(1:Nn)=0;
for nn = wnStart+1:Nn-1
pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
% pyydt(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyydt(nn));
end
toc
ItoHermiteMean=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %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*T)%Mean reverting SDE original variable true average
%yy0
theta1=1;
rng(29079137, 'twister')
paths=50000;
YY(1:paths)=x0; %Original process monte carlo.
%YYdt(1:paths)=0;
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;
%YYdt(1:paths)=YYdt(1:paths)+.5*YY(1:paths)*dtM;
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);
% YYdt(1:paths)=YYdt(1:paths)+.5*YY(1:paths)*dtM;
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:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
%ItoHermiteMeanPI=sum(yydt(wnStart:Nn).*ZProb(wnStart:Nn))
%MCMeanPI=sum(YYdt(:))/paths
disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average
MaxCutOff=50;
NoOfBins=round(8*500*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));
title(sprintf('x0 = %.4f,mu1=%.3f,beta1=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,mu1,beta1,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.');
%NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa))/2;
%[YdtDensity,IndexOutYdt,IndexMaxYdt] = MakeDensityFromSimulation_Infiniti_NEW(YYdt,paths,NoOfBins,MaxCutOff );
%plot(yydt(wnStart+1:Nn-1),pyydt(wnStart+1:Nn-1),'r',IndexOutYdt(1:IndexMaxYdt),YdtDensity(1:IndexMaxYdt),'g');
%title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
%title(sprintf('x0 = %.4f,kappa=%.3f,theta=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,kappa,theta,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
%legend({'Ito-Hermite Path Integral Density','Monte Carlo Path Integral Density'},'Location','northeast')
%str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
endWhen you run this program, you will see the last graph in the previous post. I did rescale the graph.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, please comment line 59 and line 60 of the previous program. It makes monte carlo step too large and sometimes problem might be with monte carlo (and not in analytic program) if densities do not match since monte carlo step is too large.
Here I write the lines(59 and 60) after commenting them.
%dtM=dt*8;
%TtM=Tt/8;
Or rather you can use this program where I have commented the lines.
Here I write the lines(59 and 60) after commenting them.
%dtM=dt*8;
%TtM=Tt/8;
Or rather you can use this program where I have commented the lines.
Code: Select all
function [] = FPESeriesCoeffsBesselWithVariance04()
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334
%or skype ahsan.amin2999
%In this program, I am simulating the SDE given as
%dy(t)=mu1 x(t)^beta1 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/2/2/2; % Simulation time interval.%Fodiffusions close to zero
%decrease dt for accuracy.
Tt=32*2*10;%16*2*4;%*4*4*1;%*4;%128*2*2*2; % Number of simulation levels. Terminal time= Tt*dt; //.125/32*32*16=2 year;
T=Tt*dt;
Ts=Tt;
ds(1:Ts)=dt;
%
% %Below I have done calculations for smaller step size at start.
%
% %ds(1:64)=dt/16;
% %ds(65:128)=dt/4;
% if(Tt<=4)
% Ts=(Tt*8);%+(64-16);
% ds(1:Ts)=dt/8;
% end
% if(Tt>4)&&(Tt<=20)
% Ts=(32)+((Tt-4)*2);
% ds(1:32)=dt/8;
% ds(33:Ts)=dt/2;
% end
% if(Tt>20)
% Ts=(32)+((20-4)*2)+(Tt-20);
% ds(1:32)=dt/8;
% ds(33:64)=dt/2;
% ds(65:Ts)=dt;
% end
% T
% sum(ds(1:Ts))
% Ts
%str=input('Check if time is right');
OrderA=4; %
OrderM=4; %
if(T>=1)
dtM=1/16;%.0625/1;
TtM=T/dtM;
else
dtM=T/16;%.0625/1;
TtM=T/dtM;
end
%dtM=dt*8;
%TtM=Tt/8;
dNn=.2/2; % Normal density subdivisions width. would change with number of subdivisions
Nn=47*2; % No of normal density subdivisions
NnMid=4.0;
x0=2.0; % starting value of SDE
beta1=0.65;
beta2=1.0; % Second drift term power.
%beta1=0.950;
%beta2=0.0; % Second drift term power.
gamma=.95;%50; % volatility power.
kappa=1.00;%.950; %mean reversion parameter.
theta=1.0;%mean reversion target
sigma0=.250;%Volatility value
%you can specify any general mu1 and mu2 and beta1 and beta2.
mu1=.10; %first drift coefficient.
mu2=0;%-1*kappa; % Second drift coefficient.
%mu1=theta*kappa;
%mu2=-1*kappa;
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;
w0=x0^(1-gamma)/(1-gamma);
Z(1:Nn)=(((1:Nn)-7.5/1)*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);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ExpnOrder=5;
SeriesOrder=5;
wnStart=1;
tic
for tt=1:Ts
tt
if(tt==1)
[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
c0=wMu0dt+w0;
c(1)=sigma0*sqrt(ds(tt));
c(2:10)=0.0;
w0=c0;
elseif(tt>=2)
[wMu0dt,dwMu0dtdw] = BesselDriftAndDerivativesH0(w0,YqCoeff0,Fp1,gamma,ds(tt),6);
%Below we get Z-series representation of drift term.
[b0,b] = CalculateDriftbCoeffs08A(wMu0dt,dwMu0dtdw,c,SeriesOrder);
%Below we convert Z_series of drift term to hermite representation
[bH0,bH] = ConvertZCoeffsToHCoeffs(b0,b,SeriesOrder);
[wVol0dt,dwVol0dtdw] = BesselVolAndDerivativesH1(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
%Below we get Z-series representation of coefficient of first hermite polynomial in diffusion term.
[g10,g1] = CalculateDriftbCoeffs08A(wVol0dt,dwVol0dtdw,c,SeriesOrder);
g10=g10+sigma0*sqrt(dt);
%The Z_series(called g1) below is coefficient of first hermite polynomial
%g10+g1(1)*Z+g1(2)*Z^2+g1(3)*Z^3+g1(4)*Z^4+g1(5)*Z^5
[g1H0,g1H] = ConvertZCoeffsToHCoeffs(g10,g1,SeriesOrder);
%g1Sqr is variance of Z_series g1 that makes coefficient of first
%hermite polynomial
g1Sqr=(sign(g1H0).*g1H0^2+sign(g1H(1)).*g1H(1).^2+sign(g1H(2)).*g1H(2).^2*2+ ...
sign(g1H(3)).*g1H(3).^2*6+sign(g1H(4)).*g1H(4).^2*24+sign(g1H(5)).*g1H(5).^2*120);
g11=sign(g1Sqr).*sqrt(abs(g1Sqr));
[wVol2dt,dwVol2dtdw] = BesselVolAndDerivativesH2(w0,YqCoeff0,Fp1,gamma,ds(tt),SeriesOrder+1);
%Below we get Z-series representation of coefficient of second hermite polynomial in diffusion term.
[g20,g2] = CalculateDriftbCoeffs08A(wVol2dt,dwVol2dtdw,c,SeriesOrder);
%The Z_series(called g2) below is coefficient of second hermite polynomial
%g20+g2(1)*Z+g2(2)*Z^2+g2(3)*Z^3+g2(4)*Z^4+g2(5)*Z^5
[g2H0,g2H] = ConvertZCoeffsToHCoeffs(g20,g2,SeriesOrder);
%g2Sqr is variance of Z_series g2 that makes coefficient of second
%hermite polynomial
g2Sqr=(sign(g2H0).*g2H0^2+sign(g2H(1)).*g2H(1).^2+sign(g2H(2)).*g2H(2).^2*2+ ...
sign(g2H(3)).*g2H(3).^2*6+sign(g2H(4)).*g2H(4).^2*24+sign(g2H(5)).*g2H(5).^2*120);
g22=sign(g2Sqr).*sqrt(abs(g2Sqr));
%Below are derived hermite polynomial representation of SDE variable.
[cH0,cH] = ConvertZCoeffsToHCoeffs(c0,c,SeriesOrder);
%Drift hermite coefficients are added to SDE variable hermite
%coefficients variables linearly. This is done before adding variances
cH0=cH0+bH0;
cH=cH+bH;
%Since our diffusion term has only first and second hermites, we add
%them in a squared fashion to hermite coefficients of SDE variable.
%Other than first and second hermite, coefficients of other hermite
%polynomials representing the SDE variable remain the same.
wH0=cH0;
wH(1)=sign(cH(1)+g11).*sqrt(abs(sign(cH(1)).*cH(1).^2+g1Sqr));
wH(2)=sign(cH(2)+g22).*sqrt(abs(sign(cH(2)).*cH(2).^2+g2Sqr));
wH(3:SeriesOrder)=cH(3:SeriesOrder);
[c0,c] = ConvertHCoeffsToZCoeffs(wH0,wH,SeriesOrder);%
%Above we convert from hermite coefficients representation of SDE to
%simple Z-series representation.
w0=c0;
end
end
wnStart=1;
%yy0=((1-gamma).*w0).^(1/(1-gamma));
w(1:Nn)=c0;
for nn=1:ExpnOrder
w(1:Nn)=w(1:Nn)+c(nn)*Z(1:Nn).^nn;
end
%str=input('Look at numbers')
yy0(wnStart:Nn)=((1-gamma).*w(wnStart:Nn)).^(1/(1-gamma));
yy=yy0;
Dfyy(wnStart:Nn)=0;
%Dfyydt(wnStart:Nn)=0;
for nn=wnStart+1:Nn-1
Dfyy(nn) = (yy(nn + 1) - yy(nn - 1))/(Z(nn + 1) - Z(nn - 1));
% Dfyydt(nn) = (yydt(nn + 1) - yydt(nn - 1))/(Z(nn + 1) - Z(nn - 1));
%Change of variable derivative for densities
end
%Dfyy(Nn)=Dfyy(Nn-1);
%Dfyy(1)=Dfyy(2);
%Dfyydt(Nn)=Dfyydt(Nn-1);
%Dfyydt(1)=Dfyydt(2);
pyy(1:Nn)=0;
%pyydt(1:Nn)=0;
for nn = wnStart+1:Nn-1
pyy(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyy(nn));
% pyydt(nn) = (normpdf(Z(nn),0, 1))/abs(Dfyydt(nn));
end
toc
ItoHermiteMean=sum(yy(wnStart:Nn).*ZProb(wnStart:Nn)) %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*T)%Mean reverting SDE original variable true average
%yy0
theta1=1;
rng(29079137, 'twister')
paths=50000;
YY(1:paths)=x0; %Original process monte carlo.
%YYdt(1:paths)=0;
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;
%YYdt(1:paths)=YYdt(1:paths)+.5*YY(1:paths)*dtM;
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);
% YYdt(1:paths)=YYdt(1:paths)+.5*YY(1:paths)*dtM;
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:Nn).*ZProb(wnStart:Nn)) %Original process average from coordinates
%ItoHermiteMeanPI=sum(yydt(wnStart:Nn).*ZProb(wnStart:Nn))
%MCMeanPI=sum(YYdt(:))/paths
disp('true Mean only applicble to standard SV mean reverting type models otherwise disregard');
TrueMean=theta+(x0-theta)*exp(-kappa*T)%Mean reverting SDE original variable true average
MaxCutOff=50;
NoOfBins=round(8*500*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));
title(sprintf('x0 = %.4f,mu1=%.3f,beta1=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,TM=%.4f', x0,mu1,beta1,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.');
%NoOfBins=round(1*500*gamma^2*4*sigma0/sqrt(MCMean)/(1+kappa))/2;
%[YdtDensity,IndexOutYdt,IndexMaxYdt] = MakeDensityFromSimulation_Infiniti_NEW(YYdt,paths,NoOfBins,MaxCutOff );
%plot(yydt(wnStart+1:Nn-1),pyydt(wnStart+1:Nn-1),'r',IndexOutYdt(1:IndexMaxYdt),YdtDensity(1:IndexMaxYdt),'g');
%title(sprintf('x0 = %.4f,theta=%.3f,kappa=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,theta,kappa,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
%title(sprintf('x0 = %.4f,kappa=%.3f,theta=%.2f,gamma=%.3f,sigma=%.2f,T=%.2f,dt=%.5f,M=%.4f,MCM=%.4f', x0,kappa,theta,gamma,sigma0,T,dt,ItoHermiteMeanPI,MCMeanPI));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
%legend({'Ito-Hermite Path Integral Density','Monte Carlo Path Integral Density'},'Location','northeast')
%str=input('red line is density of SDE from Ito-Hermite method, green is monte carlo.');
endYou think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, I had an anti-psychotic injection on Wednesday. 2-7 days after the injection are very difficult for me. I was very lethargic and drowsy all day. Then I slept at 2:00 pm and just woke up now at 9:30 pm. Sorry that I could not write the explanation for the program I posted.
I am not sure but I think I might know why my algorithm is not working for mean-reverting SDEs. I will try fixing it and posting an updated program with explanation later tomorrow.
Though I had an anti-psychotic injection just on wednesday, mind control agencies still continued to try to drug my food and water. I am tired of complaining about their behavior and I want to tell friends that when I do not complain about them it does not mean that they have stopped my persecution activity. Today Friday morning when I woke up, I went out for a 30 km drive around the city for the fun of it and I still got drugged water at many remote places in the city. Early in the morning 99% of the shops were closed and they drugged water at 1% of the shops that were open.
And this pattern of drugging my water with mind control chemicals has continued for all of the past week.
I am not sure but I think I might know why my algorithm is not working for mean-reverting SDEs. I will try fixing it and posting an updated program with explanation later tomorrow.
Though I had an anti-psychotic injection just on wednesday, mind control agencies still continued to try to drug my food and water. I am tired of complaining about their behavior and I want to tell friends that when I do not complain about them it does not mean that they have stopped my persecution activity. Today Friday morning when I woke up, I went out for a 30 km drive around the city for the fun of it and I still got drugged water at many remote places in the city. Early in the morning 99% of the shops were closed and they drugged water at 1% of the shops that were open.
And this pattern of drugging my water with mind control chemicals has continued for all of the past week.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
.
Here I write the method I have used to analytically simulate the noise and drifting SDEs(that take one drift term) by adding the variances. The program I posted two days ago is based on this method.
Let us suppose our Bessel SDE has two drift terms and we write the Bessel transformed SDE equation as
[$]\, dw(t) \, = \, \mu_1 \, {w(t)}^{\zeta_1} \, dt \,+ \, \mu_2 \, {w(t)}^{\zeta_2} \, dt \,+ \sigma \,dz(t) \,[$]
Please note that [$]\mu_1[$], [$]\mu_2[$], [$]\zeta_1[$], and [$]\zeta_2[$] all take specific values in Bessel transformed SDE based on values of coefficients of the SDE in original coordinates.
At every analytic simulation step, we are initially given a Z-series representation of the SDE variable [$] w(t) [$] , and this value is represented (to fifth order) as
[$] \, w(t) \, = \, c_0 \, + \, c_1 Z \, + \, c_2 \, Z^2 \, + \, c_3 \, Z^3 \, + \, c_4 \, Z^4 \, + \, c_5 \, Z^5 \, [$]
Furthermore, we have an algebraically equivalent series in terms of Hermite polynomials represented as
[$] \, w(t) \, = \,ch_0 \, + \, ch_1 H_1 \, + \, ch_2 \, H_2 \, + \, ch_3 \, H_3 \, + \, ch_4 \, H_4 \, + \, ch_5 \, H_5 \, [$]
First we find a higher order expansion of the Bessel transformed SDE. We have covered these details several times while we were working on higher order expansions of simple monte carlo simulations earlier in this thread. Suppose our expansion for drift term is of the form
[$]= [\, f_1(w)\, +\, f_2(w)\, ] \,dt + [\, f_3(w)\, +\, f_4(w)\,+\, f_5(w)\, ] \,dt^2[$]
We calculate the value of above drift term expansion and its derivatives at the median of density of the SDE and using the derivatives represent it as a Z_series as we have learnt previously in this thread several times. We represent this Z-series expansion of the drift terms as
[$][\, f_1(w)\, +\, f_2(w)\, ] \,dt + [\, f_3(w)\, +\, f_4(w)\,+\, f_5(w)\, ] \,dt^2[$]
[$] \, = \, b_0 \, + \, b_1 Z \, + \, b_2 \, Z^2 \, + \, b_3 \, Z^3 \, + \, b_4 \, Z^4 \, + \, b_5 \, Z^5 \, [$]
[$] \, = \, bh_0 \, + \, bh_1 H_1 \, + \, bh_2 \, H_2 \, + \, bh_3 \, H_3 \, + \, bh_4 \, H_4 \, + \, bh_5 \, H_5 \, [$]
Similar to higher order expansion of drift, we find high order expansions for volatility terms (As we have been doing for monte carlo and other expansions for past several years). Let us suppose this expansion is given as
[$]\sigma \sqrt{dt} \, Z_1+ \sigma \, dt^{1.5} \Big[ c_1 \, w(t)^{\alpha_1} \, + \, c_2 \, w(t)^{\alpha_2} \, \, \Big] \, Z_1[$]
[$]+\sigma^2 \, dt^{2.0} \Big[ d_1 \, w(t)^{\beta_1} \, + \, d_2 \, w(t)^{\beta_2} \, \Big] \, ({Z_1}^2 \, - 1)\, [$]
Please note that in above expansion the transition normal [$]Z_1[$] is a normal variable independent and orthogonal to previous instance of normal random number denoted as Z.
Just like we did for volatility, we will expand the coefficient terms of both transition hermites [$]Z_1[$] and [$]{Z_1}^2 \, - 1[$] in terms of Z-series. The coefficient terms of first transition hermite as expanded as
[$]\sigma \sqrt{dt} \, \,+ \sigma \, dt^{1.5} \Big[ c_1 \, w(t)^{\alpha_1} \, + \, c_2 \, w(t)^{\alpha_2} \, \, \Big][$]
[$]=g_{10} \, + \, g_{11} Z + \, g_{12} {Z}^2 \, + \, g_{13} {Z}^3 \, + \, g_{14} {Z}^4 \,+ \, g_{15} {Z}^5[$]
Which has an equivalent hermite representation as
[$]=gh_{10} \, + \, gh_{11} H_1(Z) + \, gh_{12} H_2(Z) \, + \, gh_{13} H_3(Z) \, + \, gh_{14} H_4(Z) \,+ \, gh_{15} H_5(Z)[$]
This makes our expression for first transition hermite equivalent to
[$]\sigma \sqrt{dt} \, Z_1+ \sigma \, dt^{1.5} \Big[ c_1 \, w(t)^{\alpha_1} \, + \, c_2 \, w(t)^{\alpha_2} \, \, \Big] \, Z_1[$]
[$]=\, \Big[\, gh_{10} \, + \, gh_{11} H_1(Z) + \, gh_{12} H_2(Z) \, + \, gh_{13} H_3(Z) \, + \, gh_{14} H_4(Z) \,+ \, gh_{15} H_5(Z) \, \Big]\, H_1(Z_1) [$]
In order to find one single coefficient of first transition hermite, we calculate square root of variance of its coefficient series as
[$]G_{11}=\sqrt{\, \Big[\, {gh_{10}}^2 \, + \, {gh_{11}}^2 +2 \, {gh_{12}}^2 \, + 6 \, {gh_{13}}^2 \, +24 \, {gh_{14}}^2 \,+120 \, {gh_{15}}^2 \, \Big]}\,[$]
We repeat the above step to find the square root of variance of second transition hermite. First we expand the coefficient terms of the second transition hermite in the form of a Z-series as
[$]\sigma^2 \, dt^{2.0} \Big[ d_1 \, w(t)^{\beta_1} \, + \, d_2 \, w(t)^{\beta_2} \, \Big] [$]
[$]=g_{20} \, + \, g_{21} Z + \, g_{22} {Z}^2 \, + \, g_{23} {Z}^3 \, + \, g_{24} {Z}^4 \,+ \, g_{25} {Z}^5[$]
Which has an equivalent hermite representation as
[$]=gh_{20} \, + \, gh_{21} H_1(Z) + \, gh_{22} H_2(Z) \, + \, gh_{23} H_3(Z) \, + \, gh_{24} H_4(Z) \,+ \, gh_{25} H_5(Z)[$]
This makes our expression for second transition hermite equivalent to
[$]\sigma^2 \, dt^{2.0} \Big[ d_1 \, w(t)^{\beta_1} \, + \, d_2 \, w(t)^{\beta_2} \, \Big] \, ({Z_1}^2 \, - 1)\, [$]
[$]=\, \Big[\, gh_{20} \, + \, gh_{21} H_1(Z) + \, gh_{22} H_2(Z) \, + \, gh_{23} H_3(Z) \, + \, gh_{24} H_4(Z) \,+ \, gh_{25} H_5(Z) \, \Big]\, H_2(Z_1) [$]
In order to find one coefficient of second transition hermite, we calculate square root of variance of its hermite coefficienmtseries as
[$]G_{22}=\sqrt{\, \Big[\, {gh_{20}}^2 \, + \, {gh_{21}}^2 +2 \, {gh_{22}}^2 \, + 6 \, {gh_{23}}^2 \, +24 \, {gh_{24}}^2 \,+120 \, {gh_{25}}^2 \, \Big]}\,[$]
Now we can write the simple equation to advance the hermite representation of the SDE in two steps. First we add drift hermite series to SDE hermite series linearly as
[$]ch_0 \,=\,ch_0 +bh_0[$]
and similarly for rest of the higher hrmites are updated as (n represents coefficient of nth hermite)
[$]ch_n\,=\,ch_n +bh_n[$]
And then we add first and second transition hermite to the SDE hermite representation in a squared fashion.
[$]ch_1\,=\sqrt{\, {ch_1}^2 \,+ \, {G_{11}}^2 \,}[$]
and for adding the variance of second hermite, it is done as
[$]ch_2\,=\sqrt{\, {ch_2}^2 \,+ \, {G_{22}}^2 \,}[$]
Transition variance is not added to SDE hermites of order greater than two since 3rd transition hermite in Bessel coordinates takes coefficients in order dt^2.5 which become pretty much negligible for the kind of time step we are takin in our analytic simulation.
Now we can reconstruct the hermite representation of w(t) at next time step with updated hermite coefficients as
[$] \, w(t+1) \, = \,ch_0 \, + \, ch_1 H_1 \, + \, ch_2 \, H_2 \, + \, ch_3 \, H_3 \, + \, ch_4 \, H_4 \, + \, ch_5 \, H_5 \, [$]
Please note that we have to continue to switch our series representation from Z-series to hermite series and back and forth. Z-series is found around median of the density and can be used to find a series representation of any function of SDE variable. It is Z-series that we use to calculate the coefficients of drift term expansion and expansion of coefficient terms of first and second transition hermite. We cannot do a change of variable with hermite representation(which is around the mean) that is great for addition of variances something that we cannot do with a Z-series representation.
Here I write the method I have used to analytically simulate the noise and drifting SDEs(that take one drift term) by adding the variances. The program I posted two days ago is based on this method.
Let us suppose our Bessel SDE has two drift terms and we write the Bessel transformed SDE equation as
[$]\, dw(t) \, = \, \mu_1 \, {w(t)}^{\zeta_1} \, dt \,+ \, \mu_2 \, {w(t)}^{\zeta_2} \, dt \,+ \sigma \,dz(t) \,[$]
Please note that [$]\mu_1[$], [$]\mu_2[$], [$]\zeta_1[$], and [$]\zeta_2[$] all take specific values in Bessel transformed SDE based on values of coefficients of the SDE in original coordinates.
At every analytic simulation step, we are initially given a Z-series representation of the SDE variable [$] w(t) [$] , and this value is represented (to fifth order) as
[$] \, w(t) \, = \, c_0 \, + \, c_1 Z \, + \, c_2 \, Z^2 \, + \, c_3 \, Z^3 \, + \, c_4 \, Z^4 \, + \, c_5 \, Z^5 \, [$]
Furthermore, we have an algebraically equivalent series in terms of Hermite polynomials represented as
[$] \, w(t) \, = \,ch_0 \, + \, ch_1 H_1 \, + \, ch_2 \, H_2 \, + \, ch_3 \, H_3 \, + \, ch_4 \, H_4 \, + \, ch_5 \, H_5 \, [$]
First we find a higher order expansion of the Bessel transformed SDE. We have covered these details several times while we were working on higher order expansions of simple monte carlo simulations earlier in this thread. Suppose our expansion for drift term is of the form
[$]= [\, f_1(w)\, +\, f_2(w)\, ] \,dt + [\, f_3(w)\, +\, f_4(w)\,+\, f_5(w)\, ] \,dt^2[$]
We calculate the value of above drift term expansion and its derivatives at the median of density of the SDE and using the derivatives represent it as a Z_series as we have learnt previously in this thread several times. We represent this Z-series expansion of the drift terms as
[$][\, f_1(w)\, +\, f_2(w)\, ] \,dt + [\, f_3(w)\, +\, f_4(w)\,+\, f_5(w)\, ] \,dt^2[$]
[$] \, = \, b_0 \, + \, b_1 Z \, + \, b_2 \, Z^2 \, + \, b_3 \, Z^3 \, + \, b_4 \, Z^4 \, + \, b_5 \, Z^5 \, [$]
[$] \, = \, bh_0 \, + \, bh_1 H_1 \, + \, bh_2 \, H_2 \, + \, bh_3 \, H_3 \, + \, bh_4 \, H_4 \, + \, bh_5 \, H_5 \, [$]
Similar to higher order expansion of drift, we find high order expansions for volatility terms (As we have been doing for monte carlo and other expansions for past several years). Let us suppose this expansion is given as
[$]\sigma \sqrt{dt} \, Z_1+ \sigma \, dt^{1.5} \Big[ c_1 \, w(t)^{\alpha_1} \, + \, c_2 \, w(t)^{\alpha_2} \, \, \Big] \, Z_1[$]
[$]+\sigma^2 \, dt^{2.0} \Big[ d_1 \, w(t)^{\beta_1} \, + \, d_2 \, w(t)^{\beta_2} \, \Big] \, ({Z_1}^2 \, - 1)\, [$]
Please note that in above expansion the transition normal [$]Z_1[$] is a normal variable independent and orthogonal to previous instance of normal random number denoted as Z.
Just like we did for volatility, we will expand the coefficient terms of both transition hermites [$]Z_1[$] and [$]{Z_1}^2 \, - 1[$] in terms of Z-series. The coefficient terms of first transition hermite as expanded as
[$]\sigma \sqrt{dt} \, \,+ \sigma \, dt^{1.5} \Big[ c_1 \, w(t)^{\alpha_1} \, + \, c_2 \, w(t)^{\alpha_2} \, \, \Big][$]
[$]=g_{10} \, + \, g_{11} Z + \, g_{12} {Z}^2 \, + \, g_{13} {Z}^3 \, + \, g_{14} {Z}^4 \,+ \, g_{15} {Z}^5[$]
Which has an equivalent hermite representation as
[$]=gh_{10} \, + \, gh_{11} H_1(Z) + \, gh_{12} H_2(Z) \, + \, gh_{13} H_3(Z) \, + \, gh_{14} H_4(Z) \,+ \, gh_{15} H_5(Z)[$]
This makes our expression for first transition hermite equivalent to
[$]\sigma \sqrt{dt} \, Z_1+ \sigma \, dt^{1.5} \Big[ c_1 \, w(t)^{\alpha_1} \, + \, c_2 \, w(t)^{\alpha_2} \, \, \Big] \, Z_1[$]
[$]=\, \Big[\, gh_{10} \, + \, gh_{11} H_1(Z) + \, gh_{12} H_2(Z) \, + \, gh_{13} H_3(Z) \, + \, gh_{14} H_4(Z) \,+ \, gh_{15} H_5(Z) \, \Big]\, H_1(Z_1) [$]
In order to find one single coefficient of first transition hermite, we calculate square root of variance of its coefficient series as
[$]G_{11}=\sqrt{\, \Big[\, {gh_{10}}^2 \, + \, {gh_{11}}^2 +2 \, {gh_{12}}^2 \, + 6 \, {gh_{13}}^2 \, +24 \, {gh_{14}}^2 \,+120 \, {gh_{15}}^2 \, \Big]}\,[$]
We repeat the above step to find the square root of variance of second transition hermite. First we expand the coefficient terms of the second transition hermite in the form of a Z-series as
[$]\sigma^2 \, dt^{2.0} \Big[ d_1 \, w(t)^{\beta_1} \, + \, d_2 \, w(t)^{\beta_2} \, \Big] [$]
[$]=g_{20} \, + \, g_{21} Z + \, g_{22} {Z}^2 \, + \, g_{23} {Z}^3 \, + \, g_{24} {Z}^4 \,+ \, g_{25} {Z}^5[$]
Which has an equivalent hermite representation as
[$]=gh_{20} \, + \, gh_{21} H_1(Z) + \, gh_{22} H_2(Z) \, + \, gh_{23} H_3(Z) \, + \, gh_{24} H_4(Z) \,+ \, gh_{25} H_5(Z)[$]
This makes our expression for second transition hermite equivalent to
[$]\sigma^2 \, dt^{2.0} \Big[ d_1 \, w(t)^{\beta_1} \, + \, d_2 \, w(t)^{\beta_2} \, \Big] \, ({Z_1}^2 \, - 1)\, [$]
[$]=\, \Big[\, gh_{20} \, + \, gh_{21} H_1(Z) + \, gh_{22} H_2(Z) \, + \, gh_{23} H_3(Z) \, + \, gh_{24} H_4(Z) \,+ \, gh_{25} H_5(Z) \, \Big]\, H_2(Z_1) [$]
In order to find one coefficient of second transition hermite, we calculate square root of variance of its hermite coefficienmtseries as
[$]G_{22}=\sqrt{\, \Big[\, {gh_{20}}^2 \, + \, {gh_{21}}^2 +2 \, {gh_{22}}^2 \, + 6 \, {gh_{23}}^2 \, +24 \, {gh_{24}}^2 \,+120 \, {gh_{25}}^2 \, \Big]}\,[$]
Now we can write the simple equation to advance the hermite representation of the SDE in two steps. First we add drift hermite series to SDE hermite series linearly as
[$]ch_0 \,=\,ch_0 +bh_0[$]
and similarly for rest of the higher hrmites are updated as (n represents coefficient of nth hermite)
[$]ch_n\,=\,ch_n +bh_n[$]
And then we add first and second transition hermite to the SDE hermite representation in a squared fashion.
[$]ch_1\,=\sqrt{\, {ch_1}^2 \,+ \, {G_{11}}^2 \,}[$]
and for adding the variance of second hermite, it is done as
[$]ch_2\,=\sqrt{\, {ch_2}^2 \,+ \, {G_{22}}^2 \,}[$]
Transition variance is not added to SDE hermites of order greater than two since 3rd transition hermite in Bessel coordinates takes coefficients in order dt^2.5 which become pretty much negligible for the kind of time step we are takin in our analytic simulation.
Now we can reconstruct the hermite representation of w(t) at next time step with updated hermite coefficients as
[$] \, w(t+1) \, = \,ch_0 \, + \, ch_1 H_1 \, + \, ch_2 \, H_2 \, + \, ch_3 \, H_3 \, + \, ch_4 \, H_4 \, + \, ch_5 \, H_5 \, [$]
Please note that we have to continue to switch our series representation from Z-series to hermite series and back and forth. Z-series is found around median of the density and can be used to find a series representation of any function of SDE variable. It is Z-series that we use to calculate the coefficients of drift term expansion and expansion of coefficient terms of first and second transition hermite. We cannot do a change of variable with hermite representation(which is around the mean) that is great for addition of variances something that we cannot do with a Z-series representation.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, mind control agents seem to be very annoyed with me and are playing all sort of pranks and tactics to torture me. It is 8:45 pm now and I slept earlier at 7:30 pm. During my brief sleep, they continued to make my arm numb and when I woke up, I had great difficulty moving it. I know they can do such things as they have done such things earlier several times. I still recall several years ago, I was lying on the carpet upside down with my tummy to the floor and working on the computer. I would have hardly remained in the stance for twenty minutes and they pulled all neurotransmitters out of my muscles of arms and legs. When I tried to stand up, it was almost impossible, I simply could not do that. Slowly I moved with sheer effort and then walked with extreme difficulty out of the drawing room where I was lying into the open garage and even then it took me almost half an hour to regain myself.
Just right now, they were trying something like that but I woke up prematurely to go into the wash room and my arm was almost completely numb and I had great difficulty moving it.
Later now I changed the position of cot in my room and tried sleeping again and then they started giving me serious itch in my back and it has become very difficult for me to sleep.
Again it seems that mind control crooks are very upset with me and want to use every tactic to unnerve me. I want to request friends to please force the American defense officials to play ugly and inhuman with me and end the torture.
Just right now, they were trying something like that but I woke up prematurely to go into the wash room and my arm was almost completely numb and I had great difficulty moving it.
Later now I changed the position of cot in my room and tried sleeping again and then they started giving me serious itch in my back and it has become very difficult for me to sleep.
Again it seems that mind control crooks are very upset with me and want to use every tactic to unnerve me. I want to request friends to please force the American defense officials to play ugly and inhuman with me and end the torture.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Today, I tried to solve for CEV noises and one drift SDEs in original coordinates by matching variances and it does not work. I get very similar shapes that I was getting when I tried to solve these SDEs in original coordinates by moment matching(which did not work for these SDEs).
So it seems that SDEs can be solved by matching variances only in Bessel coordinates only but once solved in Bessel framework, we can easily find the variances of the SDE and its functions in original coordinates by applying change of variable formula on our Z-series. And variances of orthogonal parts can be added in a squared fashion in original coordinates to retrieve the Variance of complete time evolution of SDE and do various analytics on it.
So it seems that SDEs can be solved by matching variances only in Bessel coordinates only but once solved in Bessel framework, we can easily find the variances of the SDE and its functions in original coordinates by applying change of variable formula on our Z-series. And variances of orthogonal parts can be added in a squared fashion in original coordinates to retrieve the Variance of complete time evolution of SDE and do various analytics on it.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
.To find the solution of density of SDE along constant CDF (or Constant Z-lines) lines, we suppose that we can express the solution of w(t) where w is the SDE/PDE variable as a power series in Z. After all the research we have done in past few months, it seems very natural that we express w(t) as a series in Z for any given time along time evolution of the density. We find the coordinates of this series at median of the density where Z=0. If we can represent w(t) at any given time as a valid series in Z, we can find all its derivatives and other smooth functions and also represent them as a series in Z. For a given evolution equation of the constant CDF lines of the density, we can represent all the variables in right hand side as a series in Z. The products, reciprocals and integrals of theses series again result as answers in another series in Z. So we find power series in Z for all variables in the right hand side of evolution equation and then do all mathematical operations on them in the form of series. The resulting answer is based on our original series at time [$]t_0[$] but generally has different series coefficients than that of series at time [$]t_0[$]. This new series that follows the application of SDE evolution equation at series coefficients at time [$]t_0[$] represents the solution of constant CDF lines of SDE density [$]w(t_1)[$] at next time level [$]t_1[$] as a power series in Z again. Now we take this power series at time [$]t_1[$] and find the series representations of its derivatives and other smooth functions and substitute them again in SDE evolution equation to find the series representation at time [$]t_2[$].
I will write the proposed power series for w(t) and its derivatives and then explain how to find the power series of smooth functions of w when power series for w is given. The series for w(t) is given as
[$]w(t)=a_0 \,+ a_1 \, Z + a_2 \, Z^2 + a_3 \, Z^3 + a_4 \, Z^4 +\, ... [$]
The first, second and third derivatives of w(t) w.r.t Z are given as
[$]\frac{\partial w(t)}{\partial Z}= a_1 \, + 2 \,a_2 \, Z +3 \, a_3 \, Z^2 +4 \, a_4 \, Z^3 +\, ... [$]
[$]\frac{\partial^2 w(t)}{\partial Z^2}= 2 \,a_2 \, +6 \, a_3 \, Z +12 \, a_4 \, Z^2 +\, ... [$]
[$]\frac{\partial^3 w(t)}{\partial Z^3}= 6 \, a_3 \, +24 \, a_4 \, Z \,+60 \, a_5 \, Z^2 +\, ... [$]
We do all the calculations of series evolution equation at median where Z=0. here all the derivatives are given only by leading coefficient associated with zeroth power of Z. So at Z=0,
[$]w(t)=a_0 \,[$]
[$]\frac{\partial w(t)}{\partial Z}= a_1 \, [$]
[$]\frac{\partial^2 w(t)}{\partial Z^2}= 2 \,a_2 \, [$]
[$]\frac{\partial^3 w(t)}{\partial Z^3}= 6 \, a_3 \, [$]
In order to represent smooth functions of w(t) as a series in Z, we suppose a form of the series and equate its coefficients with analytic derivatives w.r.t Z and this calculation is done at median that is Z=0.
So let us take drift in SDE, [$]\mu(w)[$] as a smooth function of w and represent it as a power series in Z. In what follows [$]\mu(w)[$] is considered a function operating on w just like we use symbol [$]f(w)[$] to indicate that f is a function operating on w. We suppose that functional form of [$]\mu(w)[$] and its derivatives is known analytically.
Let us suppose we know the form of power series of [$]\mu(w)[$] with the coefficients at yet unknown. This power series is
[$]\mu(w(t))=b_0 \,+ b_1 \, Z + b_2 \, Z^2 + b_3 \, Z^3 + b_4 \, Z^4 +\, ... [$]
The first, second and third derivatives of [$]\mu(w(t))[$] w.r.t Z are given as
[$]\frac{\partial \mu(w(t))}{\partial Z}= b_1 \, + 2 \,b_2 \, Z +3 \, b_3 \, Z^2 +4 \, b_4 \, Z^3 +\, ... [$]
[$]\frac{\partial^2 \mu(w(t))}{\partial Z^2}= 2 \,b_2 \, +6 \, b_3 \, Z +12 \, b_4 \, Z^2 +\, ... [$]
[$]\frac{\partial^3 \mu(w(t))}{\partial Z^3}= 6 \, b_3 \, +24 \, b4 \, Z \,+60 \, b_5 \, Z^2 +\, ... [$]
We know that values of above drift and its derivatives are given by leading coefficients of the series when Z=0. So at Z=0,
[$]\mu(w(t))=b_0 \,[$]
[$]\frac{\partial \mu(w(t))}{\partial Z}= b_1 \, [$]
[$]\frac{\partial^2 \mu(w(t))}{\partial Z^2}= 2 \,b_2 \, [$]
[$]\frac{\partial^3 \mu(w(t))}{\partial Z^3}= 6 \, b_3 \, [$]
We also know the value of [$]w(t)[$] at median which is [$]a_0[$]
but [$]\mu(w(t))=b_0 \,[$] which means that [$]b_0=\mu(a_0) \,[$]
To find the second coefficient, we equate first derivative of drift w.r.t Z at [$]w_0=a_0[$] with its coefficient. We know that
[$]\frac{\partial \mu(w(t))}{\partial Z} \,= \, \frac{\partial \mu(w(t))}{\partial w} \,\frac{\partial w(t)}{\partial Z} \,[$] which equals
[$]=\, \frac{\partial \mu(w(t))}{\partial w} \, a_1[$]
but [$]\frac{\partial \mu(w(t))}{\partial Z}= b_1 \, [$]
so
[$]b_1=\, \frac{\partial \mu(w(t))}{\partial w} \, a_1=\, \frac{\partial \mu(a_0)}{\partial w} \, a_1[$]
Similarly
[$]\frac{\partial^2 \mu(w(t))}{\partial Z^2} \,= \, \frac{\partial^2 \mu(w(t))}{\partial w^2} \, {(\frac{\partial w(t)}{\partial Z} )}^2+ \, \, \frac{\partial \mu(w(t))}{\partial w} \,\frac{\partial^2 w(t)}{\partial Z^2} \,[$]
[$]= \, \frac{\partial^2 \mu(w(t))}{\partial w^2} {a_1}^2 \, + \, \, \, \frac{\partial \mu(w(t))}{\partial w} (2a_2)[$]
We equate this with [$]\frac{\partial^2 \mu(w(t))}{\partial Z^2}= 2 \,b_2 \, [$] to find out that
[$]b_2= 1/2 \, \Big[\, \frac{\partial^2 \mu(w(t))}{\partial w^2} {a_1}^2 \, + \, \, \, \frac{\partial \mu(w(t))}{\partial w} (2a_2) \Big][$]
since above values of derivatives of w are calculated at [$]w=a_0[$], we can write
[$]b_2= 1/2 \, \Big[\, \frac{\partial^2 \mu(a_0)}{\partial w^2} {a_1}^2 \, + \, \, \, \frac{\partial \mu(a_0)}{\partial w} (2a_2) \Big][$]
similarly equating third derivative we get
[$]6 b_3= \, \frac{\partial^3 \mu(w)}{\partial w^3} {(\frac{\partial w(t)}{\partial Z} )}^3 \, + \,3 \, \frac{\partial^2 \mu(w)}{\partial w^2} \frac{\partial w(t)}{\partial Z} \frac{\partial^2 w(t)}{\partial Z^2} \, + \, \frac{\partial \mu(w)}{\partial w} \frac{\partial^3 w(t)}{\partial Z^3} \, [$]
which means
[$] b_3= 1/6 \, \Big[\, \frac{\partial^3 \mu(a_0)}{\partial w^3} {a_1}^3 \, + \,6 \, \frac{\partial^2 \mu(a_0)}{\partial w^2}\, a_1 \, a_2 + \, \frac{\partial \mu(a_0)}{\partial w} \, 6a_3 \, \Big][$]
So this way we continue to find the power series of smooth functions of w by equating the series coefficients for derivatives at Z=0 with their functional form there and using chain rule of derivatives.
I am very sure that similar series method can be used for a lot of other partial differential equations where solution remains smooth and well-behaved and series representation converges well.
In next relatively non-technical post, I will just give an idea how to relate to variables I have used in my matlab implementation to equations in this and previous post. In another later post, I will explain some enhancements I have used to get the series solution method to get to work better closer to zero.
.
https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=1290#p868484
.
In the above post, I discussed how to convert Z-series of an SDE variable to Z-series of functions of that variable.
Here are some interesting facts about Z-series of SDE variable and it functions and related values of their means and medians.
Suppose we are given Z-series of SDE variable w(t) as
[$]w(t)=a_0 \,+ a_1 \, Z + a_2 \, Z^2 + a_3 \, Z^3 + a_4 \, Z^4 +\, ... [$]
Obviously median of w(t) is given as
[$]\,Median[w(t)]\,=\,a_0\,[$]
From the copied post, you can find how to convert the above Z-series into Z-series of functions of w(t). Let us suppose we have calculated that Z-series for an arbitrary smooth function and it is given as
[$]f(w(t))=b_0 \,+ b_1 \, Z + b_2 \, Z^2 + b_3 \, Z^3 + b_4 \, Z^4 +\, ... [$]
Median of the w(t) is the constant term of the Z-series given as [$]a_0[$] while median of f(w(t)) is given as [$]b_0[$]
Both medians are directly related by the function as
[$]Median[f(w(t))] \,= \,b_0 \,= \,f(a_0) \,= \,f(Median[w(t)]) \,[$]
while for means we have
[$]Mean[w(t)]=E[w(t)]=a_0 \,+ a_2\, +\,3 \, a_4 \,+\,15 \, a_6 \,\,+ ... [$]
[$]Mean[w(t)]=Median[w(t)] \,+ a_2\, +\,3 \, a_4 \,+\,15 \, a_6 \,\,+ ... [$]
We have seen this in constant mean SDEs when variance increases and coefficients [$]a_2\,[$] , [$]a_4\,[$] , [$]a_6\,[$] increase in value and get larger, median has to decrease to keep the mean as constant.
Obviously there is no direct relationship between means of w(t) and that of its function f(w(t)). However for mean of f(w(t)), we have
[$]Mean[f(w(t))]=E[f(w(t))]=b_0 \,+ b_2\, +\,3 \, b_4 \,+\,15 \, b_6 \,\,+ ... [$]
[$]Mean[f(w(t))]=Median[f(w(t))] \,+ b_2\, +\,3 \, b_4 \,+\,15 \, b_6 \,\,+ ... [$]
So in order to find mean of functions of SDE variable w(t), we have to first find the Z-series of f(w(t)) (By method suggested in copied post) and then use this series to find the mean of f(w(t)).
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, during this week and next, we will be finding analytic solution of simple system of SDEs where asset price would be given by a Stochastic Volatility CEV SDE with one drift term and SV SDE would have mean-reverting dynamics. I plan to handle SDEs in Bessel coordinates. I hope to do this in next few days.
Early next week, I plan to solve for a basket of 5 correlated stocks having CEV dynamics with drift in a five factor model. Since we can handle many high-dimensional/basket models using several one dimensional SDEs, I hope that this new method will start to replace Monte Carlo to some extent.
Early next week, I plan to solve for a basket of 5 correlated stocks having CEV dynamics with drift in a five factor model. Since we can handle many high-dimensional/basket models using several one dimensional SDEs, I hope that this new method will start to replace Monte Carlo to some extent.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, I believe that problem of uncorrelated two SDEs SV pricing can possibly be solved in three different ways. I am trying to solve it first by matching variance but it is obviously more complicated than matching variances in simple 1D SDE in bessel coordinates. In order to properly take into account the effect of Stochastic volatility, I have expanded the asset in the form of a Z-series in two dimensions. One dimension in Z-series is in direction of asset and the other dimension in Z-series is in direction of stochastic volatility. I have also done work to convert two dimensional Z-series of the asset into two dimensional Z-series of functions of the asset. To advance the SDE, variances are added linearly (squared fashion for volatility) and for example first hermite polynomial of asset takes its coefficients in different hermite polynomials of stochastic volatility whose coefficients are updated by addition of coefficients of each stochastic volatility hermite polynomial in squared fashion and we keep track of coefficients of SV hermites (which form a series that makes coefficient of each asset hermite in two dimensional expansion).
Even though I have given brief details, I have still not been able to complete the work to solve the problem in full. Since only seeing is believing, there is still some possibility that something may not work in the method. I hope I will be able to share the results in a few days.
I also receive message that my IP is blacklisted and I cannot post the message. I then used a VPN to post the message. I have not posted anywhere else on any forum for several months. I do not send emails to anyone unknown and have not done that for several weeks.
Even though I have given brief details, I have still not been able to complete the work to solve the problem in full. Since only seeing is believing, there is still some possibility that something may not work in the method. I hope I will be able to share the results in a few days.
I also receive message that my IP is blacklisted and I cannot post the message. I then used a VPN to post the message. I have not posted anywhere else on any forum for several months. I do not send emails to anyone unknown and have not done that for several weeks.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, it seems that my first and simplest approach to analytic solution of SV SDEs does not work properly. It is close but I did not find it satisfactory enough. As I mentioned earlier that there are multiple ways to approach the problem. I am taking a different approach now but I really want the final result for the asset density to be in the form of a two dimensional Z-series. It is a lot of very interesting work and I hope to come back by the end of the week or over the weekend. I am very optimistic that soon we would have some very good results.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, I am writing this post to talk about recent mind control torture and related problems. As I had told friends earlier I still sleep in the car garage at night. About a week ago, I stopped sleeping in the garage and slept in my room for two or three days but they continued to release suffocating gases in my room and my sleep inside the room was of extremely poor quality. Then after three days of sleeping in my room, I started sleeping in the car garage again and my sleep is of relatively far better quality again.
Second major problem is that mind control agents continue to ask my family to put mind control drugs in my food. For past several weeks, almost everything I eat at home has slight amount of mind control drugs in it and food at home makes me very lethargic and it becomes very difficult to work after that. I am writing this post since two hours ago, I came from outside and I was very enthusiastic about my work but then I fried two eggs and had them with bread. There were mind control drugs either in oil, salt or eggs (They have ways to drug eggs with mind control drugs in entire markets) and just after eating food, I became extremely lethargic and all my energy and enthusiasm was gone and I dozed off for two hours and could not do any work.
I want to request good American friends to please stop these nasty crooks in American army and defense who bitterly want to continue my inhuman persecution because these crooks are desperate to enter the paradise of jewish Gold-father and they are very aware that jewish Gold-father is in no way obliged to hire them after they retire and it would happen only if they satisfactorily succeed in my persecution. Otherwise waste, these officers in American army are very aware that their future will be completely destroyed if Gold-father would not be pleased with them to give them a job in his companies.
Second major problem is that mind control agents continue to ask my family to put mind control drugs in my food. For past several weeks, almost everything I eat at home has slight amount of mind control drugs in it and food at home makes me very lethargic and it becomes very difficult to work after that. I am writing this post since two hours ago, I came from outside and I was very enthusiastic about my work but then I fried two eggs and had them with bread. There were mind control drugs either in oil, salt or eggs (They have ways to drug eggs with mind control drugs in entire markets) and just after eating food, I became extremely lethargic and all my energy and enthusiasm was gone and I dozed off for two hours and could not do any work.
I want to request good American friends to please stop these nasty crooks in American army and defense who bitterly want to continue my inhuman persecution because these crooks are desperate to enter the paradise of jewish Gold-father and they are very aware that jewish Gold-father is in no way obliged to hire them after they retire and it would happen only if they satisfactorily succeed in my persecution. Otherwise waste, these officers in American army are very aware that their future will be completely destroyed if Gold-father would not be pleased with them to give them a job in his companies.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, I was not able to work yesterday and day before yesterday on this two dimensional Z-series expansion of stochastic volatility derivatives. Mind control agents really made me unstable when I wrote my previous post and I was not able to work properly after that. I got back to work properly only today.
Here I want to give friends some idea about my work. I am using the analysis I earlier did on two dimensional fokker-planck equation. We have earlier solved one dimensional SDEs directly from analysis of one dimensional fokker-planck equation and that is when we started using Z-series in our analytics. I am just using a two dimensional Z-series to solve the equation derived from 2D Fokker-planck equation. I want it to be a general program where asset SDE can be very general with mean-reverting stochastic volatility.
I recall when I wrote the function of reciprocal of two dimensional series, mind control agents were surprised (Sometimes, they give enough expressions in my brain.) They were not expecting that I would be able to find reciprocal of the series. I had earlier written programs for functions of two dimensional series. After completing the work for functions of two dimensional series and reciprocal of 2D series, there was no conceptual problem in completing the work other than programming of some simple series algebra and once mind control agents realized that I was close to completing my work, they made me unstable with their drugs. I strongly suspect that they might have given my programs to someone in US who would claim that they have done everything on their own since I was not able to work for past two days and mind control agents were not suspecting that I would be able to restart today.
Good thing is that I am very close to completing the work and hope that I would be able to post the program tomorrow or on Sunday if mind control agents do not start any special antics again.
They were drugging food and water all along my path for morning walk earlier today. I leave for the walk a bit after 5:00 am and then walk for more than two hours. By the time I return, sun is already on fire and it is quite hot and I really need to buy some cold water or drinks during my walk. I was able to buy good water from a pharmacy where staff had no water or chillers in the customer area and the staff brought water from a fridge inside the pharmacy where customers cannot go. They had drugged all water where it was on display in chillers in customer area of pharmacies and other stores. It was so hot that I wanted cold water again and I tried water and juices at several other places and they were all drugged with mind control agencies.
Here I want to give friends some idea about my work. I am using the analysis I earlier did on two dimensional fokker-planck equation. We have earlier solved one dimensional SDEs directly from analysis of one dimensional fokker-planck equation and that is when we started using Z-series in our analytics. I am just using a two dimensional Z-series to solve the equation derived from 2D Fokker-planck equation. I want it to be a general program where asset SDE can be very general with mean-reverting stochastic volatility.
I recall when I wrote the function of reciprocal of two dimensional series, mind control agents were surprised (Sometimes, they give enough expressions in my brain.) They were not expecting that I would be able to find reciprocal of the series. I had earlier written programs for functions of two dimensional series. After completing the work for functions of two dimensional series and reciprocal of 2D series, there was no conceptual problem in completing the work other than programming of some simple series algebra and once mind control agents realized that I was close to completing my work, they made me unstable with their drugs. I strongly suspect that they might have given my programs to someone in US who would claim that they have done everything on their own since I was not able to work for past two days and mind control agents were not suspecting that I would be able to restart today.
Good thing is that I am very close to completing the work and hope that I would be able to post the program tomorrow or on Sunday if mind control agents do not start any special antics again.
They were drugging food and water all along my path for morning walk earlier today. I leave for the walk a bit after 5:00 am and then walk for more than two hours. By the time I return, sun is already on fire and it is quite hot and I really need to buy some cold water or drinks during my walk. I was able to buy good water from a pharmacy where staff had no water or chillers in the customer area and the staff brought water from a fridge inside the pharmacy where customers cannot go. They had drugged all water where it was on display in chillers in customer area of pharmacies and other stores. It was so hot that I wanted cold water again and I tried water and juices at several other places and they were all drugged with mind control agencies.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal
Re: Breakthrough in the theory of stochastic differential equations and their simulation
Friends, I completed the program for 2D SV SDEs just one and a half hour ago and have run it a few times. I get the shape of analytic density quite close to what I was getting from the simpler addition of variances method a few days ago. I am slightly concerned that the way I am constructing the final 1D Density of asset, my method to construct the density might be problematic and may be off. In next two days, I will try to construct 1D density of asset from a different method and see how it works. I might also write program to construct a 2D joint density of asset and stochastic vol from monte carlo and compare it to 2D density from the analytic method.
You think life is a secret, Life is only love of flying, It has seen many ups and downs, But it likes travel more than the destination. Allama Iqbal