NuclearPhynance Archive

Here is the new improved program for simulation of densities using Ito-Taylor density Simulation method I have worked on. It works quite well. But it is still not perfect for large terminal time with very high volatilities in the SDE. I have worked with squared orthogonal increments method and removed instabilities in the earlier versions and made the algorithm more accurate. Here is the program. It computes the density analytically and then makes a comparison with monte carlo density.





function [] = ItoTaylorDensityMRRecombining07()

%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 + mu2 x(t)^beta2 dt +sigma x(t)^gamma dz(t)

%I have not directly simulated the SDE but simulated the transformed

%Bessel process version of the SDE and then changed coordinates to retreive

%the SDE in original coordinates.

%The present program will directly simulate only the Bessel Process version of the

%SDE in transformed coordinates.

dt=.125/8; % Simulation time interval.%For diffusions close to zero

%decrease dt for accuracy.

Tt=8*4; % Number of simulation levels. Terminal time= Tt*dt; //.125/8*8*4=.5 years;

OrderA=2; %Keep at two.

OrderM=2; %Keep at two.

dNn=.2/1; % Normal density subdivisions width. would change with number of subdivisions

Nn=40; % No of normal density subdivisions

NnMidl=20;%One half left from mid of normal density

NnMidh=21;%One half right from the mid of normal density

NnMid=4.0;

x0=1.0; % starting value of SDE

beta1=0.0; % the first drift term power.

beta2=1.0; % Second drift term power.

gamma=.850; % volatility power.

kappa=.750; %mean reversion parameter.

theta=1.0; %mean reversion target

%Below you can get rid of kappa and theta and freely choose your own values

%for both drift coefficients.

mu1=1*theta*kappa; % first drift coefficient.

mu2=-1*kappa; % Second drift coefficient.

sigma0=.750;%Volatility value

alpha=1;% x^alpha is being expanded. This is currently for monte carlo only.

y(1:Nn)=x0;

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

Z(1:Nn)=(((1:Nn)-.5)*dNn-NnMid);



ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1);

ZProb(Nn)=1-normcdf(.5*Z(Nn)+.5*Z(Nn-1),0,1);

for nn=2:Nn-1

ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1)-normcdf(.5*Z(nn)+.5*Z(nn-1),0,1);

%Above calculate probability mass in each probability subdivision.

end



dt2=dt^2/2;

dtsqrt=sqrt(dt);

dtonehalf=dt.^1.5;



wnStart=1;



tic



for tt=1:Tt



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



theta1=mu1;

theta2=mu2;

theta3=.5*sigma0^2*(-gamma);

omega1=beta1-gamma;

omega2=beta2-gamma;

omega3=gamma-1;



GG(wnStart:Nn)=theta1*yy(wnStart:Nn).^omega1+theta2*yy(wnStart:Nn).^omega2+theta3*yy(wnStart:Nn).^omega3;

dGG(wnStart:Nn)=omega1*theta1*yy(wnStart:Nn).^(omega1-1)+omega2*theta2*yy(wnStart:Nn).^(omega2-1)+ ...

omega3*theta3*yy(wnStart:Nn).^(omega3-1);

DD(wnStart:Nn)=mu1*yy(wnStart:Nn).^beta1+mu2*yy(wnStart:Nn).^beta2;

dGG2(wnStart:Nn)=theta1*omega1*(omega1-1).*yy(wnStart:Nn).^(omega1-2)+ ...

theta2*omega2*(omega2-1).*yy(wnStart:Nn).^(omega2-2)+ ...

theta3*omega3*(omega3-1).*yy(wnStart:Nn).^(omega3-2);

QQ(wnStart:Nn)=sigma0^2*yy(wnStart:Nn).^(2*gamma);

VV(wnStart:Nn)=sigma0.*yy(wnStart:Nn).^gamma;

% d2GGVV(nn)=omega1*theta1*sigma0.*(omega1+gamma-1).*yy(nn).^(omega1+gamma-2)+ ...

% omega2*theta2*sigma0.*(omega2+gamma-1).*yy(nn).^(omega2+gamma-2)+ ...

% omega3*theta3*sigma0.*(omega3+gamma-1).*yy(nn).^(omega3+gamma-2);



wMu0dt(wnStart:Nn)=GG(wnStart:Nn)*dt + dGG(wnStart:Nn).*DD(wnStart:Nn).*dt2 + ...

.5*dGG2(wnStart:Nn).*QQ(wnStart:Nn).*dt2;



dw(wnStart:Nn)=sigma0*Z(wnStart:Nn).*dtsqrt+ ...

dGG(wnStart:Nn).*VV(wnStart:Nn).*(1-1/sqrt(3.0)).*Z(wnStart:Nn).*dtonehalf;%+ ...

%0*d2GGVV(nn).*VV(nn).*(1/sqrt(2.0).*(1-sqrt(2.0)/2.0)).*dt.^2.0.*(Z(nn).^2-1);



if(tt>0)

w(isnan(w)==1)=0;

wMu0dt(isnan(wMu0dt)==1)=0;

wPrev=w;

wMeanPrev=sum(ZProb(wnStart:Nn).*w(wnStart:Nn));%Calculate the mean.

w(wnStart:Nn)=wMeanPrev+wMu0dt(wnStart:Nn)+ ...

sign(w(wnStart:Nn)-wMeanPrev+dw(wnStart:Nn)).*sqrt(abs(sign(w(wnStart:Nn)-wMeanPrev).*(w(wnStart:Nn)-wMeanPrev).^2+ ...

sign(dw(wnStart:Nn)).*(dw(wnStart:Nn).^2)));



D1wl=(11/6*w(NnMidl-2)-3*w(NnMidl-3)+3/2*w(NnMidl-4)-1/3*w(NnMidl-5))/dNn;

D2wl=(35/12*w(NnMidl-2) -26/3*w(NnMidl-3)+ 19/2*w(NnMidl-4)-14/3*w(NnMidl-5)+11/12*w(NnMidl-6))/dNn.^2;

DeltaD1w=(w(NnMidh+2)-w(NnMidl-2)-5*D1wl*dNn-5.0*D2wl*dNn.^2.0/1.0)/12.50/dNn;



w(NnMidl-1)=w(NnMidl-2)+(.5*(D1wl+D2wl*dNn/1)+.5*(D1wl+D2wl*dNn/1+DeltaD1w))*dNn;

w(NnMidl)=w(NnMidl-1)+(.5*(D1wl+D2wl*dNn/1+DeltaD1w)+.5*(D1wl+D2wl*dNn/1+2*DeltaD1w))*dNn;

w(NnMidh)=w(NnMidl)+(.5*(D1wl+D2wl*dNn/1+2*DeltaD1w)+.5*(D1wl+D2wl*dNn/1+3*DeltaD1w))*dNn;

w(NnMidh+1)=w(NnMidh)+(.5*(D1wl+D2wl*dNn/1+3*DeltaD1w)+.5*(D1wl+D2wl*dNn/1+4*DeltaD1w))*dNn;





w1(1:Nn-1)=w(1:Nn-1);

w2(1:Nn-1)=w(2:Nn);

w(w1(:)>w2(:))=0;%Be careful;might not universally hold;



w(ww(nn+1))&&(EnterLoop1==1))

% w(1:nn)=0;

% wnStart=nn+1;

% EnterLoop1=0;

% end

% %Makes sure values are not smaller than zero.

% if((w(nn)<0)&&(EnterLoop2==1))

% w(1:nn)=0.0;

% wnStart=nn+1;

% EnterLoop2=0;

% end

% end

end



end



y_w(1:Nn)=0;

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

Dfy_w(wnStart:Nn)=0;

for nn=wnStart+1:Nn-1

Dfy_w(nn) = (y_w(nn + 1) - y_w(nn - 1))/(Z(nn + 1) - Z(nn - 1));

%Change of variable derivative for densities

end

fy_w(1:Nn)=0;

for nn = wnStart:Nn-1

fy_w(nn) = (normpdf(Z(nn),0, 1))/abs(Dfy_w(nn));%Origianl coordinates density

end



toc

str=input('Analytic Values calcuated; Press a key to start monte carlo');



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:Tt+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.

%Pre-compute the time and power exponent values in small multi-dimensional arrays

for k = 0 : (OrderA+1)

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) = dt^((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));

end

end

end

end

%Below call the program that calculates the Ito-Taylor coeffs.

YCoeff = ItoTaylorCoeffsNew(alpha,beta1,beta2,gamma);

%Calculate the density with monte carlo below.

rng(85109634, 'twister')

paths=100000;

YY(1:paths)=x0; %Original process monte carlo.

Random1(1:paths)=0;

for tt=1:Tt

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;

YY0=YY;

for k = 0 : OrderM

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;

if(l4<=5)



YY(1:paths)=YY(1:paths) + YCoeff(l1,l2,l3,l4) * ...

mu11(l1)*mu22(l2)*sigma22(l3)*sigma11(l4)* ...

YY0(1:paths).^Fp(l1,l2,l3,l4) .* HermiteP1(l4,1:paths) * Ft(l1,l2,l3,l4);

%Above is original coordinates monte carlo equation.

end

end

end

end

end

end

YY(YY<0)=0;

sum(YY(:))/paths %origianl coordinates monte carlo average.

sum(y_w(wnStart+1:Nn-1).*ZProb(wnStart+1:Nn-1)) %Original process average from coordinates

theta+(x0-theta)*exp(-kappa*dt*Tt)%Mean reverting SDE original variable true average



BinSize=.0075/2;%Please change this bin size accordingly. When data range is too small decrease it.

%When density range is too large increase it.

MaxCutOff=30;

[YDensity,IndexOutY,IndexMaxY] = MakeDensityFromSimulation_Infiniti(YY,paths,BinSize,MaxCutOff);

plot(y_w(wnStart+1:Nn-1),fy_w(wnStart+1:Nn-1),'r',IndexOutY(1:IndexMaxY),YDensity(1:IndexMaxY),'g');

str=input('Blue graph is from squared orthogonal incrments method and red is from squared volatility, green is monte carlo.');



end

As the friends would notice that the above program is quite accurate for reasonable volatilities but for very high volatilities, the results start to deteriorate after a year and more. On the other hand, well written PDEs give precise results even for extremely high volatilities. Therefore I decided to add PDEs to my existing research repertoire and substitute numerical solution of Fokker Planck where it is easily available and continue with the analytics. Once a precisely calculated density on a grid is available, it does not matter for further analytics if we obtained the grid from Ito-Taylor density generation method or from numerical solution of a PDE. I will continue to try to make the results of my Ito-Taylor density method sharply precise for extreme volatilities but in the mean time would also would do analytics with numerical solution of forward Fokker-Planck PDEs. Here is a post I made on Linkedin that elaborates on this idea.



https://www.linkedin.com/pulse/introducing-new-technologies-mathematical-finance-discipline-amin/



Introducing New Technologies in Mathematical Finance Discipline. A Totally New Paradigm.



I want to mention some interesting research achievements at Infiniti Derivatives Technologies in the area of stochastics leading to huge improvements in mathematical methods of derivatives valuation, hedging, trading and risk management.



1. Valuation of Path Dependent Stochastic Integrals on a PDE Grid.

It has been common historical wisdom that path dependent stochastic integrals cannot be calculated on a PDE grid. In our new research, we debunk these past notions and show that there are methods of introducing path dependence on a PDE grid. We have introduced new technologies that can calculate any stochastic path integral on a one or more dimensional PDE grid in fraction of a second after the PDE has been solved. Of course calculating path integrals on higher dimensional PDEs would require larger time but the time required to solve stochastic path integrals would strictly be far lesser than required to solve the PDE itself. Some of these stochastic integrals could be for example conditional value of path integral of squared volatility on the PDE grid for valuation of derivatives with stochastic volatility. Once the PDE has been solved conditional backward stochastic integrals like exponential of path integral of interest rates representing conditional bond prices or discount factor could also be calculated on the same PDE grid. From the PDE of interest rates, you could calculate the backward evolution of difficult stochastic path integrals depending on interest rates of whose PDE cannot be written easily. For Asian options, you could find the path integral of underlying as a forward stochastic integral on the Fokker planck equation solved by a one dimensional PDE or on a two dimensional PDE with stochastic volatility. The conditional values of these stochastic path integrals are available anywhere on the PDE grid along the backward or forward evolution of the PDE. PDE methods can only calculate backward evolution of payoffs whose explicit PDE has been written with arbitrary final data but we can calculate backward dynamic programming valuation of far more complex stochastic integrals dependent on PDE underlying of whose explicit PDE is not available but conditional of the value of the underlying as calculated in the forward PDE. We do not use any simulation or random numbers and only rely on fast analytics to achieve the above mentioned goals.



2. Calculation of Optimal Quantizers for Monte Carlo Simulations.

Once a Fokker planck PDE has been solved, we can calculate optimal quantizers in one or two dimensions in fraction of a second. You could use these optimal quantizer schemes for forward evolution of stochastic integrals or backward dynamic programming in a monte carlo setting. As opposed to existing optimal quantizer methods that are based on monte carlo discretization schemes coupled with optimization techniques, we take a Fokker planck PDE solution as a starting point which is more cheaply available and simply use very fast proprietary analytics to generate optimal quantizer grid that can later be used for extremely fast monte carlo simulations and backward dynamic programming.



If you find any of this interesting and would like to implement this solution at your institution, you can email me at anan2999yahoocom for further information and discussions. I want to mention again that none of these new technologies are in public domain.

Download Matlab Code for fast Calculation of Densities of Stochastic Integrals on a lattice using our new method.



Here, in post 793 on this thread on W****tt forum, you can download the matlab code for calculation of densities of stochastic integrals.



https://forum.W****tt.com/viewtopic.php?f=4&t=99702&start=780

Download Matlab Code for fast Calculation of Densities of Mean Reverting SDEs with high volatility and Mean Reversion using our updated method with Second Order Correction.



You can use the earlier code now for evolution of SDEs for calculation of densities of mean reverting SDEs where earlier code was failing. We have made a second order correction in the earlier code. Now the code perfectly works for fast mean reverting SDEs. Even models like heston where the density goes into zero seem to work quite reasonably well.



Here, in post 797 on my thread in technical forum on W****tt, you can download the matlab code for calculation of densities of mean reverting SDEs. https://forum.W****tt.com/viewtopic.php?f=4&t=99702&start=795