Friends, I am quickly uploading my programs for data generation by monte carlo and calculation of multidimensional hermite series. These are basic programs and I have still not done graphing for which I will post programs tomorrow.

I have run this program on three mean reverting assets and one SV process.

Description of main function followed by programs

This functions takes as arguments joint observations of NDim correlated

random variables. And returns a multidimensional hermite series expanded to Order-th polynomial. Order in input arguments means order of highest hermite polynomials in hermite series of each of the random variables.
1st variable has a one dimensional hermite series in one orthogonal normal r

andom variable.
2nd variable has a two dimensional hermite series in first two orthogonal normal random variables.
While nth variable has an n dimensional hermite series in terms of n orthogonal standard normal random variables.
The program returns a three dimensional array of hermite series coeffs in

variable bCoeffHH.
bCoeffHH is defined as

bCoeffHH(1:Order+1,1:NDim,1:NDium);
The first dimension in bCoeffHH is associated with hermite series expanded till

Order-th hermite polynomial;
The second dimension in bCoeffHH indicates mth orthogonal variable.

The third dimension is associated with the number of random variable.

so

bCoeffHH(:,3,4) indicates hermite series associated with third orthogonal Z of the 4th random variable.
Hermite Series Form used in this program is given as
Order in hermite series below is five as it is expanded to fifth hermite

polynomial.
n also varies from 1 to NDim. If NDim=4, there are four orthogonal Z_n as

in notation example below.
Y(n) = bnm(1,1,n) + bCoeffHH(2,1,n) H_1(Z_1) + bCoeffHH(3,1,n) H_2(Z_1) + bCoeffHH(4,1,n) H_3(Z_1) + bCoeffHH(5,1,n) H_4(Z_1) + bCoeffHH(6,1,n) H_5(Z_1)

+ bCoeffHH(1,2,n) + bCoeffHH(2,2,n) H_1(Z_2) + bCoeffHH(3,2,n) H_2(Z_2) + bCoeffHH(4,2,n) H_3(Z_2) + bCoeffHH(5,2,n) H_4(Z_2) + bCoeffHH(6,2,n) H_5(Z_2)

+ bCoeffHH(1,3,n) + bCoeffHH(2,3,n) H_1(Z_3) + bCoeffHH(3,3,n) H_2(Z_3) + bCoeffHH(4,3,n) H_3(Z_3) + bCoeffHH(5,3,n) H_4(Z_3) + bCoeffHH(6,3) H_5(Z_3)

+ bCoeffHH(1,4,n) + bCoeffHH(2,4,n) H_1(Z_4) + bCoeffHH(3,4,n) H_2(Z_4) + bCoeffHH(4,4,n) H_3(Z_4) + bCoeffHH(5,4,n) H_4(Z_4) + bCoeffHH(6,4,n) H_5(Z_4)

So the above hermite series representation will be repeated for all of the NDim variables. But 1st variable will have its representation in only one orthogonal hermoite and rest of the columns would be zero while 2nd variable will have its representation in two orhtogonal variables and rest of the columns would be zero and so on.

...

.

.

.

```
function []= SVMonteCarloAndAssetConditionalDensityHermitesMultiDimensional()
%In this program, we jointly simulate the system of SV and asset SDEs to
%get the joint data of various mean reverting assets and stochastic volatility.
%Asset SDE is
%dX(t) = a X(t)^alpha1 dt+ b X(t)^alpha2 dt+ sigma1 V(t)^gammaV X(t)^gammaX dZ(t)
%We usually take a=0 and b=0 in our simulations
%SV SDE is
%dV(t)=mu1 V(t)^beta1 dt+ mu2 V(t)^beta2 dt + sigma0 V(t)^gamma dZ(t)
%In mean reverting SDEs we ususally have
%mu1= kappaV * thetaV
%beta1=0
%mu2=-kappa
%beta2=0
NDim=4;%Three assets and one SV.
X10=1.0;
alpha1(1)=0;
alpha2(1)=1;
kappaX(1)=.5;
thetaX(1)=1;
a(1)=kappaX(1)*thetaX(1);
b(1)=-kappaX(1);
rho(1)=-.8;
%sigma1(nn)=1;
sigma1(1)=.25;%.225;%1;
gammaV=.5;
gammaX(1)=.95;
% %%%%%%%%`%%%%%
X20=1.0;
alpha1(2)=0;
alpha2(2)=1;
kappaX(2)=1;
thetaX(2)=1;
a(2)=kappaX(2)*thetaX(2);
b(2)=-kappaX(2);
rho(2)=-.7;
%sigma1(2)=1;
sigma1(2)=.25;%.225;%1;
gammaV=.5;
gammaX(2)=.95;
%%%%%%%%%%%%%%%%%%%%%%5
X30=1.0;
alpha1(3)=0;
alpha2(3)=1;
kappaX(3)=1.5;
thetaX(3)=1;
a(3)=kappaX(3)*thetaX(3);
b(3)=-kappaX(3);
rho(3)=-.5;
%sigma1(3)=1;
sigma1(3)=.25;%.225;%1;
gammaV=.5;
gammaX(3)=.95;
%%%%%%%%%%%%%%%%%%%5
V0=1;%1;%.35;
thetaV=1;%1;%.04;
kappaV=1.5;%1.5;
mu1=kappaV*thetaV;
mu2=-kappaV;
beta1=0;
beta2=1;
gamma=.5;
sigma0=.95%1.250%.5;
%V0=1;
%thetaV=1;
%kappaV=1.50;
mu1=kappaV*thetaV;
mu2=-kappaV;
beta1=0;
beta2=1;
%gamma=.5;
%sigma0=1.25;
dt=.03125/2;
Tt=64*3;
T=Tt*dt;
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%In this code block, we do monte carlo simulation of stochastic volatility SDEs
%And also calculate monte carlo call prices.
seed0=21872418;
seed0=82150987;
rng(seed0, 'twister')
paths=1000;
V(1:paths,1)=V0;
X1(1:paths,1)=X10;
X1(1:paths,2)=X20;
X1(1:paths,3)=X30;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
paths0=paths;
paths00=paths;
Random1(1:paths00,1:NDim-1)=0;
Random2(1:paths0,1)=0;
for tt=1:Tt
Random2=randn(size(Random2));
Random1=randn(size(Random1));
for nn=1:NDim-1
X1(1:paths,nn)=X1(1:paths,nn)+ ...
(a(nn).* X1(1:paths,nn).^alpha1(nn) + b(nn).* X1(1:paths,nn).^alpha2(nn))* dt + ...
sqrt(1-rho(nn)^2).* sigma1(nn).* V(1:paths,1).^gammaV.* X1(1:paths,nn).^gammaX(nn) .*Random1(1:paths,nn) * sqrt(dt) + ...
rho(nn).* sigma1(nn).* V(1:paths,1).^gammaV .*X1(1:paths,nn).^gammaX(nn) .*Random2(1:paths,1)*sqrt(dt) + ...
(a(nn).*alpha1(nn).* X1(1:paths,nn).^(alpha1(nn)-1)+b(nn).*alpha2(nn).* X1(1:paths,nn).^(alpha2(nn)-1)).* ...
(((a(nn).* X1(1:paths,nn).^alpha1(nn) + b(nn).* X1(1:paths,nn).^alpha2(nn)).* dt^2/2)+ ...
(sqrt(1-rho(nn)^2).* sigma1(nn)* V(1:paths,1).^gammaV.* X1(1:paths,nn).^gammaX(nn) .*Random1(1:paths,nn) .*(1-1/sqrt(3)).*dt^1.5+ ...
rho(nn).* sigma1(nn).* V(1:paths,1).^gammaV .*X1(1:paths,nn).^gammaX(nn) .*Random2(1:paths,1).*(1-1/sqrt(3)).*dt^1.5))+ ...
.5*(a(nn).*alpha1(nn).*(alpha1(nn)-1).* X1(1:paths,nn).^(alpha1(nn)-2)+b(nn).*alpha2(nn).*(alpha2(nn)-1).* X1(1:paths,nn).^(alpha2(nn)-2)).* ...
( sigma1(nn)^2.* V(1:paths,1).^(2*gammaV).* X1(1:paths,nn).^(2*gammaX(nn))) .*dt^2/2 + ...
sqrt(1-rho(nn)^2).* sigma1(nn).* V(1:paths,1).^gammaV.*gammaX(nn).* X1(1:paths,nn).^(gammaX(nn)-1).* ...
((a(nn).* X1(1:paths,nn).^alpha1(nn) + b(nn)* X1(1:paths,nn).^alpha2(nn)).*Random1(1:paths,nn) * 1/sqrt(3).* dt^1.5 + ...
sqrt(1-rho(nn)^2).* sigma1(nn).* V(1:paths,1).^gammaV.* X1(1:paths,nn).^gammaX(nn) .*(Random1(1:paths,nn).^2-1) * dt/2 + ...
rho(nn)* sigma1(nn)* V(1:paths,1).^gammaV .*X1(1:paths,nn).^gammaX(nn) .*Random1(1:paths,nn).*Random2(1:paths,1)*dt/2)+ ...
.5*sqrt(1-rho(nn)^2)* sigma1(nn)* V(1:paths,1).^gammaV.*gammaX(nn).*(gammaX(nn)-1).* X1(1:paths,nn).^(gammaX(nn)-2).* ...
(sigma1(nn)^2* V(1:paths,1).^(2*gammaV).* X1(1:paths,nn).^(2*gammaX(nn)) .*Random1(1:paths,nn).*1/sqrt(3).*dt^1.5)+ ...
sqrt(1-rho(nn)^2)* sigma1(nn)*gammaV.* V(1:paths,1).^(gammaV-1).* X1(1:paths,nn).^(gammaX(nn)).* ...
((mu1.*V(1:paths,1).^beta1 + mu2.*V(1:paths,1).^beta2).*Random1(1:paths,nn) * 1/sqrt(3).* dt^1.5 + ...
sigma0*V(1:paths,1).^gamma.*Random1(1:paths,nn).*Random2(1:paths,1)*dt/2)+ ...
.5*sqrt(1-rho(nn)^2)* sigma1(nn)*gammaV.*(gammaV-1).* V(1:paths,1).^(gammaV-2).* X1(1:paths,nn).^(gammaX(nn)).* ...
(sigma0^2*V(1:paths,1).^(2*gamma).*Random1(1:paths,nn)*1/sqrt(3)*dt^1.5)+ ...
sqrt(1-rho(nn)^2)* sigma1(nn)*gammaV.* V(1:paths,1).^(gammaV-1).*gammaX(nn).* X1(1:paths,nn).^(gammaX(nn)-1).* ...
rho(nn).* sigma1(nn).* V(1:paths,1).^gammaV .*X1(1:paths,nn).^gammaX(nn) .*sigma0.*V(1:paths,1).^gamma.*Random1(1:paths,nn)*1/sqrt(3)*dt^1.5+ ...
rho(nn).* sigma1(nn).* V(1:paths,1).^gammaV.*gammaX(nn).* X1(1:paths,nn).^(gammaX(nn)-1).* ...
((a(nn).* X1(1:paths,nn).^alpha1(nn) + b(nn).* X1(1:paths,nn).^alpha2(nn)).*Random2(1:paths,1) * 1/sqrt(3).* dt^1.5 + ...
sqrt(1-rho(nn)^2).* sigma1(nn).* V(1:paths,1).^gammaV.* X1(1:paths,nn).^gammaX(nn) .*Random1(1:paths,nn).*Random2(1:paths,1) * dt/2 + ...
rho(nn).* sigma1(nn).* V(1:paths,1).^gammaV .*X1(1:paths,nn).^gammaX(nn) .*(Random2(1:paths,1).^2-1)*dt/2)+ ...
.5*rho(nn).* sigma1(nn).* V(1:paths,1).^gammaV.*gammaX(nn).*(gammaX(nn)-1).* X1(1:paths,nn).^(gammaX(nn)-2).* ...
(sigma1(nn).^2* V(1:paths,1).^(2*gammaV).* X1(1:paths,nn).^(2*gammaX(nn)) .*Random2(1:paths,1).*1/sqrt(3).*dt^1.5)+ ....
rho(nn).* sigma1(nn).*gammaV.* V(1:paths,1).^(gammaV-1).* X1(1:paths,nn).^(gammaX(nn)).* ...
((mu1.*V(1:paths,1).^beta1 + mu2.*V(1:paths,1).^beta2).*Random2(1:paths,1) * 1/sqrt(3).* dt^1.5 + ...
sigma0*V(1:paths,1).^gamma.*(Random2(1:paths,1).^2-1)*dt/2)+ ...
.5*rho(nn).* sigma1(nn).*gammaV.*(gammaV-1).* V(1:paths,1).^(gammaV-2).* X1(1:paths,nn).^(gammaX(nn)).* ...
sigma0^2.*V(1:paths,1).^(2*gamma).*Random2(1:paths,1) * 1/sqrt(3).* dt^1.5+ ...
rho(nn).* sigma1(nn).*gammaV.* V(1:paths,1).^(gammaV-1).*gammaX(nn).* X1(1:paths,nn).^(gammaX(nn)-1).* ...
rho(nn).* sigma1(nn).* V(1:paths,1).^gammaV .*X1(1:paths,nn).^gammaX(nn) .*sigma0.*V(1:paths,1).^gamma.*Random2(1:paths,1)*1/sqrt(3)*dt^1.5;
end
V(1:paths,1)=V(1:paths,1)+ ...
(mu1.*V(1:paths,1).^beta1 + mu2.*V(1:paths,1).^beta2)*dt + ...
sigma0*V(1:paths,1).^gamma .*Random2(1:paths,1)*sqrt(dt) + ...
(mu1.*beta1*V(1:paths,1).^(beta1-1) + mu2.*beta2.*V(1:paths,1).^(beta2-1)).* ...
((mu1.*V(1:paths,1).^beta1 + mu2.*V(1:paths,1).^beta2)*dt^2/2 + ...
sigma0*V(1:paths,1).^gamma .*Random2(1:paths,1)*(1-1/sqrt(3))*dt^1.5) + ...
.5*(mu1.*beta1.*(beta1-1).*V(1:paths,1).^(beta1-2) + mu2.*beta2.*(beta2-1).*V(1:paths,1).^(beta2-2)).* ...
sigma0^2.*V(1:paths,1).^(2*gamma).*dt^2/2 + ...
sigma0*gamma*V(1:paths,1).^(gamma-1) .* ...
((mu1.*V(1:paths,1).^beta1 + mu2.*V(1:paths,1).^beta2).*Random2(1:paths,1).*1/sqrt(3)*dt^1.5 + ...
sigma0.*V(1:paths,1).^gamma .*(Random2(1:paths,1).^2-1)*dt/2) + ...
.5*sigma0*gamma*(gamma-1).*V(1:paths,1).^(gamma-2) .* ...
sigma0^2.*V(1:paths,1).^(2*gamma) .*Random2(1:paths,1).*1/sqrt(3)*dt^1.5;
end
% X1=X;
% V1=V;
%
%
X1(X1<0)=0;
V(V<0)=0;
X1=real(X1);
V=real(V);
Yin(1:paths,1)=V(1:paths,1);
Yin(1:paths,2:4)=X1(1:paths,1:3);
NDim=4;
Order=5;
[bCoeffHH] = CalculateMultivariateHermiteSeriesAndDensities01(Yin,NDim,Order)
end
```

.

.

.

```
function [bCoeffHH] = CalculateMultivariateHermiteSeriesAndDensities01(Yin,NDim,Order)
%This functions takes as arguments joint observations of NDim correlated
%random variables.
%And returns a multidimensional hermite series expanded to Order-th
%polynomial. Order in input arguments means order of highest hermite
%polynomials in hermite series of each of the random variables.
%1st variable has a one dimensional hermite series in one orthogonal normal
%random variable. 2nd variable has a two dimensional hermite series in
%first two orthogonal normal random variables. While nth variable has an n
%dimensional hermite series in terms of n orthogonal standard normal random
%variables.%
%The program returns a three dimensional array of hermite series coeffs in
%variable bCoeffHH.
%bCoeffHH is defined as
%bCoeffHH(1:Order+1,1:NDim,1:NDium);
%The first dimension in bCoeffHH is associated with hermite series expanded till
%Order-th hermite polynomial;
%The second dimension in bCoeffHH indicates mth orthogonal variable.
%The third dimension is associated with the number of random variable. so
%bCoeffHH(:,3,4) indicates hermite series associated with third orthogonal
%Z of the 4th random variable.
%We will first standardize the data of all random variables and then
%calculate univariate hermite series of all variables.
% Hermite Series Form used in this program is given as
% Order in hermite series below is five as it is expanded to sixth hermite
% polynomial.
% n also varies from 1 to NDim. If NDim=4, there are four orthogonal Z_n as
% in notation example below.
% Y(n) = bnm(1,1,n) + bCoeffHH(2,1,n) H_1(Z_1) + bCoeffHH(3,1,n) H_2(Z_1) + bCoeffHH(4,1,n) H_3(Z_1) + bCoeffHH(5,1,n) H_4(Z_1) + bCoeffHH(6,1,n) H_5(Z_1)
% + bCoeffHH(1,2,n) + bCoeffHH(2,2,n) H_1(Z_2) + bCoeffHH(3,2,n) H_2(Z_2) + bCoeffHH(4,2,n) H_3(Z_2) + bCoeffHH(5,2,n) H_4(Z_2) + bCoeffHH(6,2,n) H_5(Z_2)
% + bCoeffHH(1,3,n) + bCoeffHH(2,3,n) H_1(Z_3) + bCoeffHH(3,3,n) H_2(Z_3) + bCoeffHH(4,3,n) H_3(Z_3) + bCoeffHH(5,3,n) H_4(Z_3) + bCoeffHH(6,3) H_5(Z_3)
% + bCoeffHH(1,4,n) + bCoeffHH(2,4,n) H_1(Z_4) + bCoeffHH(3,4,n) H_2(Z_4) + bCoeffHH(4,4,n) H_3(Z_4) + bCoeffHH(5,4,n) H_4(Z_4) + bCoeffHH(6,4) H_5(Z_4)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%In the code block below, we standardize the independent variable Y and dependent variable X.
%We do all later calculations of various Z-series and cross moments in
%standardized framework and then at the end of the program, we convert 00
%the calculated coefficients of 2D hermite series/Z-series of dependent
%variable X into original non-standardized framework.
Ndata=length(Yin);
Hz=HermitePoly(Order);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for nn=1:NDim
MeanYin(nn)=sum(Yin(1:Ndata,nn))/Ndata;
Ym(1:Ndata,nn)=Yin(1:Ndata,nn)-MeanYin(nn);
YmVar(nn)=sum(Ym(1:Ndata,nn).^2)/Ndata;
Y(1:Ndata,nn)=Ym(1:Ndata,nn)/sqrt(YmVar(nn));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Code block below uses polynomial regression to find univariate Z-series
%coefficients of independent variable Y.
%To know more about polynomial regression read the following and associated
%posts: https://forum.wilmott.com/viewtopic.php?t=99702&start=2010#p877143
%You can get the text file MomentMatchParams08.txt from the forum post below
%https://forum.wilmott.com/viewtopic.php?f=4&t=99702&start=2040#p877863
Mtable=readtable('C:\Users\Lenovo\Documents\MATLAB\MATLAB1\MomentMatchParams08.txt');
%ZI is moment matched grid along the normal density with number of grid
%points equal to data elements in Y given by Ndata. The grid is constructed
%so that probability mass in each grid cell is 1/Ndata so this matches the
%observation probability of Ndata data elements. Due to equal probability
%mass in each grid cell, we call it equiprobable grid.
[ZI] = MatchMomentParametersOfIdealZhalf04(Ndata,Mtable);
[W] = EvalHermiteArrayH0(ZI,Order); %We will continue to use the same
%design matrix W many times.
%Below, we calculate 1 Dim Hermite Series of each variable.
for nn=1:NDim
Ys=sort(Y(1:Ndata,nn));
Ymu(1:Ndata,1)=Ys(1:Ndata);
coeff=inv(W'*W)*(W'*Ymu);
Coeff(1:Order+1,nn)= coeff(1:Order+1,1);
Z0= CalculateZgivenXAndHSeriesCoeffs(Y(1:Ndata,nn),Coeff(1:Order+1,nn),Order);
Zy(1:Ndata,nn)=Z0(1:Ndata);
end
bCoeffHH(1:Order+1,1,1)=Coeff(1:Order+1,1);
bCoeffHH(1,1,1:NDim)=Coeff(1,1:NDim);
ZOrtho(1:Ndata,1)=Zy(1:Ndata,1);
%In the nth iteration of the loop below, we calculate hermite correlations
%of nth random varaible Y(n) with all of
%the previous n-1 orthogonal Z's, Z_1, Z_2, ..,Z_{n-1}. We use all these
%hermite correlations with previous orthogonal Z's called ZOrtho in the
%program and use them to find value of Y(n) that is correlated with
%previous n-1 orthogonal Z's. We remove this correlation to find
%decorrelated version of Y(n) and then regress this decorrelated
%version to find Coefficients of decorrelated Y(n).
for nn=2:NDim
Ycorr(1:Ndata,nn)=0;
for mm=1:nn-1
Z1(:)=ZOrtho(:,mm);
Z2(:)=Zy(:,nn);
[Wh] = EvalHermiteArrayH0(Z1,Order);
[CorrH] = CalculateCorrelationBivariateHermiteCH0(Z1,Z2,Ndata,Order);
for hh=1:Order
Ycorr(1:Ndata,nn)=Ycorr(1:Ndata,nn)+Coeff(hh+1,nn).*CorrH(hh).*Wh(1:Ndata,hh+1);
bCoeffHH(hh+1,mm,nn)=Coeff(hh+1,nn).*CorrH(hh);
end
end
Ydecorr(1:Ndata,nn)=Y(1:Ndata,nn)-Ycorr(1:Ndata,nn);
Ys=sort(Ydecorr(1:Ndata,nn));
Ymu(1:Ndata,1)=Ys(1:Ndata);
coeff=inv(W'*W)*(W'*Ymu);
Z0= CalculateZgivenXAndHSeriesCoeffs(Ydecorr(1:Ndata,nn),coeff(1:Order+1,1),Order);
ZOrtho(1:Ndata,nn)=Z0(1:Ndata);
bCoeffHH(1:Order+1,nn,nn)=coeff(1:Order+1,1);
end
bCoeffHH
str=input("Look at bCoeffHH")
for nn=1:NDim
[TVar(nn),Var1] = CalculateHermiteSeriesVar(bCoeffHH(:,:,nn),Order,NDim);
TVar(nn)
Var1
end
str=input("Look at variances opposed to normalized version")
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Reverse the standardization of each variable.
for nn=1:NDim
bCoeffHH(1:Order+1,:,nn)=bCoeffHH(1:Order+1,:,nn).*sqrt(YmVar(nn));
end
end
```

.

.

.

```
function [Hz] = HermitePoly(Order)
Hz(1:Order+1,1:Order+1)=0.0;
Hz(1,1)=1.0;
for k=1:Order
for m=1:k
Hz(k+1,m+1)=Hz(k,m);
if (m>1)
Hz(k+1,m-1)=Hz(k+1,m-1)-(m-1)*Hz(k,m);
end
end
end
end
```

.

.

.

```
function [ZnNew] = MatchMomentParametersOfIdealZhalf04(Ndata,Mtable)
[Zn] = StandardNormalIdealDesity02(Ndata,1);
[alpha0,beta0,gamma0] = ReadParamsFile(Ndata,Mtable);
ZnNew(1:Ndata)=Zn(1:Ndata).*(alpha0+beta0*(Zn(1:Ndata).^2-1)+gamma0*(Zn(1:Ndata).^4-6*Zn(1:Ndata).^2+3));
Moment2=sum(ZnNew(1:Ndata).^2)*1/Ndata;
Moment4=sum(ZnNew(1:Ndata).^4)*1/Ndata;
Moment6=sum(ZnNew(1:Ndata).^6)*1/Ndata;
%Moment2
%Moment4
%Moment6
%display("fitted Moments");
%str=input("Look at moment values")
end
```

.

.

.

```
function [alpha0,beta0,gamma0] = ReadParamsFile(Ndata,Mtable)
A=Mtable;
NSubDivs=A.(1);
alpha=A.(2);
beta=A.(3);
gamma=A.(4);
% NSubDivs
% alpha
% beta
% gamma
%
% A.(1)(1)
%Ndata;
nn=1;
if(Ndata<=NSubDivs(length(NSubDivs)))
while (NSubDivs(nn)<=Ndata)
if(NSubDivs(nn)==Ndata)
alpha0=alpha(nn);
beta0=beta(nn);
gamma0=gamma(nn);
NSubDivs(nn)
display("I am in part--1");
elseif((NSubDivs(nn)<Ndata)&&(NSubDivs(nn+1)>Ndata))
alpha1=alpha(nn-1);
alpha2=alpha(nn);
alpha3=alpha(nn+1)
alpha4=alpha(nn+2);
beta1=beta(nn-1);
beta2=beta(nn);
beta3=beta(nn+1)
beta4=beta(nn+2);
gamma1=gamma(nn-1);
gamma2=gamma(nn);
gamma3=gamma(nn+1)
gamma4=gamma(nn+2);
x1=NSubDivs(nn-1);
x2=NSubDivs(nn);
x3=NSubDivs(nn+1);
x4=NSubDivs(nn+2);
x0=Ndata;
N=4;
[alpha0] = InterpolateOrderN6(N,x0,x1,x2,x3,x4,0,0,alpha1,alpha2,alpha3,alpha4,0,0);
[beta0] = InterpolateOrderN6(N,x0,x1,x2,x3,x4,0,0,beta1,beta2,beta3,beta4,0,0);
[gamma0] = InterpolateOrderN6(N,x0,x1,x2,x3,x4,0,0,gamma1,gamma2,gamma3,gamma4,0,0);
%alpha0=alpha(nn)*(Ndata-NSubDivs(nn))/(NSubDivs(nn+1)-NSubDivs(nn))+alpha(nn+1)*(NSubDivs(nn+1)-Ndata)/(NSubDivs(nn+1)-NSubDivs(nn));
%beta0=beta(nn)*(Ndata-NSubDivs(nn))/(NSubDivs(nn+1)-NSubDivs(nn))+beta(nn+1)*(NSubDivs(nn+1)-Ndata)/(NSubDivs(nn+1)-NSubDivs(nn));
%gamma0=gamma(nn)*(Ndata-NSubDivs(nn))/(NSubDivs(nn+1)-NSubDivs(nn))+gamma(nn+1)*(NSubDivs(nn+1)-Ndata)/(NSubDivs(nn+1)-NSubDivs(nn));
display("I am in part--2");
end
nn=nn+1;
end
else
display("data out of range")
end
end
```

.

.

.

```
function [y0] = InterpolateOrderN6(N,x0,x1,x2,x3,x4,x5,x6,y1,y2,y3,y4,y5,y6)
if(N==6)
y0=(x0-x2).*(x0-x3).*(x0-x4).*(x0-x5).*(x0-x6)./((x1-x2).*(x1-x3).*(x1-x4).*(x1-x5).*(x1-x6)).*y1+ ...
(x0-x1).*(x0-x3).*(x0-x4).*(x0-x5).*(x0-x6)./((x2-x1).*(x2-x3).*(x2-x4).*(x2-x5).*(x2-x6)).*y2+ ...
(x0-x1).*(x0-x2).*(x0-x4).*(x0-x5).*(x0-x6)./((x3-x1).*(x3-x2).*(x3-x4).*(x3-x5).*(x3-x6)).*y3+ ...
(x0-x1).*(x0-x2).*(x0-x3).*(x0-x5).*(x0-x6)./((x4-x1).*(x4-x2).*(x4-x3).*(x4-x5).*(x4-x6)).*y4+ ...
(x0-x1).*(x0-x2).*(x0-x3).*(x0-x4).*(x0-x6)./((x5-x1).*(x5-x2).*(x5-x3).*(x5-x4).*(x5-x6)).*y5+ ...
(x0-x1).*(x0-x2).*(x0-x3).*(x0-x4).*(x0-x5)./((x6-x1).*(x6-x2).*(x6-x3).*(x6-x4).*(x6-x5)).*y6;
end
if(N==5)
y0=(x0-x2).*(x0-x3).*(x0-x4).*(x0-x5)./((x1-x2).*(x1-x3).*(x1-x4).*(x1-x5)).*y1+ ...
(x0-x1).*(x0-x3).*(x0-x4).*(x0-x5)./((x2-x1).*(x2-x3).*(x2-x4).*(x2-x5)).*y2+ ...
(x0-x1).*(x0-x2).*(x0-x4).*(x0-x5)./((x3-x1).*(x3-x2).*(x3-x4).*(x3-x5)).*y3+ ...
(x0-x1).*(x0-x2).*(x0-x3).*(x0-x5)./((x4-x1).*(x4-x2).*(x4-x3).*(x4-x5)).*y4+ ...
(x0-x1).*(x0-x2).*(x0-x3).*(x0-x4)./((x5-x1).*(x5-x2).*(x5-x3).*(x5-x4)).*y5;
end
if(N==4)
y0=(x0-x2).*(x0-x3).*(x0-x4)./((x1-x2).*(x1-x3).*(x1-x4)).*y1+ ...
(x0-x1).*(x0-x3).*(x0-x4)./((x2-x1).*(x2-x3).*(x2-x4)).*y2+ ...
(x0-x1).*(x0-x2).*(x0-x4)./((x3-x1).*(x3-x2).*(x3-x4)).*y3+ ...
(x0-x1).*(x0-x2).*(x0-x3)./((x4-x1).*(x4-x2).*(x4-x3)).*y4;
end
if(N==3)
y0=(x0-x2).*(x0-x3)./((x1-x2).*(x1-x3)).*y1+ ...
(x0-x1).*(x0-x3)./((x2-x1).*(x2-x3)).*y2+ ...
(x0-x1).*(x0-x2)./((x3-x1).*(x3-x2)).*y3;
end
if(N==2)
y0=(x0-x2)./((x1-x2)).*y1+ ...
(x0-x1)./((x2-x1)).*y2;
end
%y0=(x0-x3)/(x2-x3)*y2+(x0-x2)/(x3-x2)*y3;
end
```

.

.

.

```
function [Wh] = EvalHermiteArrayH0(Z,Order)
%Evaluates first Order-th hermites at the array of Z's. If input array of Z
%has a size of 100 and Order is 4, it will return a 100 X 5 array.
Hz=HermitePoly(Order);
Wh(1:length(Z),1:Order+1)=0;
for nn=1:Order+1
Wh(:,nn)=0;
for mm=1:nn
Wh(:,nn)=Wh(:,nn)+Hz(nn,mm).*Z(:).^(mm-1);
end
end
end
```

.

.

.

```
function [Z] = CalculateZgivenXAndHSeriesCoeffs(X,Coeff,Order)
%Finds Z given an array of X-values and Coeffs of hermite series.
Z(1:length(X))=0;
Hz=HermitePoly(Order);
median=0;
for nn=1:Order+1
median=median+Coeff(nn).*Hz(nn,1);
end
for nn=1:length(X)
Z(nn)=CalculateZgivenXAndHermiteSeriesBisection(X(nn),Coeff,median,Hz,Order);
end
end
```

.

.

.

```
function [Z]=CalculateZgivenXAndHermiteSeriesBisection(Y,coeff,median,Hz,Order)
%[Z]=CalculateZgivenXAndZSeriesBisection2C5(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
%EvalHermiteSeriesH0(Coeff,Z,Order,Hz)
Z=(Y-coeff(1))/coeff(2);
BracketSize=.5;
if(Y>=median)
Ztry=.5;
ZtryBefore=0;
NotBracketed=1;
while((NotBracketed)&&(Ztry<=8))
[Ytry,dYtrydZ]= EvalHermiteSeriesH0(coeff,Ztry,Order,Hz);
if((Y<= Ytry) || (dYtrydZ < 0))
Za=ZtryBefore;
Zb=Ztry;
NotBracketed=0;
end
ZtryBefore=Ztry;
Ztry=Ztry+BracketSize;
end
if(NotBracketed==1)
Za=8;
Zb=12;
end
for mm=1:20
Zmid=.5*Za+.5*Zb;
[Ymid,dYmiddZ]= EvalHermiteSeriesH0(coeff,Zmid,Order,Hz);
if((Y< Ymid)|| (dYmiddZ<0))
Zb=Zmid;
elseif(Y> Ymid)
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<median)
Ztry=-.5;
ZtryBefore=0;
NotBracketed=1;
while((NotBracketed)&&(Ztry>=-8))
[Ytry,dYtrydZ]= EvalHermiteSeriesH0(coeff,Ztry,Order,Hz);
if((Y>= Ytry) || (dYtrydZ < 0))
Za=Ztry;
Zb=ZtryBefore;
NotBracketed=0;
end
ZtryBefore=Ztry;
Ztry=Ztry-BracketSize;
end
if(NotBracketed==1)
Za=-12;
Zb=-8;
end
for mm=1:20
Zmid=.5*Za+.5*Zb;
[Ymid,dYmiddZ]= EvalHermiteSeriesH0(coeff,Zmid,Order,Hz);
if((Y> Ymid)|| (dYmiddZ<0))
Za=Zmid;
elseif(Y< Ymid)
Zb=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 [CorrH] = CalculateCorrelationBivariateHermiteCH0(Z1,Z2,paths,HH)
CorrH(1:HH)=0;
CovH(1:HH)=0;
Var1H(1:HH)=0;
Var2H(1:HH)=0;
[X1h] = EvalHermiteArrayH0(Z1,HH);
[X2h] = EvalHermiteArrayH0(Z2,HH);
for hh=1:HH
CovH(hh)=CovH(hh)+sum(X1h(:,hh+1).*X2h(:,hh+1))/paths;
Var1H(hh)=Var1H(hh)+sum(X1h(:,hh+1).*X1h(:,hh+1))/paths;
Var2H(hh)=Var2H(hh)+sum(X2h(:,hh+1).*X2h(:,hh+1))/paths;
CorrH(hh)=CovH(hh)./sqrt(Var1H(hh))./sqrt(Var2H(hh));
end
end
```

.

.

.

```
function [Xh,dXh] = EvalHermiteSeriesH0(Coeff,Z,Order,Hz)
%Hz=HermitePoly(Order);
%Xh=Coeff(1);
Xh=0;
dXh=0;
for nn=1:Order+1
for mm=1:nn
Xh=Xh+Coeff(nn).*Hz(nn,mm).*Z(:).^(mm-1);
end
for mm=1:nn-1
dXh=dXh+Coeff(nn).*(nn-1)*Hz(nn-1,mm).*Z(:).^(mm-1);
end
end
end
```

.

.

.

```
function [TVar,Var] = CalculateHermiteSeriesVar(bcoeff,Order,NDim)
for nn=1:NDim
Var(nn)=0;
for mm=1:Order+1
Var(nn)=Var(nn)+bcoeff(mm,nn).^2.*factorial(mm-1);
end
end
TVar=sum(Var(:));
end
```

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