amin


Total Posts: 281 
Joined: Aug 2005 


Here is the new improved program for simulation of densities using ItoTaylor 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/ahsanamin0a53334 %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^(1gamma)/(1gamma); Z(1:Nn)=(((1:Nn).5)*dNnNnMid);
ZProb(1)=normcdf(.5*Z(1)+.5*Z(2),0,1); ZProb(Nn)=1normcdf(.5*Z(Nn)+.5*Z(Nn1),0,1); for nn=2:Nn1 ZProb(nn)=normcdf(.5*Z(nn)+.5*Z(nn+1),0,1)normcdf(.5*Z(nn)+.5*Z(nn1),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)=((1gamma)*w(wnStart:Nn)).^(1/(1gamma));
theta1=mu1; theta2=mu2; theta3=.5*sigma0^2*(gamma); omega1=beta1gamma; omega2=beta2gamma; omega3=gamma1; 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).^(omega11)+omega2*theta2*yy(wnStart:Nn).^(omega21)+ ... omega3*theta3*yy(wnStart:Nn).^(omega31); DD(wnStart:Nn)=mu1*yy(wnStart:Nn).^beta1+mu2*yy(wnStart:Nn).^beta2; dGG2(wnStart:Nn)=theta1*omega1*(omega11).*yy(wnStart:Nn).^(omega12)+ ... theta2*omega2*(omega21).*yy(wnStart:Nn).^(omega22)+ ... theta3*omega3*(omega31).*yy(wnStart:Nn).^(omega32); QQ(wnStart:Nn)=sigma0^2*yy(wnStart:Nn).^(2*gamma); VV(wnStart:Nn)=sigma0.*yy(wnStart:Nn).^gamma; % d2GGVV(nn)=omega1*theta1*sigma0.*(omega1+gamma1).*yy(nn).^(omega1+gamma2)+ ... % omega2*theta2*sigma0.*(omega2+gamma1).*yy(nn).^(omega2+gamma2)+ ... % omega3*theta3*sigma0.*(omega3+gamma1).*yy(nn).^(omega3+gamma2); 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).*(11/sqrt(3.0)).*Z(wnStart:Nn).*dtonehalf;%+ ... %0*d2GGVV(nn).*VV(nn).*(1/sqrt(2.0).*(1sqrt(2.0)/2.0)).*dt.^2.0.*(Z(nn).^21); 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(NnMidl2)3*w(NnMidl3)+3/2*w(NnMidl4)1/3*w(NnMidl5))/dNn; D2wl=(35/12*w(NnMidl2) 26/3*w(NnMidl3)+ 19/2*w(NnMidl4)14/3*w(NnMidl5)+11/12*w(NnMidl6))/dNn.^2; DeltaD1w=(w(NnMidh+2)w(NnMidl2)5*D1wl*dNn5.0*D2wl*dNn.^2.0/1.0)/12.50/dNn;
w(NnMidl1)=w(NnMidl2)+(.5*(D1wl+D2wl*dNn/1)+.5*(D1wl+D2wl*dNn/1+DeltaD1w))*dNn; w(NnMidl)=w(NnMidl1)+(.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:Nn1)=w(1:Nn1); w2(1:Nn1)=w(2:Nn); w(w1(:)>w2(:))=0;%Be careful;might not universally hold; w(w<0)=0.0; for nn=1:Nn if(w(nn)<=0) wnStart=nn; end end % EnterLoop1=1; % EnterLoop2=1; % for nn=Nn1:1:wnStart % %Makes sure values decrease monotonically close to zero. % if((w(nn)>w(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) = ((1gamma)*w(wnStart:Nn)).^(1/(1gamma)); Dfy_w(wnStart:Nn)=0; for nn=wnStart+1:Nn1 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:Nn1 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^(k1); end if mu1 ~= 0 mu11(k)=mu1^(k1); end if mu2 ~= 0 mu22(k)=mu2^(k1); end if sigma0~=0 sigma22(k)=sigma0^(2*(k1)); 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. %Precompute the time and power exponent values in small multidimensional 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^((l11) + (l21) + (l31) + .5* (l41)); Fp(l1,l2,l3,l4) = (alpha + (l11) * beta1 + (l21) * beta2 + (l31) * 2* gamma + (l41) * gamma ...  (l11)  (l21)  2* (l31)  (l41)); end end end end %Below call the program that calculates the ItoTaylor 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).^21; 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:Nn1).*ZProb(wnStart+1:Nn1)) %Original process average from coordinates theta+(x0theta)*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:Nn1),fy_w(wnStart+1:Nn1),'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 



amin


Total Posts: 281 
Joined: Aug 2005 


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 ItoTaylor density generation method or from numerical solution of a PDE. I will continue to try to make the results of my ItoTaylor density method sharply precise for extreme volatilities but in the mean time would also would do analytics with numerical solution of forward FokkerPlanck PDEs. Here is a post I made on Linkedin that elaborates on this idea.
https://www.linkedin.com/pulse/introducingnewtechnologiesmathematicalfinancedisciplineamin/
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. 

