These are very initial programs for review by friends and many changes are required. I want to improve the programs in several ways. I also have some ideas about an alternative algorithm for construction of analytic density and associated 2D Hermite series. But just sharing the initial programs.

I have given parameters for maximum values for asset and volatility and please use them freely to not let the grid become large and keep it manageable. This can be changed once a non-uniform grid is introduced. If the paths are not too large, keep the grid coarse for proper construction of numerical density.

When volatility goes smoothly into zero and a lot of mass is at zero, there is a mismatch between analytic and numerical density at zero. I am currently using CEV power .95 for volatility. I am sure this can be improved later.

I hope that this is just start of the work and there is a lot of good research work to follow.

Here is the program.

```
function [] = SVMonteCarloDensityWith2DHermitesWilmott02()
%Copyright Ahsan Amin. Infiniti derivatives Technologies.
%Please fell free to connect on linkedin: linkedin.com/in/ahsan-amin-0a53334
%Volatility and asset equations are
%Vt(0)=V0;
%dVt=(mu1*Vt^beta1 + mu2*Vt^beta2)*dt+sigma0 * Vt^gamma * dz1(t)
%Xt(0)=X0
%dXt=(a*Xt^alpha1+b*Xt^alpha2)*dt+sigma1* Vt^gammaV * Xt^gammaX *dz2(t);
%<dz1(t),dz2(t)>=rho *dt
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
X0=1.0;
alpha1=0;
alpha2=1;
kappaX=0;
thetaX=1.0;
a=kappaX*thetaX;
b=-kappaX;
rho=0;
sigma1=1;
gammaV=.5;
gammaX=.85;
%%%%%%%%`%%%%%
V0=.25;
thetaV=.06;
kappaV=1.50;
mu1=kappaV*thetaV;
mu2=-kappaV;
beta1=0;
beta2=1;
gamma=.95;
sigma0=.80;
%%%%%%%%%%%%
dt=.03125/2;
Tt=64*4;
T=Tt*dt;
%%%%%%%%%%%%%%%%%%%%%%%%5
%rng(29079137, 'twister')
rng(15898837, 'twister')
paths=500000;
V(1:paths)=V0; %Original process monte carlo.
X(1:paths)=X0;
Random1(1:paths)=0;
Random2(1:paths)=0;
for tt=1:Tt
Random1=randn(size(Random1));
Random2=randn(size(Random2));
X(1:paths)=X(1:paths)+ ...
(a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dt + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) * sqrt(dt) + ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*sqrt(dt) + ...
(a*alpha1* X(1:paths).^(alpha1-1)+b*alpha2* X(1:paths).^(alpha2-1)).* ...
(((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2)* dt^2/2)+ ...
(sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths) *(1-1/sqrt(3)).*dt^1.5+ ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random2(1:paths)*(1-1/sqrt(3)).*dt^1.5))+ ...
.5*(a*alpha1*(alpha1-1)* X(1:paths).^(alpha1-2)+b*alpha2*(alpha2-1).* X(1:paths).^(alpha2-2)).* ...
( sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX)) *dt^2/2 + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random1(1:paths) * 1/sqrt(3).* dt^1.5 + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*(Random1(1:paths).^2-1) * dt/2 + ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths)*dt/2)+ ...
.5*sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
(sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random1(1:paths).*1/sqrt(3).*dt^1.5)+ ...
sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random1(1:paths) * 1/sqrt(3).* dt^1.5 + ...
sigma0*V(1:paths).^gamma.*Random1(1:paths).*Random2(1:paths)*dt/2)+ ...
.5*sqrt(1-rho^2)* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
(sigma0^2*V(1:paths).^(2*gamma).*Random1(1:paths)*1/sqrt(3)*dt^1.5)+ ...
sqrt(1-rho^2)* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random1(1:paths)*1/sqrt(3)*dt^1.5+ ...
rho* sigma1* V(1:paths).^gammaV.*gammaX.* X(1:paths).^(gammaX-1).* ...
((a* X(1:paths).^alpha1 + b* X(1:paths).^alpha2).*Random2(1:paths) * 1/sqrt(3).* dt^1.5 + ...
sqrt(1-rho^2)* sigma1* V(1:paths).^gammaV.* X(1:paths).^gammaX .*Random1(1:paths).*Random2(1:paths) * dt/2 + ...
rho* sigma1* V(1:paths).^gammaV .*X(1:paths).^gammaX .*(Random2(1:paths).^2-1)*dt/2)+ ...
.5*rho* sigma1* V(1:paths).^gammaV.*gammaX.*(gammaX-1).* X(1:paths).^(gammaX-2).* ...
(sigma1^2* V(1:paths).^(2*gammaV).* X(1:paths).^(2*gammaX) .*Random2(1:paths).*1/sqrt(3).*dt^1.5)+ ...
rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).* X(1:paths).^(gammaX).* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths) * 1/sqrt(3).* dt^1.5 + ...
sigma0*V(1:paths).^gamma.*(Random2(1:paths).^2-1)*dt/2)+ ...
.5*rho* sigma1*gammaV.*(gammaV-1).* V(1:paths).^(gammaV-2).* X(1:paths).^(gammaX).* ...
sigma0^2.*V(1:paths).^(2*gamma).*Random2(1:paths) * 1/sqrt(3).* dt^1.5+ ...
rho* sigma1*gammaV.* V(1:paths).^(gammaV-1).*gammaX.* X(1:paths).^(gammaX-1).* ...
rho.* sigma1.* V(1:paths).^gammaV .*X(1:paths).^gammaX .*sigma0.*V(1:paths).^gamma.*Random2(1:paths)*1/sqrt(3)*dt^1.5;
V(1:paths)=V(1:paths)+ ...
(mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dt + ...
sigma0*V(1:paths).^gamma .*Random2(1:paths)*sqrt(dt) + ...
(mu1.*beta1*V(1:paths).^(beta1-1) + mu2.*beta2.*V(1:paths).^(beta2-1)).* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2)*dt^2/2 + ...
sigma0*V(1:paths).^gamma .*Random2(1:paths)*(1-1/sqrt(3))*dt^1.5) + ...
.5*(mu1.*beta1.*(beta1-1).*V(1:paths).^(beta1-2) + mu2.*beta2.*(beta2-1).*V(1:paths).^(beta2-2)).* ...
sigma0^2.*V(1:paths).^(2*gamma).*dt^2/2 + ...
sigma0*gamma*V(1:paths).^(gamma-1) .* ...
((mu1.*V(1:paths).^beta1 + mu2.*V(1:paths).^beta2).*Random2(1:paths).*1/sqrt(3)*dt^1.5 + ...
sigma0.*V(1:paths).^gamma .*(Random2(1:paths).^2-1)*dt/2) + ...
.5*sigma0*gamma*(gamma-1).*V(1:paths).^(gamma-2) .* ...
sigma0^2.*V(1:paths).^(2*gamma) .*Random2(1:paths).*1/sqrt(3)*dt^1.5;
end
X(X<0)=0;
V(V<0)=0;
X=real(X);
V=real(V);
SVolMeanAnalytic=thetaV+exp(-kappaV*T)*(V0-thetaV)
SVolMean=sum(V(1:paths))/paths
AssetMeanAnalytic=thetaX+exp(-kappaX*T)*(X0-thetaX)
AssetMean=sum(X(1:paths))/paths
str=input("Look at monte carlo Mean of X and V")
% NoOfBins=300;
% MaxCutOff=300;
% [XDensity0,IndexOutX0,IndexMaxX0] = MakeDensityFromSimulation_Infiniti_NEW(X,paths,NoOfBins,MaxCutOff );
% plot(IndexOutX0(1:IndexMaxX0),XDensity0(1:IndexMaxX0),'g');
%
% str=input("Look at 1D graph of X density");
NoOfBinsV=75;
NoOfBinsX=150;
MaxCutOffX=7;
MaxCutOffV=0.65;
%[XYProbDensity2D,XYProbMass2D,XProbMass,VProbMass,IndexMaxX,IndexMaxV,Xn,Vm,Zx,Zv]= Find2DDensityAndCDFFromSimulation_Infiniti(X,V,paths,NoOfBinsX,NoOfBinsV);
[XYProbDensity2D,XYProbMass2D,XProbMass,VProbMass,IndexMaxX,IndexMaxV,Xn,Vm,Zx,Zv,BinSizeX,BinSizeV]= Find2DDensityAndCDFFromSimulation_Infiniti04(X,V,paths,NoOfBinsX,NoOfBinsV,MaxCutOffX,MaxCutOffV);
%[XYProbDensity2D,XYProbMass2D,XProbMass,YProbMass,IndexMaxX,IndexMaxY,Xn,Ym,Zx,Zy,BinSizeX,BinSizeY]
%Please look at matlab file of above function for etailed explanation of
%all its input and output parameters.
surf(Xn,Vm,transpose(XYProbDensity2D))
xlabel('Asset'), ylabel('Stochastic Volatility'), zlabel('Probability Density 2D')
title('Graph of 2D Density of Asset and SV from Parzen Density')
MeanX=sum(X(1:paths))/paths
MeanV=sum(V(1:paths))/paths
str=input("Look at 3D graphs of the numerical density")
%Below calculate Coeffs of 2D Hermite products.
%Oxn is order of Zx hermites. ranges from 1 to 8 in indexing or (0 to 7 in reality)
%Ovm is order of Zv hermites. ranges from 1 to 8 in indexing or (0 to 7 in reality)
%HermiteH(Oxn,Zx(nn,mm)) is a small function that caclulates Oxn order hermite at Zx(nn,mm) value
%HermiteH(Ovm,Zv(mm)) caclulates Ovm order hermite at Zy(mm) value
%Below Eq calculates 2D Coeffs with inner product of coordinates of 2D uniform grid with
%product of X and V hermites.
C(1:8,1:8)=0;
for Oxn=1:8
for Ovm=1:8
for nn=1:(NoOfBinsX+1)
for mm=1:(NoOfBinsV+1)
C(Oxn,Ovm)=C(Oxn,Ovm)+1/factorial(Oxn-1)/factorial(Ovm-1).*Xn(nn).*Vm(mm).*HermiteH(Oxn-1,Zx(nn,mm)).*HermiteH(Ovm-1,Zv(mm)).*XYProbMass2D(nn,mm);
%C(Oxn,Ovm)=C(Oxn,Ovm)+1/factorial(Oxn-1)/factorial(Ovm-1).*Xn(nn).*HermiteH(Oxn-1,Zx(nn,mm)).*HermiteH(Ovm-1,Zv(mm)).*XYProbMass2D(nn,mm);
end
end
end
end
C
%If C are 2D hermite Coeffs as calculated in above equation
%Hermite Series representation of 2D random variable would be given by
%double sum equation below as
%Sum_{n=0}^{n=7} Sum_{m=0}^{m=7} C(n+1,m+1) He(n,Zx) He(m,Zy);
str=input("Look at Coefficients of 2D Hermite products")
%Below we calculate 1D Hermite Series for volatility.
%We only need to calculate the mean (zeroth hermite) aH0 of the volatility
%random variable.
aH0=0;
aH(1:7)=0;
for nn=1:IndexMaxV
aH0 = aH0+Vm(nn).*VProbMass(nn);
aH(1)=aH(1)+Vm(nn).*Zv(nn).*VProbMass(nn);
aH(2)=aH(2)+.5*Vm(nn).*(Zv(nn).^2-1).*VProbMass(nn);
aH(3)=aH(3)+1/6*Vm(nn).*(Zv(nn).^3-3*Zv(nn)).*VProbMass(nn);
aH(4)=aH(4)+1/24*Vm(nn).*(Zv(nn).^4-6*Zv(nn).^2+3).*VProbMass(nn);
aH(5)=aH(5)+1/120*Vm(nn).*(Zv(nn).^5-10*Zv(nn).^3+15*Zv(nn)).*VProbMass(nn);
aH(6)=aH(6)+1/720*Vm(nn).*(Zv(nn).^6-15*Zv(nn).^4+45*Zv(nn).^2-15).*VProbMass(nn);
aH(7)=aH(7)+1/720/7*Vm(nn).*(Zv(nn).^7-21*Zv(nn).^5+105*Zv(nn).^3-105*Zv(nn)).*VProbMass(nn);
end
[a0,a] = ConvertHCoeffsToZCoeffs(aH0,aH,7);
a0
a
str=input("Look at volatility coefficients ----a0,a(1:7)");
%Below we calculate the Hermite series of asset by integrating over the
%volatility to get marginal series of asset. When we integrate over
%volatility, all volatility hermites(alone or as products with X-hermites)
%in the 2D asset series go to zero(they all have mean zero).
%So basically all volatility terms go to zero and we are left with
%one dimensional array containing terms from 2D series where
%volatility hermites in the product had order zero. This order zero is
%represented by volatility index one in the 2d array of coeffs C(n,m);
%So our 1D asset series is from the line C(1:8,1) where volatilit hermites
%are zero. This is very interesting.
bH0=C(1,1);
bH(1:7)=C(2:8,1);
%Above bH0 is asset Hermite series coeff of zeroth order hermite and
%bH(1:7) represent hermites from 1st to 7th order.
bH0=bH0/aH0;
bH=bH/aH0;
%Above, we need to normalize the asset hermite series with expected value
%of volatility since we integrated over the volatility. Just like in
%Conditioning when we integrate over a variable to get marginal density and
%then we have to divide with the integral of the integrated density as
%well. Similalry, Here we have to divide with the mean/expected value of
%volatility hermite series(which is integral of the volatility hermite series)
%to normalize our asset equation hermite coefficients.
%Below we convert from hermite series to Z-series of the asset for
%plotting.
[b0,b] = ConvertHCoeffsToZCoeffs(bH0,bH,7);
bH0
bH
b0
b
str=input("Look at 1D Coeffs of X")
[pXY2D] = Calculate2VariablesOn2DGridFrom2DZSeries(Xn,Vm,NoOfBinsX,NoOfBinsV,C,a0,a,BinSizeX,BinSizeV);
Xn0(1:NoOfBinsX-1)=Xn(2:NoOfBinsX);
Vm0(1:NoOfBinsV-1)=Vm(2:NoOfBinsV);
pXY2D0(1:NoOfBinsX-1,1:NoOfBinsV-1)=pXY2D(2:NoOfBinsX,2:NoOfBinsV);
%surf(Xn,Vm,transpose(pXY2D))
surf(Xn0,Vm0,transpose(pXY2D0))
xlabel('Asset'), ylabel('Stochastic Volatility'), zlabel('Probability Density 2D')
title('Graph of 2D Density of Asset and SV From Z-series')
str=input("Look at graph of Analytic 2D density");
XYProbDensity2D0(1:NoOfBinsX-1,1:NoOfBinsV-1)=XYProbDensity2D(2:NoOfBinsX,2:NoOfBinsV);
surf(Xn0,Vm0,transpose(pXY2D0-XYProbDensity2D0))
xlabel('Asset'), ylabel('Stochastic Volatility'), zlabel('Difference of Probability Density 2D')
title('Graph of 2D Difference of Parzen Density and Z-series Analytic Density of Asset and SV on the same grid')
str=input("Look at difference of Analytic and Numerical density graphs")
end
```

.

.

.

```
function [XYProbDensity2D,XYProbMass2D,XProbMass,YProbMass,IndexMaxX,IndexMaxY,Xn,Ym,Zx,Zy,BinSizeX,BinSizeY] = Find2DDensityAndCDFFromSimulation_Infiniti02(X,Y,Paths,NoOfBinsX,NoOfBinsY,MaxCutOffX,MaxCutOffY)
%This function processes joint observations of two stochastic variables X and Y from
%common monte carlo paths and makes a two dimensional density of X and Y.
%Explanation of Input Parameters
%X and Y are 1D arrays so that nth element of each makes a joint
%observation of nth monte carlo path X(n) and Y(n) form joint nth observation of two stochastic variables.
%paths ---- total monte carlo paths.
%NoOfBinsX ---- Uniform Grid size in X for density construction. You can alter
%this. But you will have to make sure to alter grid sizes in main program after
%this function has returned results.
%NoOfBinsY ---- Uniform Grid size in Y for density construction. You can alter
%this.
%Explanation of output Parameters
%XYProbDensity2D --- is 2D numerically constructed Parzen density on input
%grid in X and Y.
%XYProbMass2D, ------ is integrated version of 2D parzen density that gives
%two D probability mass on input grid in X and Y.
%XProbMass ---- is probability mass of 1D marginal density in X on NoOfBinsX Grid.
%YProbMass ---- is probability mass of 1D marginal density in Y on NoOfBinsY Grid.
%IndexMaxX ----- This is size of grid in X and corresponds to NoOfBinsX and
%is NoOfBinsX+1. Naturally This is also size of array that returns
%X-values on the X-axis of 2D uniform grid.
%IndexMaxY ----- This is size of grid in X and corresponds to NoOfBinsY and
%is NoOfBinsY+1. Naturally This is also size of array that returns
%Y-values on Y-axis of the 2D uniform grid.
%Xn--- is an array of size IndexMaxX that contains values of X at the
%center of each grid cell on X-axis on 2D uniform grid.
%Ym--- is an array of size IndexMaxY that contains values of Y at the
%center of each grid cell on Y-axis on 2D uniform grid.
%Zx----- this is a 2D array. This returns the value of standard normal associated with Xn for every
%Ym grid cell. Given a certain grid cell Ym, we calculate CDF along all X grid
%cells with the same Ym. This is derived from inversion of univariate CDF on X-grid
%or X-axis associated with each Y grid cell. Zx(n,m) is a 2D array.
%Zy----- this is the value of standard normal associated with Ym at center
%of each grid cell on Y-axis. This is derived from inversion of univariate CDF on Y-grid
%or Y-axis. This is 1D array. This is used for vol which is independent of
%X.
meanX=sum(X(:))/Paths;
meanY=sum(Y(:))/Paths;
sdX=sqrt(sum((X(:)-meanX).^2));
sdY=sqrt(sum((Y(:)-meanY).^2));
Xmin=99999;
Xmin1=999998;
Xmin2=999997;
Xmax=0;
Xmax1=1*10^(-12);
Xmax2=1*10^(-11);
Ymin=999999;
Ymin1=999998;
Ymin2=999997;
Ymax=0;
Ymax1=1*10^(-12);
Ymax2=1*10^(-11);
%meanX=0;
%meanY=0;
for p=1:Paths
if(isreal(X(p))==0)
X(p)
str=input("Look at X(p)");
end
if(isreal(Y(p))==0)
Y(p)
str=input("Look at Y(p)");
end
if(X(p)>MaxCutOffX)
X(p)=MaxCutOffX;
end
if(Y(p)>MaxCutOffY)
Y(p)=MaxCutOffY;
end
if(Xmin>real(X(p)))
Xmin2=Xmin1;
Xmin1=Xmin;
Xmin=real(X(p));
elseif(Xmax<real(X(p)))
Xmax2=Xmax1;
Xmax1=Xmax;
Xmax=real(X(p));
end
if(Ymin>real(Y(p)))
Ymin2=Ymin1;
Ymin1=Ymin;
Ymin=real(Y(p));
elseif(Ymax<real(Y(p)))
Ymax2=Ymax1;
Ymax1=Ymax;
Ymax=real(Y(p));
end
end
%Xmin
%Xmin1
%Xmin2
%Xmax
%Xmax1
%Xmax2
%str=input('Look at Xmin and Xmax');
%IndexMax=NoOfBins+1;
BinSizeX=(Xmax2-Xmin2)/NoOfBinsX;
BinSizeY=(Ymax2-Ymin2)/NoOfBinsY;
%IndexMax=floor((Xmax-Xmin)/BinSize+.5)+1
IndexMaxX=floor((Xmax2-Xmin2)/BinSizeX+.5)+1;
IndexMaxY=floor((Ymax2-Ymin2)/BinSizeY+.5)+1;
BinSizeX
BinSizeY
IndexMaxX
IndexMaxY
Xmax2
Xmax1
Xmax
Xmin2
Xmin1
Xmin
Ymax2
Ymax1
Ymax
Ymin2
Ymin1
Ymin
str=input("Look at max-min parameters of numerical density")
%IndexOut(1:IndexMax)=Xmin+(0:(IndexMax-1))*BinSize;
Xn(1:IndexMaxX)=Xmin2+(0:(IndexMaxX-1))*BinSizeX;
Ym(1:IndexMaxY)=Ymin2+(0:(IndexMaxY-1))*BinSizeY;
XYProbMass2D(1:IndexMaxX,1:IndexMaxY)=0.0;
XYProbDensity2D(1:IndexMaxX,1:IndexMaxY)=0.0;
XProbMass(1:IndexMaxX)=0;
YProbMass(1:IndexMaxY)=0;
for p=1:Paths
indexX=real(floor(real(X(p)-Xmin2)/BinSizeX+.5)+1);
%indexXp=real(X(p)-Xmin2)/BinSizeX+1;
if(real(indexX)<1)
indexX=1;
end
if(real(indexX)>IndexMaxX)
indexX=IndexMaxX;
end
indexY=real(floor(real(Y(p)-Ymin2)/BinSizeY+.5)+1);
%indexYp=real(Y(p)-Ymin2)/BinSizeY;
if(real(indexY)<1)
indexY=1;
end
if(real(indexY)>IndexMaxY)
indexY=IndexMaxY;
end
%XDensity(index)=XDensity(index)+1.0/Paths/BinSize;
XProbMass(indexX)=XProbMass(indexX)+1.0/Paths;
YProbMass(indexY)=YProbMass(indexY)+1.0/Paths;
XYProbMass2D(indexX,indexY)=XYProbMass2D(indexX,indexY)+1.0/Paths;
XYProbDensity2D(indexX,indexY)=XYProbDensity2D(indexX,indexY)+1.0/Paths/BinSizeX/BinSizeY;
end
%BinSizeX
%BinSizeY
%str=input("Look at BinSizes");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%Below calculate Zx by inversion of 1D X-CDF associated with each volatility grid cell
%so that resulting array is 2D.
XProbSum(1:IndexMaxX,1:IndexMaxY)=0;
XCDF(1:IndexMaxX,1:IndexMaxY)=0;
Zx(1:IndexMaxX,1:IndexMaxY)=0;
for mm=1:IndexMaxY
XProbSum(1,mm)=XYProbMass2D(1,mm)*.5;
if(YProbMass(mm)>0)
XCDF(1,mm)=XProbSum(1,mm)./YProbMass(mm);
else
XCDF(1,mm)=0;
end
if((XCDF(1,mm)>0) && (XCDF(1,mm)<1))
Zx(1,mm)=norminv(XCDF(1,mm));
elseif (XCDF(1,mm)<.000000001)
Zx(1,mm)=-5;
elseif (XCDF(1,mm)>.999999999)
Zx(1,mm)=5;
end
for nn=2:IndexMaxX
XProbSum(nn,mm)=XProbSum(nn-1,mm)+XYProbMass2D(nn-1,mm)*.5+XYProbMass2D(nn,mm)*.5;
if(YProbMass(mm)>0)
XCDF(nn,mm)=XProbSum(nn,mm)./YProbMass(mm);
else
XCDF(nn,mm)=0;
end
if((XCDF(nn,mm)>0) && (XCDF(nn,mm)<1))
Zx(nn,mm)=norminv(XCDF(nn,mm));
elseif (XCDF(nn,mm)<.000000001)
Zx(nn,mm)=-5;
elseif (XCDF(nn,mm)>.999999999)
Zx(nn,mm)=5;
end
if (Zx(nn,mm)<-5)
Zx(nn,mm)=-5;
end
if (Zx(nn,mm)>5)
Zx(nn,mm)=5;
end
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
YCDF(1:IndexMaxY)=0;
Zy(1:IndexMaxY)=0;
YCDF(1)=YProbMass(1)*.5;
Zy(1)=norminv(YCDF(1));
for nn=2:IndexMaxY
YCDF(nn)=YCDF(nn-1)+YProbMass(nn-1)*.5+YProbMass(nn)*.5;
Zy(nn)=norminv(YCDF(nn));
if (Zy(nn)<-5)
Zy(nn)=-5;
end
if (Zy(nn)>5)
Zy(nn)=5;
end
end
end
```

.

.

.

```
function [pXY2D] = Calculate2VariablesOn2DGridFrom2DZSeries(XX,YY,NoOfBinsX,NoOfBinsY,C,a0,a,BinSizeX,BinSizeY)
%YY are Y-coordinates of 2D Grid. In our case this is stochastic
%volatility and is univariate and independent of asset so it can be
%inverted as a univariate Z-series. its Z-series is given by a0 and a.
[Zy] = CalculateZgivenXAndZSeriesC7(YY,a0,a);
for mm=1:NoOfBinsY+1
if(Zy(mm)>5.0)
Zy(mm)=5.0;
end
if(Zy(mm)<-5.0)
Zy(mm)=-5.0;
end
end
Zy
str=input("Look at Zy");
bH(1:NoOfBinsY+1,1:7)=0;
bH0(1:NoOfBinsY+1)=0;
for mm=1:NoOfBinsY+1
%bH0(mm)=C(1,1);
for Ovm=1:8
for Oxn=2:8
bH(mm,Oxn-1)=bH(mm,Oxn-1)+C(Oxn,Ovm).*HermiteH(Ovm-1,Zy(mm));
end
bH0(mm)=bH0(mm)+C(1,Ovm).*HermiteH(Ovm-1,Zy(mm));
end
end
for mm=1:NoOfBinsY+1
bH0(mm)=bH0(mm)./YY(mm);
for Oxn=1:7
bH(mm,Oxn)=bH(mm,Oxn)./YY(mm);
end
end
for mm=1:NoOfBinsY+1
[b0,b]=ConvertHCoeffsToZCoeffs(bH0(mm),bH(mm,1:7),7);
[Zx0] = CalculateZgivenXAndZSeriesC7(XX,b0,b);
Zx(1:NoOfBinsX+1,mm)=Zx0(1:NoOfBinsX+1);
end
for mm=1:NoOfBinsY+1
for nn=1:NoOfBinsX+1
if(Zx(nn,mm)>5.0)
Zx(nn,mm)=5.0; %possible error block
end
if(Zx(nn,mm)<-5.0)
Zx(nn,mm)=-5.0; %possible error block
end
% if(Zx(nn+1,mm)<Zx(nn,mm))
% Zx(nn+1,mm)=Zx(nn,mm); %possible error block
% end
end
end
% Hermite2DSeries(1:NoOfBinsX+1,1:NoOfBinsY+1)=0;
% %
% for nn=1:NoOfBinsX+1
% for mm=1:NoOfBinsY+1
% for Oxn=1:8
% for Ovm=1:8
% Hermite2DSeries(nn,mm)=Hermite2DSeries(nn,mm)+C(Oxn,Ovm).*HermiteH(Oxn-1,Zx(nn,mm)).*HermiteH(Ovm-1,Zy(mm));
% end
% end
% end
% end
%Below we make calculations of density in original coordinates by doing
%change of probability derivative with respect to gaussian density.
DinvfXX(1:NoOfBinsX+1,1:NoOfBinsY+1)=0;
for mm=1:NoOfBinsY+1
for nn=2:NoOfBinsX
DinvfXX(nn,mm) = (Zx(nn + 1,mm) - Zx(nn - 1,mm))/(XX(nn + 1) - XX(nn - 1));
end
DinvfXX(NoOfBinsX+1,mm)=DinvfXX(NoOfBinsX,mm);
DinvfXX(1,mm)=DinvfXX(2,mm);
end
DinvfYY(1:NoOfBinsY+1)=0;
for mm=2:NoOfBinsY
DinvfYY(mm) = (Zy(mm+1) - Zy(mm-1))/(YY(mm + 1) - YY(mm - 1));
end
DinvfYY(NoOfBinsY+1)=DinvfYY(NoOfBinsY);
DinvfYY(1)=DinvfYY(2);
pXY2D(1:NoOfBinsX+1,1:NoOfBinsY+1)=0;
for mm=2:NoOfBinsY
for nn=2:NoOfBinsX
pXY2D(nn,mm) = (normpdf(Zx(nn,mm),0, 1))*abs(DinvfXX(nn,mm)).*(normpdf(Zy(mm),0, 1))*abs(DinvfYY(mm));
end
end
end
```

.

.

.

```
function [Z] = CalculateZgivenXAndZSeriesC7(X,c0,c)
Z(1:length(X))=0;
for nn=1:length(X)
%X(nn)
%c0
Z(nn)=CalculateZgivenXAndZSeriesBisection2C7(X(nn),c0,c);
%Z(nn)
end
end
```

.

.

.

```
function [Z]=CalculateZgivenXAndZSeriesBisection2C7(Y,c0,c)
%This program takes a value of Y as input and returns the corresponding
%value of Z given the parameters of the Z-series.
%It first see if the value of Y is greater than median or smaller than
%median. If the value of Y is greater than median Z has to be greater than
%zero and if the value of Y is smaller than median Z has to be negative.
% Once we know whether Z is greater than Z or smaller than Zero, we divide
% the program into two repeated if loops, one for Z greater than zero and
%one for Z smaller than zero.
%If Z is positive, we have to assign it a valid band of upper and
%lower values between which it can be found.
%This is done sequentially by calculating Y(.5) at Z=.5, if Y is smaller
%than Y(.5), we know that true value of Y will lie between 0 and .5
%so we assign a lower band value Za=0 and upper band value Zb=.5
%If however Y is greater than Y(.5), we know that true Z should be larger
%than .5 and we again evaluate Y(1) at Z=1, If Y is smaller than Y(1),
%we know that true value of Z should be between .5 and 1, and we assign a
%lower band value Za=.5 and upper band value of Zb=1. Through repeated
%elseif we continue until we find an appropriate band value of Z.
%Once we have determined the upper band value of Zb and lower band value of
%Za, We find Zmid=(Za+Zb)/2. We then evaluate Y(Zmid) and if Y < Y(Zmid) we
%move the upper band to Zb=Zmid while keeping Za at older value. However if
% Y> Y(Zmid), we assign the lower band value equal to Zmid as Za=Zmid. In a
% few iterations, this bisection procedure gets extremely close to true
% value of Z for which Y(Z)=Y
if(Y>c0)
%Z=.5;
if(Y< EvaluateZSeriesC7(c0,c,.5))
Za=0;
Zb=.5;
elseif(Y< EvaluateZSeriesC7(c0,c,1))
Za=.5;
Zb=1;
elseif(Y< EvaluateZSeriesC7(c0,c,1.5))
Za=1;
Zb=1.5;
elseif(Y< EvaluateZSeriesC7(c0,c,2))
Za=1.5;
Zb=2;
elseif(Y< EvaluateZSeriesC7(c0,c,2.5))
Za=2;
Zb=2.5;
elseif(Y< EvaluateZSeriesC7(c0,c,3))
Za=2.5;
Zb=3;
elseif(Y< EvaluateZSeriesC7(c0,c,3.5))
Za=3;
Zb=3.5;
elseif(Y< EvaluateZSeriesC7(c0,c,4))
Za=3.5;
Zb=4;
elseif(Y< EvaluateZSeriesC7(c0,c,4.5))
Za=4;
Zb=4.5;
elseif(Y< EvaluateZSeriesC7(c0,c,5))
Za=4.5;
Zb=5;
elseif(Y< EvaluateZSeriesC7(c0,c,5.5))
Za=5;
Zb=5.5;
elseif(Y< EvaluateZSeriesC7(c0,c,6))
Za=5.5;
Zb=6;
else
Za=6;
Zb=12;
end
for mm=1:35
Zmid=.5*Za+.5*Zb;
if(Y< EvaluateZSeriesC7(c0,c,Zmid))
Zb=Zmid;
elseif(Y> EvaluateZSeriesC7(c0,c,Zmid))
Za=Zmid;
end
Z=.5*Za+.5*Zb; %This can be taken out of the loop and does not
%necessarily have to be calculated again and again. Has to be
%calculated at least once after the loop to determine the true value of
%Z.
end
end
if(Y<c0)
if(Y> EvaluateZSeriesC7(c0,c,-.5))
Za=-.5;
Zb=0;
elseif(Y> EvaluateZSeriesC7(c0,c,-1))
Za=-1;
Zb=-.5;
elseif(Y> EvaluateZSeriesC7(c0,c,-1.5))
Za=-1.5;
Zb=-1;
elseif(Y> EvaluateZSeriesC7(c0,c,-2))
Za=-2;
Zb=-1.5;
elseif(Y> EvaluateZSeriesC7(c0,c,-2.5))
Za=-2.5;
Zb=-2.0;
elseif(Y> EvaluateZSeriesC7(c0,c,-3))
Za=-3;
Zb=-2.5;
elseif(Y> EvaluateZSeriesC7(c0,c,-3.5))
Za=-3.5;
Zb=-3.0;
elseif(Y> EvaluateZSeriesC7(c0,c,-4))
Za=-4;
Zb=-3.5;
elseif(Y> EvaluateZSeriesC7(c0,c,-4.5))
Za=-4.5;
Zb=-4.0;
elseif(Y> EvaluateZSeriesC7(c0,c,-5))
Za=-4.5;
Zb=-5.0;
elseif(Y> EvaluateZSeriesC7(c0,c,-5.5))
Za=-5.0;
Zb=-5.5;
elseif(Y> EvaluateZSeriesC7(c0,c,-6.0))
Za=-5.5;
Zb=-6.0;
else
Za=-12.0;
Zb=-6.0;
end
for mm=1:35
Zmid=.5*Za+.5*Zb;
if(Y< EvaluateZSeriesC7(c0,c,Zmid))
Zb=Zmid;
elseif(Y> EvaluateZSeriesC7(c0,c,Zmid))
Za=Zmid;
end
Z=.5*Za+.5*Zb; %This can be taken out of the loop and does not
%necessarily have to be calculated again and again. Has to be
%calculated at least once after the loop to determine the true value of
%Z
end
end
end
```

.

.

.

```
function [X] = EvaluateZSeriesC7(c0,c,Z)
X=c0+c(1)*Z+c(2)*Z.^2+c(3)*Z.^3+c(4)*Z.^4+c(5)*Z.^5+c(6)*Z.^6+c(7)*Z.^7;
end
```

.

.

.

This is the last few outputs you should see on your screen when you run this program.

a0 =

0.051856848795143

a =

Columns 1 through 3

0.025308447061333 0.006232951335150 0.001147500450485

Columns 4 through 6

0.000877636380746 0.000394818257614 -0.000020803817895

Column 7

-0.000017810888272

Look at volatility coefficients ----a0,a(1:7)

bH0 =

0.999230567144020

bH =

Columns 1 through 3

0.608033135481096 0.159292825114088 0.027749297555118

Columns 4 through 6

0.002484662411384 -0.001193699130006 -0.000894232267067

Column 7

-0.000332249594428

b0 =

0.860805213270085

b =

Columns 1 through 3

0.541765963280562 0.104144398627779 0.004800081440268

Columns 4 through 6

0.015898146417386 0.005783542352976 -0.000894232267067

Column 7

-0.000332249594428

Look at 1D Coeffs of X

str =

[]

Zy =

Columns 1 through 3

-2.634004451734654 -2.167319894339016 -1.506064279914426

Columns 4 through 6

-0.874060740963614 -0.403909689433931 -0.028130067155871

Columns 7 through 9

0.288271345467365 0.559033179619291 0.791351162428327

Columns 10 through 12

0.991148471650376 1.164052401240042 1.315111681913550

Columns 13 through 15

1.448541791520256 1.567709385250055 1.675248979525350

Columns 16 through 18

1.773211661733512 1.863199918814644 1.946475835640740

Columns 19 through 21

2.024043532168434 2.096710658392112 2.165133997455996

Columns 22 through 24

2.229853354067018 2.291316897368233 2.349900274937681

Columns 25 through 27

2.405921167512133 2.459650482451252 2.511321048172249

Columns 28 through 30

2.561134433439292 2.609266346429649 2.655870948154188

Columns 31 through 33

2.701084328444267 2.745027330245648 2.787807863038324

Columns 34 through 36

2.829522812062351 2.870259626048210 2.910097647218208

Columns 37 through 39

2.949109233573836 2.987360712962982 3.024913200300944

Columns 40 through 42

3.061823303280107 3.098143736897327 3.133923863621021

Columns 43 through 45

3.169210173226020 3.204046713879507 3.238475484708033

Columns 46 through 48

3.272536798591318 3.306269623215485 3.339711907850869

Columns 49 through 51

3.372900903174013 3.405873481293384 3.438666464011476

Columns 52 through 54

3.471316967836174 3.503862775753078 3.536342747749586

Columns 55 through 57

3.568797284948232 3.601268866150349 3.633802681411908

Columns 58 through 60

3.666447395451542 3.699256085390516 3.732287414874008

Columns 61 through 63

3.765607132278092 3.799290020535409 3.833422487899952

Columns 64 through 66

3.868106089015782 3.903462431546359 3.939640211632650

Columns 67 through 69

3.976825640718744 4.015258517880284 4.055258214437345

Columns 70 through 72

4.097268263758451 4.141938997338002 4.190297468179779

Columns 73 through 75

4.244153386120161 4.307333394208399 4.391834105372254

Column 76

5.000000000000000

Look at Zy

Look at graph of Analytic 2D density

Look at difference of Analytic and Numerical density graphs

.

.

.

These are some output graphs in no particular order

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