Breakthrough in the theory of stochastic differential equations and their simulation

Amin · May 23rd, 2026, 8:54 pm

.
.
Friends, I tried to write a good modular and re-usable code for the hermite to Gaussian correlation by Newton part of the program but it is too late at the moment so I will post a new fully commented version of the code tomorrow. Its basically very very simple and I hope most friends would like using this code.

.
.

function []= ConditionalDensityHermitesNvariateNewton()

%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
%The present program works only when drift is zero (for pure noises)
%meaning mu1=0 and mu2=0;

Order=7;
%NDim=4;%Three assets and one SV.
%NMomentsY=Order+1;
%NMomentsX=Order+1;

%%%%%%%%%%%%%%%%%%%5
dt=.03125/2/8;
Tt=64*2*8;
T=Tt*dt;
seed0=12130649;
%seed0=94210876;
%seed0=14210876;
%seed0=53612297
%seed0=97459782
rng(seed0, 'twister')
paths=50000;
TT1=100*8/2; %Transition distribution start
TT2=256*8/2; %Transition distribution end

%%%%%%%%%%%%%%%%%%%5

NDims=3;
% 
V0(1)=1.000;%.32;
V00(1)=V0(1);

%SDE1
 thetaV(1)=1.0;%.045;%1;%.04;
 kappaV(1)=0.000;%1.5;%1.5;
 gamma(1)=.75;%.75
 sigma0(1)=.5;%.35


 
%SDE2%
%thetaV(1)=.90;%.045;%1;%.04;
%kappaV(1)=1.000;%1.5;%1.5;
%gamma(1)=.75;%.75
%sigma0(1)=.85;%.35

%SDE3%
%thetaV(1)=1.350;%.045;%1;%.04;
%kappaV(1)=2.000;%1.5;%1.5;
%gamma(1)=.75;%.75
%sigma0(1)=.95;%.35


mu1(1)=kappaV(1)*thetaV(1);
mu2(1)=-kappaV(1);
beta1(1)=0;
beta2(1)=1;



V0(2)=1.000;%.32;
V00(2)=V0(2);

%SDE1
thetaV(2)=1.0;%.045;%1;%.04;
 kappaV(2)=0.000;%1.5;%1.5;
 gamma(2)=.65;%.75
 sigma0(2)=.55;%.35

%SDE2 
% thetaV(2)=1.0;%.045;%1;%.04;
% kappaV(2)=2.000;%1.5;%1.5;
% gamma(2)=.65;%.75
% sigma0(2)=.95;%.35

%SDE3
%thetaV(2)=.850;%.045;%1;%.04;
%kappaV(2)=.35000;%1.5;%1.5;
%gamma(2)=.65;%.75
%sigma0(2)=.55;%.35


mu1(2)=kappaV(2)*thetaV(2);
mu2(2)=-kappaV(2);
beta1(2)=0;
beta2(2)=1;



V0(3)=1.000;%.33;
V00(3)=V0(3);

%SDE1
thetaV(3)=1.0;%.045;%1;%.04;
 kappaV(3)=0.000;%1.5;%1.5;
 gamma(3)=.85;%.75
 sigma0(3)=.45;%.35

%SDE2 
% thetaV(3)=.750;%.045;%1;%.04;
% kappaV(3)=0.5000;%1.5;%1.5;
% gamma(3)=.95;%.75
% sigma0(3)=.45;%.35


%SDE3 
%thetaV(3)=1.150;%.045;%1;%.04;
%kappaV(3)=1.5000;%1.5;%1.5;
%gamma(3)=.95;%.75
%sigma0(3)=.75;%.35

mu1(3)=kappaV(3)*thetaV(3);
mu2(3)=-kappaV(3);
beta1(3)=0;
beta2(3)=1;


Corr1(1,1)=1;
Corr1(1,2)=.7;
Corr1(1,3)=.8;
Corr1(2,2)=1;
Corr1(2,1)=.7;
Corr1(2,3)=.5;
Corr1(3,1)=.8;
Corr1(3,2)=.5;
Corr1(3,3)=1;

%87%Corr1(1:3,1:3)=0.0;

R=chol(Corr1);

 % R(1:3,1:3)=0;
 % R(1,1)=1;
 % R(2,2)=1;
 % R(3,3)=1;


V0(4)=1.000;%.44;
V00(4)=V0(4);
thetaV(4)=1.0;%.045;%1;%.04;
kappaV(4)=1.000;%1.5;%1.5;
mu1(4)=kappaV(4)*thetaV(4);
mu2(4)=-kappaV(4);
beta1(4)=0;
beta2(4)=1;
gamma(4)=.95;%.75
sigma0(4)=.55;%.45


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
for nn=1:NDims
V(1:paths,nn)=V0(nn);
end
%V(1:paths,1:NDims)=V0(1,1:NDims);

Random1(1:paths/2,1:NDims)=0;

for tt=1:Tt
    
    Random2(1:paths/2,1:NDims)=randn(size(Random1));
    Random2(paths/2+1:paths,1:NDims)=-Random2(1:paths/2,1:NDims);

    %Random2=R'*Random2;
    random1=Random2*R(:,1);
    random2=Random2*R(:,2);
    random3=Random2*R(:,3);
    
    Random2(:,1)=random1;
    Random2(:,2)=random2;
    Random2(:,3)=random3;
    
    
    for nn=1:NDims
    
    V(1:paths,nn)=V(1:paths,nn)+ ...
        (mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn))*dt + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*Random2(1:paths,nn)*sqrt(dt) + ...
        (mu1(nn).*beta1(nn)*V(1:paths,nn).^(beta1(nn)-1) + mu2(nn).*beta2(nn).*V(1:paths,nn).^(beta2(nn)-1)).* ...
        ((mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn))*dt^2/2 + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*Random2(1:paths,nn)*(1-1/sqrt(3))*dt^1.5) + ...
        .5*(mu1(nn).*beta1(nn).*(beta1(nn)-1).*V(1:paths,nn).^(beta1(nn)-2) + mu2(nn).*beta2(nn).*(beta2(nn)-1).*V(1:paths,nn).^(beta2(nn)-2)).* ...    
        sigma0(nn)^2.*V(1:paths,nn).^(2*gamma(nn)).*dt^2/2 + ...
        sigma0(nn).*gamma(nn).*V(1:paths,nn).^(gamma(nn)-1) .* ...
        ((mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn)).*Random2(1:paths,nn).*1/sqrt(3)*dt^1.5 + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*(Random2(1:paths,nn).^2-1)*dt/2) + ...
        .5*sigma0(nn).*gamma(nn).*(gamma(nn)-1).*V(1:paths,nn).^(gamma(nn)-2) .* ...
        sigma0(nn).^2.*V(1:paths,nn).^(2*gamma(nn)) .*Random2(1:paths,nn).*1/sqrt(3)*dt^1.5;

    V(V(:,nn)<=0)=.0001;

    
    if(tt==TT1)
        Xin(1:paths,nn)=real(V(1:paths,nn));
    end
        
    if(tt==TT2)
        
        Xin2(1:paths,nn)=real(V(1:paths,nn));
    end

    end    
    
end


Yin(1:paths)=0.0;

for pp=1:paths
for mm=1:NDims
    Yin(pp)=Yin(pp)+Xin2(pp,mm);%-Xin(pp,mm);
%-Xin(pp,mm);
end
end


str=input("I have reached stone 1");


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% 



%save("C:\Users\Lenovo\Documents\MATLAB\MATLAB1\YXfile.mat","Yin","Xin")


%str=input("Data has been saved");




% load("C:\Users\Lenovo\Documents\MATLAB\MATLAB1\YXfile.mat","Yin","Xin")
% 
% str=input("Data has been loaded");
% 
% 
% 
% 
 Order=7;
 clf;
[BetaCorrH,YCoeffH,XCoeffH,HCorrXX,Zy,Zx] = CalculateHermiteCorrelations0NA04(Xin,Yin,paths,NDims);
YCoeffZ(1:Order+1)= ConvertHermiteSeriesToZSeries01(YCoeffH(1:Order+1),Order);

%HCorrXX(:,:)=0.0;

BetaCorrH


%CorrXX1st
str=input("I have reached Stone2");


%[cn,dn,bnv,bbnv,ch] = Calculate2DZSeriesNDimCorrected(YCoeffH,BetaCorrH,HCorrXX,NDims,Order);
%cn=bnv;
[cn,ch]=CalculateGaussianCorrelationFromHermiteCorrelationsNDim(YCoeffH,XCoeffH,BetaCorrH,HCorrXX,NDims,Order);
% 
% 
% 
% save("C:\Users\Lenovo\Documents\MATLAB\MATLAB1\Zyxfile.mat","Zy","Zx","cn","HCorrXX")
% 
% 
% str=input("Data has been saved");

%load("C:\Users\Lenovo\Documents\MATLAB\MATLAB1\Zyxfile.mat","Zy","Zx","cn","HCorrXX")

%str=input("Data has been loaded");

size(Zy)
size(Zx)
size(cn)


%OptimizeCorrFromHermitesIterative(Zy,Zx,cn,HCorrXX,NDims,paths,Order);

%OptimizeCorrFromHermitesIterative02(Zy,Zx,cn,HCorrXX,NDims,paths,Order);

%OptimizeCorrFromHermitesNewton02(Zy,Zx,cn,HCorrXX,NDims,paths,Order);

[cnOptim] = OptimizeCorrFromHermitesNewton(Zy,Zx,cn,HCorrXX,NDims,paths,Order);

cn=cnOptim;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5555

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Set1
V10(1)=.65;%.65;%
V10(2)=1.2;
V10(3)=.85;
V10(4)=1.85;




% % %Set2
V10(1)=.75;
V10(2)=1.2;%1.2;
V10(3)=.7;
% %V10(4)=.65
% % 
% % %Set3
  V10(1)=1.25;
  V10(2)=1.35;%1.2;
  V10(3)=1.25;
%   V10(4)=1.4
% % 
% %   
  V10(1)=.85;
  V10(2)=.81;%1.2;
  V10(3)=.77;
% %  V10(4)=.7

 

%V10(1)=.65;
%V10(2)=.51;%1.2;
%V10(3)=.73;

%V10(1)=1.65;
%V10(2)=1.51;%1.2;
%V10(3111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111)=1.63;

for nn=1:NDims
V(1:paths,nn)=V10(nn);
end
%V(1:paths,1:NDims)=V0(1,1:NDims);

Random1(1:paths/2,1:NDims)=0;

for tt=1:TT2-TT1
    
    Random2(1:paths/2,1:NDims)=randn(size(Random1));
    Random2(paths/2+1:paths,1:NDims)=-Random2(1:paths/2,1:NDims);

    random1=Random2*R(:,1);
    random2=Random2*R(:,2);
    random3=Random2*R(:,3);
    
    Random2(:,1)=random1;
    Random2(:,2)=random2;
    Random2(:,3)=random3;
    
    
    
    for nn=1:NDims
    
    V(1:paths,nn)=V(1:paths,nn)+ ...
        (mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn))*dt + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*Random2(1:paths,nn)*sqrt(dt) + ...
        (mu1(nn).*beta1(nn)*V(1:paths,nn).^(beta1(nn)-1) + mu2(nn).*beta2(nn).*V(1:paths,nn).^(beta2(nn)-1)).* ...
        ((mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn))*dt^2/2 + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*Random2(1:paths,nn)*(1-1/sqrt(3))*dt^1.5) + ...
        .5*(mu1(nn).*beta1(nn).*(beta1(nn)-1).*V(1:paths,nn).^(beta1(nn)-2) + mu2(nn).*beta2(nn).*(beta2(nn)-1).*V(1:paths,nn).^(beta2(nn)-2)).* ...    
        sigma0(nn)^2.*V(1:paths,nn).^(2*gamma(nn)).*dt^2/2 + ...
        sigma0(nn).*gamma(nn).*V(1:paths,nn).^(gamma(nn)-1) .* ...
        ((mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn)).*Random2(1:paths,nn).*1/sqrt(3)*dt^1.5 + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*(Random2(1:paths,nn).^2-1)*dt/2) + ...
        .5*sigma0(nn).*gamma(nn).*(gamma(nn)-1).*V(1:paths,nn).^(gamma(nn)-2) .* ...
        sigma0(nn).^2.*V(1:paths,nn).^(2*gamma(nn)) .*Random2(1:paths,nn).*1/sqrt(3)*dt^1.5;

    V(V(:,nn)<0)=.0001;

    
    %if(tt==TT1)
    %    XXin(1:paths,nn)=V(1:paths,nn);
    %end
        
    if(tt==TT2-TT1)
        
        XXin2(1:paths,nn)=V(1:paths,nn);
    end

    end    
    
end


Yin2(1:paths)=0.0;

for pp=1:paths
for mm=1:NDims
    Yin2(pp)=Yin2(pp)+XXin2(pp,mm);%-V10(mm);
end
end



randn()

str=input("I have reached stone 4");



NoOfBins=200;
MaxCutOff=20;

[XDensity,IndexOut,IndexMax] = MakeDensityFromSimulation_Infiniti_NEW(Yin2,paths,NoOfBins,MaxCutOff );

for qq=1:8
    
    DataMoments(qq)=sum(Yin2(:).^qq)/paths;
    
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

for nn=1:NDims
    [Zx0(nn)] = CalculateZgivenXAndHSeriesCoeffs(V10(nn),XCoeffH(1:Order+1,nn),Order);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
Zx(1:Order)=0; 

for hh=1:Order
    for nn=1:NDims
        %Zx(hh) = Zx(hh) + bnv(nn+1,hh+1).^1.*Zx0(nn);
        Zx(hh) = Zx(hh) + cn(nn,hh).^1.*Zx0(nn);
    end
end

bh0=YCoeffH(1);

[W] = EvalHermiteArrayH0(Zx0,Order);

for hh=1:Order
    for nn=1:NDims
        bh0 = bh0 + ch(nn+1,hh+1).^1.*W(nn,hh+1);
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


a0=YCoeffZ(1);
a1=YCoeffZ(2);
a2=YCoeffZ(3);
a3=YCoeffZ(4);
a4=YCoeffZ(5);
a5=YCoeffZ(6);
a6=YCoeffZ(7);
a7=YCoeffZ(8);


d1=cn(4,1);
d2=cn(4,2);
d3=cn(4,3);
d4=cn(4,4);
d5=cn(4,5);
d6=cn(4,6);
d7=cn(4,7);


%b0=(a0+a1.* Zx(1) +a2.* Zx(2).^2+a3.* Zx(3).^3+a4* Zx(4).^4+a5 *Zx(5).^5+a6* Zx(6).^6+a7* Zx(7).^7);
b1=(a1* d1+2* a2* Zx(2)* d2 +3* a3* Zx(3).^2* d3+4* a4* Zx(4).^3* d4+5* a5* Zx(5).^4* d5+6* a6* Zx(6).^5* d6+7* a7* Zx(7).^6* d7);
b2=(a2* d2.^2+3* a3* Zx(3)* d3^2 +6* a4* Zx(4).^2* d4^2+10* a5* Zx(5).^3* d5^2+15* a6* Zx(6).^4* d6^2+21* a7* Zx(7).^5* d7^2);
b3=(a3* d3^3+4* a4 *Zx(4)* d4^3 +10* a5* Zx(5).^2* d5^3+20* a6* Zx(6).^3* d6^3+35* a7* Zx(7).^4* d7^3);
b4=(a4* d4^4+5* a5* Zx(5)* d5^4 +15* a6 *Zx(6).^2* d6^4+35* a7* Zx(7).^3 *d7^4);
b5=(a5* d5^5+6* a6* Zx(6)* d6^5 +21* a7 *Zx(7).^2* d7^5);
b6=(a6* d6^6+7* a7* Zx(7)* d7^6 );
b7=(a7* d7^7);


%bh0=b0+b2+3*b4+15*b6;
bh1=b1+3*b3+15*b5+105*b7;
bh2=b2+6*b4+45*b6;
bh3=b3+10*b5+105*b7;
bh4=b4+15*b6;
bh5=b5+21*b7;
bh6=b6;
bh7=b7;


bh(1)=bh0;
bh(2)=bh1;
bh(3)=bh2;
bh(4)=bh3;
bh(5)=bh4;
bh(6)=bh5;
bh(7)=bh6;
bh(8)=bh7;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

YCoeffz(1:Order+1)= ConvertHermiteSeriesToZSeries01(bh(1:Order+1),Order);


[YMomentsModel]=CalculateMomentsOfZSeries(YCoeffz(1),YCoeffz(2:Order+1),Order,Order+1);
 
YMomentsModel
 
DataMoments
 
str=input("Look at comparison of moments");

[cMu1,YCentralMoments]=ConvertRawMomentsToCentralMoments(YMomentsModel,Order+1);
[cMu2,DataCentralMoments]=ConvertRawMomentsToCentralMoments(DataMoments,Order+1);


YCentralMoments

DataCentralMoments

str=input("Look at central moments");

% Order0=7;
% Order1=13;
% NData=paths;
% [CondMomentYgivX,BetaCorrHxy] = CalculateConditionalMomentYgivXNDims(Xin,Yin,NData,Zx0,NDims,Order0,Order1)
% 
% [VarConditional]=FindSecondConditionalMoment(YCoeffH,BetaCorrH,Zx0)
% 
% 
% BetaCorrHxy
% BetaCorrH
% 
% str=input("Look at conditionalVariance");

clf;

PlotZSeriesDensity(YCoeffz(1),YCoeffz(2:Order+1),'b')


hold on

plot(IndexOut(1:IndexMax),XDensity(1:IndexMax),'r');

%PlotHermiteSeriesDensityAndRvGraph(YCoeffh(1),YCoeffh1(2:yOrder+1),'g')
title("Graph of Density of sum of Three Correlated SDE Variables conditional on starting values of the three SDE random variables.")
%title(sprintf('V00 = %.3f; kappa=%.3f, theta=%.3f,gamma=%.3f,sigma=%.3f;V10=%.3f at TT1=%.3f;TT2=%.3f;Paths=%d ',V00, kappaV, thetaV,gamma,sigma0, V10,(TT1*dt),(TT2*dt),paths));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
legend({'Analytical Density','Numerical Density1'}, ...
    'Location','northeast')

str=input("This is the first version without altering the moments of conditional density")

hold off 
 





end

.
.
.
.
.

function [cnOptim] = OptimizeCorrFromHermitesNewton(Zy,Zx,cnn,HCorrXX,NDims,paths,Order)

%Order;Max Order of Hermites Calibration is seven.


ccn(1:NDims+1,1:Order)=cnn(1:NDims+1,2:Order+1);

cn(1:NDims,1:Order)=ccn(2:NDims+1,1:Order);

cn(NDims+1,1:Order)=ccn(1,1:Order);

[WHy] = EvalHermiteArrayH1(Zy,Order);

[EZ] = EZZ2();

for nH=1:Order
    cnIn(1:NDims+1,1:nH)=0.0;
    if(nH==1)
        cnIn(1:NDims+1,1)=cn(1:NDims+1,1);
    else
        cnIn(1:NDims+1,1:nH-1)=cnOptim(1:NDims+1,1:nH-1);
        cnIn(1:NDims+1,nH)=cn(1:NDims+1,nH);
    end
    
    WHyn(1:paths,1)=WHy(1:paths,nH);
    
    Iterations=20;
    
    [cnbest,InitialObj(nH),FinalObj(nH)] = OptimizenthHermiteGaussianWithNewton(cnIn,Zx,WHyn,HCorrXX,nH,EZ,NDims,paths, Iterations);

    cnOptim(1:NDims+1,nH)=cnbest(1:NDims+1,1);
    
    nH
end


InitialObj
FinalObj
str=input("Look at initial and Final Objs");

end

.
.
.
.

function [Wh1] = EvalHermiteArrayH1(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 4 array without the first constant.

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

Wh1(:,1:Order)=Wh(:,2:Order+1);

end

.
.
.

function [bnbest,InitialObj,FinalObj] = OptimizenthHermiteGaussianWithNewton(cn,Zx,WHyn,HCorrXX,nH,EZ,NDims,paths,Iterations)

%Why0 in array of nth Hermite Values

KK=NDims;

bn(1:NDims+1,1)=cn(1:NDims+1,nH);

[dHermite] = EvaluatePrevGaussPowersForHCorrs(nH,cn,Zx,NDims,paths,EZ);

WHy0=WHyn+dHermite;


var=0.0;
for kk1=1:KK
    for kk2=1:KK    
        var=var+bn(kk1,1).*bn(kk2,1).*HCorrXX(kk1,kk2,1);
    end
end

var=var+bn(KK+1,1).^2;

bn

var

bn(:,1)=bn(:,1)/sqrt(var);

%str=input("have reached 1st Stone")        

%[EZ] = EZZ2();

GG1(1:paths,1)=0.0;
for kk=1:KK
    GG1(1:paths,1)=GG1(1:paths,1)+bn(kk,1).*Zx(1:paths,kk);
end

Zeta0=bn(KK+1,1);

[GaussPH,GaussPh] = EvaluateHermitePolyForGaussCorrs(GG1,Zeta0,EZ,nH,paths);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Starting ObjBest and calculations here
%cnbest here
bnbest(1:NDims+1,1)=bn(1:NDims+1);

ObjSqrdBestNewton=sum((GaussPH(:,1)-WHy0(:,1)).^2)/paths;

InitialObj=ObjSqrdBestNewton;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

for iter =1:Iterations

    dGG1(1:paths,1:NDims)=0.0;
    for kk=1:KK
        dGG1(1:paths,kk)=Zx(1:paths,kk);
    end
    
    bn=bnbest;
    
    GG1(1:paths,1)=0.0;
    for kk=1:KK
        GG1(1:paths,1)=GG1(1:paths,1)+bn(kk,1).*Zx(1:paths,kk);
    end

    Zeta0=bn(KK+1,1);
    
    bbn=0;
    bbn(1:NDims,1)=bn(1:NDims,1);

    
    [GaussPH,GaussPh,DGaussPH] = EvaluateHermitePolyAndDerivsForGaussCorrs(GG1,dGG1,Zeta0,EZ,nH,NDims,paths);


    %size(DGaussPH')
    %size(DGaussPH)
    %size(GaussPH)
    %size(WHy0(:,1))
    
    %str=input("Look at sizes");
    
    bbn=bbn+.1*(inv(DGaussPH'*DGaussPH)*(DGaussPH'*(WHy0(:,1)-GaussPH)));
    
    %DGaussPH'*DGaussPH
    
    %str=input("Look at matrix product");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Midfunction ObjBest and calculations here
%cnbest here

var=0.0;
for kk1=1:KK
    for kk2=1:KK    
        var=var+bbn(kk1,1).*bbn(kk2,1).*HCorrXX(kk1,kk2,1);
    end
end

bbn(KK+1,1)=sqrt(1-var);

var=var+bbn(KK+1,1).^2;
bn
bbn

var

bbn=bbn/sqrt(var);



GG1(1:paths,1)=0.0;
for kk=1:KK
    GG1(1:paths,1)=GG1(1:paths,1)+bbn(kk,1).*Zx(1:paths,kk);
end


Zeta0=bbn(KK+1,1);

[GaussPH,GaussPh] = EvaluateHermitePolyForGaussCorrs(GG1,Zeta0,EZ,nH,paths);

ObjSqrdNewton=sum((GaussPH(:,1)-WHy0(:,1)).^2)/paths;
ObjSqrdBestNewton
if(ObjSqrdNewton<ObjSqrdBestNewton)
    ObjSqrdBestNewton=ObjSqrdNewton;
    bnbest=bbn;
else
    %%%%change-step or Newton parameters
end
    
ObjSqrdBestNewton
ObjSqrdNewton
iter

%str=input("Look at the solution");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5    

end

FinalObj=ObjSqrdBestNewton;

end

.
.
.

function [dHermite] = EvaluatePrevGaussPowersForHCorrs(nH,ccn,Zx,NDims,paths,EZ)

if(nH==1)
    dHermite(1:paths,1)=0.0;
end

if(nH==2)
    dHermite(1:paths,1)=1.0;
end

if(nH==3)
    nH1=1;
    [GaussPH1] = EvaluateGaussianPowerForHCorrs(ccn,nH1,Zx,NDims,paths,EZ);
    dHermite(1:paths,1)=3*GaussPH1(1:paths,1);
end

if(nH==4)
     nH1=2;
    [GaussPH1] = EvaluateGaussianPowerForHCorrs(ccn,nH1,Zx,NDims,paths,EZ);
     dHermite(1:paths,1)=6*GaussPH1(1:paths,1)-3;
end

if(nH==5)
    nH1=1;
    [GaussPH1] = EvaluateGaussianPowerForHCorrs(ccn,nH1,Zx,NDims,paths,EZ);
    
    nH2=3;
    [GaussPH2] = EvaluateGaussianPowerForHCorrs(ccn,nH2,Zx,NDims,paths,EZ);
    dHermite(1:paths,1)=-(-10*GaussPH2(1:paths,1)+15*GaussPH1(1:paths,1));
end

if(nH==6)
    nH1=2;
    [GaussPH1] = EvaluateGaussianPowerForHCorrs(ccn,nH1,Zx,NDims,paths,EZ);
    
    nH2=4;
    [GaussPH2] = EvaluateGaussianPowerForHCorrs(ccn,nH2,Zx,NDims,paths,EZ);
    dHermite(1:paths,1)=-(-15*GaussPH2(1:paths,1)+45*GaussPH1(1:paths,1)-15);
end

if(nH==7)
    nH1=1;
    [GaussPH1] = EvaluateGaussianPowerForHCorrs(ccn,nH1,Zx,NDims,paths,EZ);
    
    nH2=3;
    [GaussPH2] = EvaluateGaussianPowerForHCorrs(ccn,nH2,Zx,NDims,paths,EZ);
    
    nH3=5;
    [GaussPH3] = EvaluateGaussianPowerForHCorrs(ccn,nH3,Zx,NDims,paths,EZ);
    
    dHermite(1:paths,1)=-(-21*GaussPH3(1:paths,1)+105*GaussPH2(1:paths,1)-105*GaussPH1(1:paths,1));
end


end

.
.
.

function [GaussPH] = EvaluateGaussianPowerForHCorrs(ccn,nH,Zx,NDims,paths,EZ)

KK=NDims;

GG1(1:paths,1)=0.0;

for kk=1:KK
    GG1(1:paths,1)=GG1(1:paths,1)+ccn(kk,nH).*Zx(1:paths,kk);
end
Zeta0=ccn(KK+1,nH);
[GaussPH,GaussPh] = EvaluateHermitePolyForGaussCorrs(GG1,Zeta0,EZ,nH,paths);

end

.
.
.

function [GaussPH,GaussPh] = EvaluateHermitePolyForGaussCorrs(GG1,Zeta0,EZ,HH,paths)
GaussPH(1:paths,1)=0.0;
GaussPh(1:paths,1:HH+1)=0.0;

for hh=0:HH
    GaussPh(1:paths,hh+1)=factorial(HH)/factorial(hh)/factorial(HH-hh).*GG1(1:paths).^(HH-hh).*Zeta0.^(hh).*EZ(hh+1);
    GaussPH(1:paths,1)=GaussPH(1:paths,1)+GaussPh(1:paths,hh+1);
end
end

.
.
.

function [GaussPH,GaussPh,DGaussPH] = EvaluateHermitePolyAndDerivsForGaussCorrs(GG1,dGG1,Zeta0,EZ,HH,NDims,paths)


GaussPH(1:paths,1)=0.0;
GaussPh(1:paths,1:HH+1)=0.0;
DGaussPh(1:paths,1:HH+1)=0.0;
DGaussPh0(1:paths,1:HH+1)=0.0;
%dGG1(1:paths,1:NDims)=0;
DGaussPH(1:paths,1:NDims)=0.0;
for hh=0:HH
    GaussPh(1:paths,hh+1)=factorial(HH)/factorial(hh)/factorial(HH-hh).*GG1(1:paths).^(HH-hh).*Zeta0.^(hh).*EZ(hh+1);
    GaussPH(1:paths,1)=GaussPH(1:paths,1)+GaussPh(1:paths,hh+1);
    DGaussPh(1:paths,hh+1)=factorial(HH)/factorial(hh)/factorial(HH-hh).*(HH-hh).*dGG1(1:paths).^(HH-hh-1).*Zeta0.^(hh).*EZ(hh+1);
    DGaussPh0(1:paths,hh+1)=factorial(HH)/factorial(hh)/factorial(HH-hh).*GG1(1:paths).^(HH-hh).*hh.*Zeta0.^(hh-1).*EZ(hh+1);
    DGaussPH(1:paths,1:NDims)=DGaussPH(1:paths,1:NDims)+DGaussPh(1:paths,hh+1).*dGG1(1:paths,1:NDims);
    %DGaussPH(1:paths,NDims+1)=DGaussPH(1:paths,NDims+1)+DGaussPh0(1:paths,hh+1);
end

%DGaussPh0
%DGaussPh0(1:paths,2)
%Zeta0

%str=input("Look at DGauss Ph0");


end

.
.
.

function [EZ] = EZZ2()

EZ(1)=1;
EZ(2)=0;
EZ(3)=1;
for nn=4:128 %Increase for larger products.
    if rem(nn,2)==0
        EZ(nn)=0;
    else
        EZ(nn)=EZ(nn-2)*(nn-2);
        EZ(nn);
    end
end
% EZ(64)
% EZ(64).*EZ(64)

end

.
.
.

Amin · May 24th, 2026, 9:07 am

Friends, I will come back with a detailed latex post and commented code later tonight.

But very quickly for friends trying to play with the program, I have changed the order of variables and indexing and it is different in new output variable cnOptim that comes out of Newton optimization. Please let me explain here. In the expansion of nth power of Gaussian, if there are K correlated variables, the correlation coefficients are denoted as

\[{(\,\rho_{1,m}\,Z_1\,+\,\rho_{2,m}\,Z_2\,+\,\ldots\,+\,\rho_{K,m}\,Z_K\,+\,\zeta_m)}^m\]

earlier in previous function that was written weeks ago, cn they were arranged as

\[cn(1,m+1)\,=\,\zeta_m\]
and
\[cn(k+1,m+1)\,=\,\rho_{k,m}\]

however, in the new out put variable cnOptim, zeta that was at first row, has gone to last row. the indexing for cnOptim is given as (Please note that K is the last index related to explainingvaribels)

\[cnOptim(k,m)\,=\,\rho_{k,m}\]
and
\[cnOptim(K+1,m)\,=\,\zeta_m\]

As you can see , In my later program, I have changed indices of variables to reflect the change in indexing of cnOptim. If you want to plug my code in your program, you would have to slightly change this indexing.

Especially again zeta variable that was in first row earlier has gone to the last (K+1)th row (resulting in one step upward shift in rows associated with each correlation coefficient). And indexing of correlation coefficients and zeta of mth power that was earlier at column m+1 has now gone to column m.

Please check how I have changed indexing and position of zeta_m in my new program as I altered them to comply with the new changes and therefore the program works well. If you do not want to do anything, you might want to copy the last part of my program into your code where I have changed everything so that new changes in indexing continue to work well with slightly modified density construction program.

My program also mentions starting value of Objective function when optimization started and ending value of objective function when Newton optimization ended and the difference is mostly huge in case of larger hermites. But for difficult problem, you might want different starting values for Newton since I suppose it is a Gradient descent type algorithm.

You might try extremely small steps of Newton with very large number of iterations for very difficult problems but starting seed has to be good.

I am sure in coming weeks friends will do all sort of new things and experiments with the algorithm.

Amin · May 24th, 2026, 4:41 pm

Friends, I am commenting the new program and will post a commented version in a few hours. In the meantime, I want to mention how to calculate the model values of hermite polynomials that are matched with data values of hermite polynomials. And also describe how to calculate its derivatives for Newton.

But first please look at the details of Newton method in a previous post 3098 here: https://forum.wilmott.com/viewtopic.php?t=99702&start=3090#p890259

In that post I had described that I want to match second power/moment of hermite polynomials simultaneously with its first power but I have not dealt with the second moment in this program I posted yesterday. It works only by matching the values of hermite polynomials by optimizing over the correlation coefficients. It still needs to be seen if matching squared values adds to robustness of the regression and I would continue to see that in coming week. The current program just matches the values of first powers of hermites directly. I will continue to play (over next week) with the idea of including second power later and try to see if it improves the optimization in any way.

Here I want to describe how to calculate model values of hermite polynomials for each observation.

First we optimize the only highest power of Gaussian in each hermite and lower order Gaussians (in higher hermites) are calculated from the correlation coefficients associated with lower order hermites. For example consider fifth hermite given below in terms of its Gaussians expansion as

\[H_5(Z_{y_p})\,=\,{\left( \rho_{1,5}\, Z_{p,1}\,+\, \rho_{2,5}\, Z_{p,2}\,+\,\ldots\,+\,\rho_{K,5}\, Z_{p,K}\,+\,\zeta_5\,Z_{\tilde{y}}\right)}^5\,\\-10\,{\left( \rho_{1,3}\, Z_{p,1}\,+\, \rho_{2,3}\, Z_{p,2}\,+\,\ldots\,+\,\rho_{K,3}\, Z_{p,K}\,+\,\zeta_3\,Z_{\tilde{y}}\right)}^3\\+15\,{\left( \rho_{1,1}\, Z_{p,1}\,+\, \rho_{2,1}\, Z_{p,2}\,+\,\ldots\,+\,\rho_{K,1}\, Z_{p,K}\,+\,\zeta_1\,Z_{\tilde{y}}\right)}\]

Basically when reach optimization of fifth hermite, we have already calculated the correlation coefficients in third power multinomial and first power multinomial from third and first hermite respectively already. We just subtract the first and third Gaussian power multinomials from the data values of fifth hermite and are left with finding the correlation coefficients in the fifth Gaussian power multinomial only.

Now I will describe how to calculate mth power of the Gaussian multinomial for a general mth hermite.

\[\,{\left( \rho_{1,m}\, Z_{p,1}\,+\, \rho_{2,m}\, Z_{p,2}\,+\,\ldots\,+\,\rho_{K,m}\, Z_{p,K}\,+\,\zeta_m\,Z_{\tilde{y}}\right)}^m\,\\=\, {\left(G\,+\,\,\zeta_m\,Z_{\tilde{y}}\right)}^m\]
where
\[\,G\,=\,\left( \rho_{1,m}\, Z_{p,1}\,+\, \rho_{2,m}\, Z_{p,2}\,+\,\ldots\,+\,\rho_{K,m}\, Z_{p,K}\,\right)\]

Please know that all variables and parameters required to calculate G are already known either from data or their trial values in Newton iterations. So we can basically just calculate one G and use binomial expansion as

\[\,{\left( \rho_{1,m}\, Z_{p,1}\,+\, \rho_{2,m}\, Z_{p,2}\,+\,\ldots\,+\,\rho_{K,m}\, Z_{p,K}\,+\,\zeta_m\,Z_{\tilde{y}}\right)}^m\,\\=\, {\left(G\,+\,\,\zeta_m\,Z_{\tilde{y}}\right)}^m \\
=\,\sum\limits_{n=0}^{n=m}\,\frac{m!}{(m-n)!\,n!}G^{m-n}\,{\left(\zeta_m\,\right)}^{n}\,E\left[{\,Z_{\tilde{y}}}^{n}\right]\]

I take expectation over powers of \(Z_{\tilde{y}}\) and the above calculation of expectation is done over all data paths to calculate model values of the mth power of Gaussian that is regressed on mth Hermite. This value is different for every path since G is different for all paths in the data.

Similarly for Newton derivatives of above expectation with respect to a general correlation coefficient \(\rho_{k,m}\), we will use the simple formula

\[\,\sum\limits_{n=0}^{n=m}\,\frac{m!}{(m-n)!\,n!}\,\left(m-n\right)\,\,G^{m-n-1}\,\frac{dG}{d\rho_{k,m}}\,{\left(\zeta_m\,\right)}^{n}\,E\left[{\,Z_{\tilde{y}}}^{n}\right]\]

This last formula is how I calculated sensitivities in Newton method. \(\zeta_m\) was a constrained parameter in the model and was not directly optimized and its value was found from constraint that total variance of the Gaussian term goes to one i.e.

\[\,{\left( \rho_{1,m}\, Z_{p,1}\,+\, \rho_{2,m}\, Z_{p,2}\,+\,\ldots\,+\,\rho_{K,m}\, Z_{p,K}\,+\,\zeta_m\,Z_{\tilde{y}}\right)}^2\,=\,1\,\]

Please also note that \(Z_{\tilde{y}}\) is uncorrelated with all of the K explaining variables.
In the above expansion of second power, cross correlations between different explaining \(Z_x\)'s also appears. I just used first order correlation from data for all the correlations of the form

\[\rho_{k_1,k_2}\,=\,\sum\limits_{p=1}^{paths}\,Z_{p,k_1}\,\,Z_{p,k_2}\]

I had the temptation that above cross-correlation should be different for each mth Gaussian power expansion but then I was thinking that when substitute the data of \(Z_{k_1}\) and \(Z_{k_2}\), they must have the correlation described in above formula.

Amin · May 24th, 2026, 6:42 pm

Here is the matlab code with light comments to explain the logic of Newton method. These are only functions related to over-constrained Newton root search.
.
.
.

function [cnOptim] = OptimizeCorrFromHermitesNewton(Zy,Zx,cnn,HCorrXX,NDims,paths,Order)

%Order;Max Order of Hermites Calibration is seven.
%Zx is paths X K array for joint values of Zx arranged for "paths" number of observations
%Zy is paths X 1 array for Zy associated with each obsewrvation.
%cnn is the original input of correlation coefficients used to derive seed for Newton methjod.
%cnn is in format of the old program and its order and indexing is changed
%in this NEwton program.

ccn(1:NDims+1,1:Order)=cnn(1:NDims+1,2:Order+1); %Shift the indices back by one index.
%so that index matches with order of associated hermite polynomial.

cn(1:NDims,1:Order)=ccn(2:NDims+1,1:Order);%The first row in each hermite column earlier was 
%associated zeta but zeta is shifted from first row to (K+1)th row.

cn(NDims+1,1:Order)=ccn(1,1:Order);


%The function below calculates the first hermites till "Order" from data of
%Zy. WHy is (paths X K) array. kth hermite associated with pth observation appears in
%WHy(p,k)
[WHy] = EvalHermiteArrayH1(Zy,Order);

%The program below calculates moments of standard normal. Odd moments are
%zero. 
[EZ] = EZZ2();

%The loop below repeats itself for optimization of correlation coefficients 
%associated with highest power in each hermite. It repeatrs over all hermites one by one. 
for nH=1:Order
    cnIn(1:NDims+1,1:nH)=0.0;
    if(nH==1)
        cnIn(1:NDims+1,1)=cn(1:NDims+1,1);
    else
        %cnOptim is already optimized value of lower powers of Gaussians
        %calculated in earlier loops of optimization over smaller hermites.
        %This is already known from previous iterations of this loop.
        cnIn(1:NDims+1,1:nH-1)=cnOptim(1:NDims+1,1:nH-1);
        %The value of cnIn below is initial guess for correlation
        %coefficients associated with optimization of current hermite
        cnIn(1:NDims+1,nH)=cn(1:NDims+1,nH);
    end
    
    %Below nH-th hermite observation values are assigned to one dimensional array for optimization target in Newton. 
    WHyn(1:paths,1)=WHy(1:paths,nH);
    
    %Total iterations of Newton root search.
    Iterations=50;
    %Below is the function that actually does the Multi-dimensional Newton.
    [cnbest,InitialObj(nH),FinalObj(nH)] = OptimizenthHermiteGaussianWithNewton(cnIn,Zx,WHyn,HCorrXX,nH,EZ,NDims,paths, Iterations);

    %Below is optimum value of correlation coefficients associated with
    %nH-th hermite after optimization. This will be used by the program as
    %an input to optimization of higher hermites
    cnOptim(1:NDims+1,nH)=cnbest(1:NDims+1,1);
    
    nH
end


InitialObj %Obj function at start of optimization of correlation coefficients for nH-th hermite
FinalObj%Obj function after optimization of correlation coefficients for nH-th hermite
str=input("Look at initial and Final Objs");

end

.
.
.

function [bnbest,InitialObj,FinalObj] = OptimizenthHermiteGaussianWithNewton(cn,Zx,WHyn,HCorrXX,nH,EZ,NDims,paths,Iterations)

%Why0 in array of nth Hermite Values

KK=NDims;

bn(1:NDims+1,1)=cn(1:NDims+1,nH);
%Below, Calculate the model contribution of lower Gaussian powers that have
%already been calculated through associated lower hermites towards current
%hermite polynomial.
[dHermite] = EvaluatePrevGaussPowersForHCorrs(nH,cn,Zx,NDims,paths,EZ);
%Remove the model contribution of lower Gaussian powers from data value of
%current hermite.
WHy0=WHyn+dHermite;

%Below calcualte variance
var=0.0;
for kk1=1:KK
    for kk2=1:KK    
        var=var+bn(kk1,1).*bn(kk2,1).*HCorrXX(kk1,kk2,1);
    end
end

var=var+bn(KK+1,1).^2;

bn

var
%Below, make sure that variance is one.
bn(:,1)=bn(:,1)/sqrt(var);

%str=input("have reached 1st Stone")        


%Below calculate the combined value of K correlated variables as described
%on Wilmott post and Calculate G.
GG1(1:paths,1)=0.0;
for kk=1:KK
    GG1(1:paths,1)=GG1(1:paths,1)+bn(kk,1).*Zx(1:paths,kk);
end

Zeta0=bn(KK+1,1);
%Evaluate the nHth Gaussian power below.
[GaussPH,GaussPh] = EvaluateHermitePolyForGaussCorrs(GG1,Zeta0,EZ,nH,paths);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Starting ObjBest and calculations here
%cnbest here
bnbest(1:NDims+1,1)=bn(1:NDims+1);

ObjSqrdBestNewton=sum((GaussPH(:,1)-WHy0(:,1)).^2)/paths;

InitialObj=ObjSqrdBestNewton;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%Below, enter the Newton loop.
for iter =1:Iterations

    dGG1(1:paths,1:NDims)=0.0;
    for kk=1:KK
        dGG1(1:paths,kk)=Zx(1:paths,kk);
    end
    
    bn=bnbest;
    
    GG1(1:paths,1)=0.0;
    for kk=1:KK
        GG1(1:paths,1)=GG1(1:paths,1)+bn(kk,1).*Zx(1:paths,kk);
    end

    Zeta0=bn(KK+1,1);
    
    bbn=0;
    bbn(1:NDims,1)=bn(1:NDims,1);

    
    [GaussPH,GaussPh,DGaussPH] = EvaluateHermitePolyAndDerivsForGaussCorrs(GG1,dGG1,Zeta0,EZ,nH,NDims,paths);


    %size(DGaussPH')
    %size(DGaussPH)
    %size(GaussPH)
    %size(WHy0(:,1))
    
    %str=input("Look at sizes");
    
    %I am using a multiple of .1 in Newton equation below. this might have
    %to be further decreased for some cases. Below is the Newton
    %over-constrained iteration.
    bbn=bbn+.1*((DGaussPH'*DGaussPH)\(DGaussPH'*(WHy0(:,1)-GaussPH)));
    
    %DGaussPH'*DGaussPH
    
    %str=input("Look at matrix product");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Midfunction ObjBest and calculations here
%cnbest here

var=0.0;
for kk1=1:KK
    for kk2=1:KK    
        var=var+bbn(kk1,1).*bbn(kk2,1).*HCorrXX(kk1,kk2,1);
    end
end

bbn(KK+1,1)=sqrt(1-var);

var=var+bbn(KK+1,1).^2;
bn
bbn

var

bbn=bbn/sqrt(var);



GG1(1:paths,1)=0.0;
for kk=1:KK
    GG1(1:paths,1)=GG1(1:paths,1)+bbn(kk,1).*Zx(1:paths,kk);
end


Zeta0=bbn(KK+1,1);

[GaussPH,GaussPh] = EvaluateHermitePolyForGaussCorrs(GG1,Zeta0,EZ,nH,paths);

ObjSqrdNewton=sum((GaussPH(:,1)-WHy0(:,1)).^2)/paths;
ObjSqrdBestNewton
if(ObjSqrdNewton<ObjSqrdBestNewton)
    ObjSqrdBestNewton=ObjSqrdNewton;
    bnbest=bbn;
else
    %%%%change-step or Newton parameters
end
    
ObjSqrdBestNewton
ObjSqrdNewton
iter

%str=input("Look at the solution");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5    

end

FinalObj=ObjSqrdBestNewton;

end

.
.
.

function [GaussPH] = EvaluateGaussianPowerForHCorrs(ccn,nH,Zx,NDims,paths,EZ)

KK=NDims;

GG1(1:paths,1)=0.0;
%Below calculate the binomial expansion related to nH-th power of the sum
%of correlated Gaussians
for kk=1:KK
    GG1(1:paths,1)=GG1(1:paths,1)+ccn(kk,nH).*Zx(1:paths,kk);
end
Zeta0=ccn(KK+1,nH);
[GaussPH,GaussPh] = EvaluateHermitePolyForGaussCorrs(GG1,Zeta0,EZ,nH,paths);

end

.
.
.

function [dHermite] = EvaluatePrevGaussPowersForHCorrs(nH,ccn,Zx,NDims,paths,EZ)
%This function clculates the contribution of lower Gaussian powers of sum
%of K correlated Zx's towards the value of higher hermites. Please note
%that sign is reversed so they are basically added in the calling function
%to subtract them.
if(nH==1)
    dHermite(1:paths,1)=0.0;
end

if(nH==2)
    dHermite(1:paths,1)=1.0;
end

if(nH==3)
    nH1=1;
    [GaussPH1] = EvaluateGaussianPowerForHCorrs(ccn,nH1,Zx,NDims,paths,EZ);
    dHermite(1:paths,1)=3*GaussPH1(1:paths,1);
end

if(nH==4)
     nH1=2;
    [GaussPH1] = EvaluateGaussianPowerForHCorrs(ccn,nH1,Zx,NDims,paths,EZ);
     dHermite(1:paths,1)=6*GaussPH1(1:paths,1)-3;
end

if(nH==5)
    nH1=1;
    [GaussPH1] = EvaluateGaussianPowerForHCorrs(ccn,nH1,Zx,NDims,paths,EZ);
    
    nH2=3;
    [GaussPH2] = EvaluateGaussianPowerForHCorrs(ccn,nH2,Zx,NDims,paths,EZ);
    dHermite(1:paths,1)=-(-10*GaussPH2(1:paths,1)+15*GaussPH1(1:paths,1));
end

if(nH==6)
    nH1=2;
    [GaussPH1] = EvaluateGaussianPowerForHCorrs(ccn,nH1,Zx,NDims,paths,EZ);
    
    nH2=4;
    [GaussPH2] = EvaluateGaussianPowerForHCorrs(ccn,nH2,Zx,NDims,paths,EZ);
    dHermite(1:paths,1)=-(-15*GaussPH2(1:paths,1)+45*GaussPH1(1:paths,1)-15);
end

if(nH==7)
    nH1=1;
    [GaussPH1] = EvaluateGaussianPowerForHCorrs(ccn,nH1,Zx,NDims,paths,EZ);
    
    nH2=3;
    [GaussPH2] = EvaluateGaussianPowerForHCorrs(ccn,nH2,Zx,NDims,paths,EZ);
    
    nH3=5;
    [GaussPH3] = EvaluateGaussianPowerForHCorrs(ccn,nH3,Zx,NDims,paths,EZ);
    
    dHermite(1:paths,1)=-(-21*GaussPH3(1:paths,1)+105*GaussPH2(1:paths,1)-105*GaussPH1(1:paths,1));
end


end

.
.
.

function [GaussPH,GaussPh,DGaussPH] = EvaluateHermitePolyAndDerivsForGaussCorrs(GG1,dGG1,Zeta0,EZ,HH,NDims,paths)
%this function calculates the binomial expansion terms and their combined
%expectations. The binaomial expansion is over G and zeta Z-tilde{y} as
%described in Wilmott post. This function also calculates the sensitivities
%of the Gaussian power terms expanded as binomial with respect to
%correlation parameters.

GaussPH(1:paths,1)=0.0;
GaussPh(1:paths,1:HH+1)=0.0;
DGaussPh(1:paths,1:HH+1)=0.0;
DGaussPh0(1:paths,1:HH+1)=0.0;
%dGG1(1:paths,1:NDims)=0;
DGaussPH(1:paths,1:NDims)=0.0;
for hh=0:HH
    GaussPh(1:paths,hh+1)=factorial(HH)/factorial(hh)/factorial(HH-hh).*GG1(1:paths).^(HH-hh).*Zeta0.^(hh).*EZ(hh+1);
    GaussPH(1:paths,1)=GaussPH(1:paths,1)+GaussPh(1:paths,hh+1);
    DGaussPh(1:paths,hh+1)=factorial(HH)/factorial(hh)/factorial(HH-hh).*(HH-hh).*dGG1(1:paths).^(HH-hh-1).*Zeta0.^(hh).*EZ(hh+1);
    DGaussPh0(1:paths,hh+1)=factorial(HH)/factorial(hh)/factorial(HH-hh).*GG1(1:paths).^(HH-hh).*hh.*Zeta0.^(hh-1).*EZ(hh+1);
    DGaussPH(1:paths,1:NDims)=DGaussPH(1:paths,1:NDims)+DGaussPh(1:paths,hh+1).*dGG1(1:paths,1:NDims);
    %DGaussPH(1:paths,NDims+1)=DGaussPH(1:paths,NDims+1)+DGaussPh0(1:paths,hh+1);
end

%DGaussPh0
%DGaussPh0(1:paths,2)
%Zeta0

%str=input("Look at DGauss Ph0");


end

Amin · Yesterday, 9:50 am

Friends, I noticed that cases of negative correlations and zero correlations between explaining variables are not being properly captured by the model. I think this is because, I naively set the cross correlations to first order cross-correlations as mentioned at the end of previous post. It was working so well for positive correlation cases that I thought this was probably right.

I think I know how to fix the problem and, in another day or two, I hope to come with a new version of matlab program that fixes this problem.

Amin · Yesterday, 1:16 pm

Friends, It is always like that whenever I do some good work or good research, they increase torture on me manifold. Now I am trying to work to fix the program and mind control agents have focused high frequency lasers on my temple and it causes very sharp pain in my temple. They have previously done that several times in the past and there were many times when they continued to focus lasers on my temple for 5-6 days in a row and I suffered from agony of very sharp pain. I tried to plead to mind control agents to not focus lasers on my temple but they respond by increasing the torture.

More than ten years ago when I started doing my research, I used to think that doing good research would decrease torture on me but now I am always prepared since I know whenever I do good research, they always, as a rule. respond by increasing the torture.

Please protest to scumbags of defense and ask them to learn to respect dignity of good people.

Amin · Yesterday, 7:21 pm

Friends, I have made a few changes to the program algorithm. And now zero correlations and mixed correlations are just as precise as positive correlations. But for the case in which all explaining variables have negative correlations, the results are still off.

I found that conditional mean from separate regressions for each hermite is quite less precise than conditional mean from the single regression in which all X variables based hermites are in a single least square regression against Y. Therefore I want to find Newton guess coefficients also from coefficients found from one single least square regression. I also want to normalize this regression with square-root of variance so we would have a very good estimate of correlations as a start. In the earlier algorithm, we were taking seventh root and sixth root of coefficients of sixth and seventh hermite polynomial but this was basically sixth and seventh root of covariances and not correlations. We needed to normalize coefficients of hermites of explaining X variables in our equation of Y with variance of Y. We have been taking roots of covariances earlier rather than roots of correlations. I will change this in the new program.

There is huge improvement on the program I posted two days ago and I hope that I would be able to fix the all negative correlations case as well sometime later tonight.

Amin · Yesterday, 11:13 pm

Friends, a lot of work needs to be done but I am posting some preliminary programs. Conditional Mean is especially improved where it was bad.
,
,
,

function []= ConditionalDensityHermitesNvariateNewton()

%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
%The present program works only when drift is zero (for pure noises)
%meaning mu1=0 and mu2=0;

Order=7;
%NDim=4;%Three assets and one SV.
%NMomentsY=Order+1;
%NMomentsX=Order+1;

%%%%%%%%%%%%%%%%%%%5
dt=.03125/2/8;
Tt=64*2*8;
T=Tt*dt;
seed0=12130649;
%seed0=94210876;
%seed0=14210876;
%seed0=53612297
seed0=97459782
rng(seed0, 'twister')
paths=50000;
TT1=100*8/2; %Transition distribution start
TT2=256*8/2; %Transition distribution end

%%%%%%%%%%%%%%%%%%%5

NDims=3;
% 
V0(1)=1.000;%.32;
V00(1)=V0(1);

%SDE1
 thetaV(1)=1.0;%.045;%1;%.04;
 kappaV(1)=0.000;%1.5;%1.5;
 gamma(1)=.75;%.75
 sigma0(1)=.5;%.35


 
%SDE2%
%thetaV(1)=.90;%.045;%1;%.04;
%kappaV(1)=1.000;%1.5;%1.5;
%gamma(1)=.75;%.75
%sigma0(1)=.85;%.35

%SDE3%
%thetaV(1)=1.350;%.045;%1;%.04;
%kappaV(1)=2.000;%1.5;%1.5;
%gamma(1)=.75;%.75
%sigma0(1)=.95;%.35


mu1(1)=kappaV(1)*thetaV(1);
mu2(1)=-kappaV(1);
beta1(1)=0;
beta2(1)=1;



V0(2)=1.000;%.32;
V00(2)=V0(2);

%SDE1
thetaV(2)=1.0;%.045;%1;%.04;
 kappaV(2)=0.000;%1.5;%1.5;
 gamma(2)=.65;%.75
 sigma0(2)=.5;%.35

%SDE2 
% thetaV(2)=1.0;%.045;%1;%.04;
% kappaV(2)=2.000;%1.5;%1.5;
% gamma(2)=.65;%.75
% sigma0(2)=.95;%.35

%SDE3
%thetaV(2)=.850;%.045;%1;%.04;
%kappaV(2)=.35000;%1.5;%1.5;
%gamma(2)=.65;%.75
%sigma0(2)=.55;%.35


mu1(2)=kappaV(2)*thetaV(2);
mu2(2)=-kappaV(2);
beta1(2)=0;
beta2(2)=1;



V0(3)=1.000;%.33;
V00(3)=V0(3);

%SDE1
thetaV(3)=1.0;%.045;%1;%.04;
 kappaV(3)=0.000;%1.5;%1.5;
 gamma(3)=.95;%.75
 sigma0(3)=.45;%.35

%SDE2 
% thetaV(3)=.750;%.045;%1;%.04;
% kappaV(3)=0.5000;%1.5;%1.5;
% gamma(3)=.95;%.75
% sigma0(3)=.45;%.35


%SDE3 
%thetaV(3)=1.150;%.045;%1;%.04;
%kappaV(3)=1.5000;%1.5;%1.5;
%gamma(3)=.95;%.75
%sigma0(3)=.75;%.35

mu1(3)=kappaV(3)*thetaV(3);
mu2(3)=-kappaV(3);
beta1(3)=0;
beta2(3)=1;



Corr1(1,1)=1;
Corr1(1,2)=.7;
Corr1(1,3)=.8;
Corr1(2,2)=1;
Corr1(2,1)=.7;
Corr1(2,3)=.5;
Corr1(3,1)=.8;
Corr1(3,2)=.5;
Corr1(3,3)=1;


% Corr1(1,1)=1;
% Corr1(1,2)=-.7/3;
% Corr1(1,3)=.8/3;
% Corr1(2,2)=1;
% Corr1(2,1)=-.7/3;
% Corr1(2,3)=-.5/3;
% Corr1(3,1)=.8/3;
% Corr1(3,2)=-.5/3;
% Corr1(3,3)=1;



% Corr1(1,1)=1;
% Corr1(1,2)=-.27;
% Corr1(1,3)=-.28;
% Corr1(2,2)=1;
% Corr1(2,1)=-.27;
% Corr1(2,3)=-.15;
% Corr1(3,1)=-.28;
% Corr1(3,2)=-.15;
% Corr1(3,3)=1;
%87%Corr1(1:3,1:3)=0.0;

R=chol(Corr1);
% 
%   R(1:3,1:3)=0;
%   R(1,1)=1;
%   R(2,2)=1;
%   R(3,3)=1;


V0(4)=1.000;%.44;
V00(4)=V0(4);
thetaV(4)=1.0;%.045;%1;%.04;
kappaV(4)=1.000;%1.5;%1.5;
mu1(4)=kappaV(4)*thetaV(4);
mu2(4)=-kappaV(4);
beta1(4)=0;
beta2(4)=1;
gamma(4)=.95;%.75
sigma0(4)=.55;%.45


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
for nn=1:NDims
V(1:paths,nn)=V0(nn);
end
%V(1:paths,1:NDims)=V0(1,1:NDims);

Random1(1:paths/2,1:NDims)=0;

for tt=1:Tt
    
    Random2(1:paths/2,1:NDims)=randn(size(Random1));
    Random2(paths/2+1:paths,1:NDims)=-Random2(1:paths/2,1:NDims);

    %Random2=R'*Random2;
    random1=Random2*R(:,1);
    random2=Random2*R(:,2);
    random3=Random2*R(:,3);
    
    Random2(:,1)=random1;
    Random2(:,2)=random2;
    Random2(:,3)=random3;
    
    
    for nn=1:NDims
    
    V(1:paths,nn)=V(1:paths,nn)+ ...
        (mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn))*dt + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*Random2(1:paths,nn)*sqrt(dt) + ...
        (mu1(nn).*beta1(nn)*V(1:paths,nn).^(beta1(nn)-1) + mu2(nn).*beta2(nn).*V(1:paths,nn).^(beta2(nn)-1)).* ...
        ((mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn))*dt^2/2 + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*Random2(1:paths,nn)*(1-1/sqrt(3))*dt^1.5) + ...
        .5*(mu1(nn).*beta1(nn).*(beta1(nn)-1).*V(1:paths,nn).^(beta1(nn)-2) + mu2(nn).*beta2(nn).*(beta2(nn)-1).*V(1:paths,nn).^(beta2(nn)-2)).* ...    
        sigma0(nn)^2.*V(1:paths,nn).^(2*gamma(nn)).*dt^2/2 + ...
        sigma0(nn).*gamma(nn).*V(1:paths,nn).^(gamma(nn)-1) .* ...
        ((mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn)).*Random2(1:paths,nn).*1/sqrt(3)*dt^1.5 + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*(Random2(1:paths,nn).^2-1)*dt/2) + ...
        .5*sigma0(nn).*gamma(nn).*(gamma(nn)-1).*V(1:paths,nn).^(gamma(nn)-2) .* ...
        sigma0(nn).^2.*V(1:paths,nn).^(2*gamma(nn)) .*Random2(1:paths,nn).*1/sqrt(3)*dt^1.5;

    V(V(:,nn)<=0)=.0001;

    
    if(tt==TT1)
        Xin(1:paths,nn)=real(V(1:paths,nn));
    end
        
    if(tt==TT2)
        
        Xin2(1:paths,nn)=real(V(1:paths,nn));
    end

    end    
    
end


Yin(1:paths)=0.0;

for pp=1:paths
for mm=1:NDims
    Yin(pp)=Yin(pp)+Xin2(pp,mm);%-Xin(pp,mm);
%-Xin(pp,mm);
end
end


str=input("I have reached stone 1");


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% 



%save("C:\Users\Lenovo\Documents\MATLAB\MATLAB1\YXfile.mat","Yin","Xin")


%str=input("Data has been saved");




% load("C:\Users\Lenovo\Documents\MATLAB\MATLAB1\YXfile.mat","Yin","Xin")
% 
% str=input("Data has been loaded");
% 
% 
% 
% 
 Order=7;
 clf;
[BetaCorrH,YCoeffH,XCoeffH,HCorrXX,Zy,Zx] = CalculateHermiteCorrelations0NA04(Xin,Yin,paths,NDims);
YCoeffZ(1:Order+1)= ConvertHermiteSeriesToZSeries01(YCoeffH(1:Order+1),Order);

%HCorrXX(:,:)=0.0;

BetaCorrH

HCorrXX
%CorrXX1st
str=input("I have reached Stone2");


%[cn,dn,bnv,bbnv,ch] = Calculate2DZSeriesNDimCorrected(YCoeffH,BetaCorrH,HCorrXX,NDims,Order);
%cn=bnv;
[cn,ch]=CalculateGaussianCorrelationFromHermiteCorrelationsNDimNew(YCoeffH,XCoeffH,BetaCorrH,HCorrXX,NDims,Order);
% 
% 
% 
% save("C:\Users\Lenovo\Documents\MATLAB\MATLAB1\Zyxfile.mat","Zy","Zx","cn","HCorrXX")
% 
% 
% str=input("Data has been saved");

%load("C:\Users\Lenovo\Documents\MATLAB\MATLAB1\Zyxfile.mat","Zy","Zx","cn","HCorrXX")

%str=input("Data has been loaded");

%size(Zy)
%size(Zx)
%size(cn)
%[HCorrXY,HCorrXX,YCoeffHNorm,XCoeffHNorm]=CalculateHermiteCorrelationCoefficientsForGaussianCorrs(Xin,Yin,paths,NDims);



%[rhoxy,rhoxx] = CalculateGaussianCorrelationFromHermiteCorrelations001(YCoeffH,XCoeffH,HCorrXY,HCorrXX,NDims,Order)
[rhoxy,rhoxx] = CalculateGaussianCorrelationFromHermiteCorrelations001(YCoeffH,XCoeffH,BetaCorrH,HCorrXX,NDims,Order)
%rhoxy
%rhoxx

%str=input("Look at these correlations")

%OptimizeCorrFromHermitesIterative(Zy,Zx,cn,HCorrXX,NDims,paths,Order);

%OptimizeCorrFromHermitesIterative02(Zy,Zx,cn,HCorrXX,NDims,paths,Order);

%OptimizeCorrFromHermitesNewton02(Zy,Zx,cn,HCorrXX,NDims,paths,Order);

% HCorrXX0=HCorrXX;
% for nn1=1:NDims
%     for nn2=1:NDims
%         VariableHIn(2:Order+1)=HCorrXX(nn1,nn2,1:Order);
%         VariableZOut(1:Order+1)= ConvertHermiteSeriesToZSeries01(VariableHIn(1:Order+1),Order);
%         HCorrXXZ(nn1,nn2,1:Order)=VariableZOut(2:Order+1);
%     end
% end



%[cnOptim] = OptimizeCorrFromHermitesNewton(Zy,Zx,cn,HCorrXX,NDims,paths,Order);
[cnOptim] = OptimizeCorrFromHermitesNewton(Zy,Zx,rhoxy,rhoxx(:,:,2:8),NDims,paths,Order);

cn=cnOptim;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5555

[CoeffR,XCoeffH,ZI] = RegressYOnX(Yin,Xin,paths,NDims,Order);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Set1
V10(1)=.65;%.65;%
V10(2)=1.2;
V10(3)=.85;
%V10(4)=1.85;




% % %Set2
 V10(1)=.75;
 V10(2)=1.2;%1.2;
 V10(3)=.7;
% %V10(4)=.65
% % 
% % %Set3
  V10(1)=1.25;
  V10(2)=1.35;%1.2;
  V10(3)=1.25;
%   V10(4)=1.4
% % 
% %   
%   V10(1)=.85;
%   V10(2)=.81;%1.2;
%   V10(3)=.77;
% %  V10(4)=.7

 

%V10(1)=.65;
%V10(2)=.51;%1.2;
%V10(3)=.73;

%V10(1)=1.65;
%V10(2)=1.51;%1.2;
%V10(3111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111)=1.63;

for nn=1:NDims
V(1:paths,nn)=V10(nn);
end
%V(1:paths,1:NDims)=V0(1,1:NDims);

Random1(1:paths/2,1:NDims)=0;

for tt=1:TT2-TT1
    
    Random2(1:paths/2,1:NDims)=randn(size(Random1));
    Random2(paths/2+1:paths,1:NDims)=-Random2(1:paths/2,1:NDims);

    random1=Random2*R(:,1);
    random2=Random2*R(:,2);
    random3=Random2*R(:,3);
    
    Random2(:,1)=random1;
    Random2(:,2)=random2;
    Random2(:,3)=random3;
    
    
    
    for nn=1:NDims
    
    V(1:paths,nn)=V(1:paths,nn)+ ...
        (mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn))*dt + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*Random2(1:paths,nn)*sqrt(dt) + ...
        (mu1(nn).*beta1(nn)*V(1:paths,nn).^(beta1(nn)-1) + mu2(nn).*beta2(nn).*V(1:paths,nn).^(beta2(nn)-1)).* ...
        ((mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn))*dt^2/2 + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*Random2(1:paths,nn)*(1-1/sqrt(3))*dt^1.5) + ...
        .5*(mu1(nn).*beta1(nn).*(beta1(nn)-1).*V(1:paths,nn).^(beta1(nn)-2) + mu2(nn).*beta2(nn).*(beta2(nn)-1).*V(1:paths,nn).^(beta2(nn)-2)).* ...    
        sigma0(nn)^2.*V(1:paths,nn).^(2*gamma(nn)).*dt^2/2 + ...
        sigma0(nn).*gamma(nn).*V(1:paths,nn).^(gamma(nn)-1) .* ...
        ((mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn)).*Random2(1:paths,nn).*1/sqrt(3)*dt^1.5 + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*(Random2(1:paths,nn).^2-1)*dt/2) + ...
        .5*sigma0(nn).*gamma(nn).*(gamma(nn)-1).*V(1:paths,nn).^(gamma(nn)-2) .* ...
        sigma0(nn).^2.*V(1:paths,nn).^(2*gamma(nn)) .*Random2(1:paths,nn).*1/sqrt(3)*dt^1.5;

    V(V(:,nn)<0)=.0001;

    
    %if(tt==TT1)
    %    XXin(1:paths,nn)=V(1:paths,nn);
    %end
        
    if(tt==TT2-TT1)
        
        XXin2(1:paths,nn)=V(1:paths,nn);
    end

    end    
    
end


Yin2(1:paths)=0.0;

for pp=1:paths
for mm=1:NDims
    Yin2(pp)=Yin2(pp)+XXin2(pp,mm);%-V10(mm);
end
end



randn()

str=input("I have reached stone 4");



NoOfBins=200;
MaxCutOff=20;

[XDensity,IndexOut,IndexMax] = MakeDensityFromSimulation_Infiniti_NEW(Yin2,paths,NoOfBins,MaxCutOff );

for qq=1:8
    
    DataMoments(qq)=sum(Yin2(:).^qq)/paths;
    
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

for nn=1:NDims
    [Zx0(nn)] = CalculateZgivenXAndHSeriesCoeffs(V10(nn),XCoeffH(1:Order+1,nn),Order);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55

WW(1,1:Order+1)=EvalHermiteArrayH0(Zx0(1),Order);
for mm=2:NDims
    WW(1,(mm-1)*Order+2:mm*Order+1)=EvalHermiteArrayH1(Zx0(mm),Order);
end    




     YProjected=0.0;
    for mm=1:NDims*Order+1
    YProjected=YProjected+WW(1,mm).*CoeffR(mm,1);    
     
    end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
Zx(1:Order)=0; 

for hh=1:Order
    for nn=1:NDims
        %Zx(hh) = Zx(hh) + bnv(nn+1,hh+1).^1.*Zx0(nn);
        Zx(hh) = Zx(hh) + cn(nn,hh).^1.*Zx0(nn);
    end
end

bh0=YCoeffH(1);

[W] = EvalHermiteArrayH0(Zx0,Order);

for hh=1:Order
    for nn=1:NDims
        bh0 = bh0 + ch(nn+1,hh+1).^1.*W(nn,hh+1);
    end
end
bh00=bh0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


a0=YCoeffZ(1);
a1=YCoeffZ(2);
a2=YCoeffZ(3);
a3=YCoeffZ(4);
a4=YCoeffZ(5);
a5=YCoeffZ(6);
a6=YCoeffZ(7);
a7=YCoeffZ(8);


d1=cn(4,1);
d2=cn(4,2);
d3=cn(4,3);
d4=cn(4,4);
d5=cn(4,5);
d6=cn(4,6);
d7=cn(4,7);


%b0=(a0+a1.* Zx(1) +a2.* Zx(2).^2+a3.* Zx(3).^3+a4* Zx(4).^4+a5 *Zx(5).^5+a6* Zx(6).^6+a7* Zx(7).^7);
b1=(a1* d1+2* a2* Zx(2)* d2 +3* a3* Zx(3).^2* d3+4* a4* Zx(4).^3* d4+5* a5* Zx(5).^4* d5+6* a6* Zx(6).^5* d6+7* a7* Zx(7).^6* d7);
b2=(a2* d2.^2+3* a3* Zx(3)* d3^2 +6* a4* Zx(4).^2* d4^2+10* a5* Zx(5).^3* d5^2+15* a6* Zx(6).^4* d6^2+21* a7* Zx(7).^5* d7^2);
b3=(a3* d3^3+4* a4 *Zx(4)* d4^3 +10* a5* Zx(5).^2* d5^3+20* a6* Zx(6).^3* d6^3+35* a7* Zx(7).^4* d7^3);
b4=(a4* d4^4+5* a5* Zx(5)* d5^4 +15* a6 *Zx(6).^2* d6^4+35* a7* Zx(7).^3 *d7^4);
b5=(a5* d5^5+6* a6* Zx(6)* d6^5 +21* a7 *Zx(7).^2* d7^5);
b6=(a6* d6^6+7* a7* Zx(7)* d7^6 );
b7=(a7* d7^7);


%bh0=b0+b2+3*b4+15*b6;
bh1=b1+3*b3+15*b5+105*b7;
bh2=b2+6*b4+45*b6;
bh3=b3+10*b5+105*b7;
bh4=b4+15*b6;
bh5=b5+21*b7;
bh6=b6;
bh7=b7;


bh(1)=YProjected;
bh(2)=bh1;
bh(3)=bh2;
bh(4)=bh3;
bh(5)=bh4;
bh(6)=bh5;
bh(7)=bh6;
bh(8)=bh7;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

YCoeffz(1:Order+1)= ConvertHermiteSeriesToZSeries01(bh(1:Order+1),Order);


[YMomentsModel]=CalculateMomentsOfZSeries(YCoeffz(1),YCoeffz(2:Order+1),Order,Order+1);
 
YMomentsModel
 
DataMoments
 
str=input("Look at comparison of moments");

[cMu1,YCentralMoments]=ConvertRawMomentsToCentralMoments(YMomentsModel,Order+1);
[cMu2,DataCentralMoments]=ConvertRawMomentsToCentralMoments(DataMoments,Order+1);


YCentralMoments

DataCentralMoments

bh00
YProjected
str=input("Look at central moments");

% Order0=7;
% Order1=13;
% NData=paths;
% [CondMomentYgivX,BetaCorrHxy] = CalculateConditionalMomentYgivXNDims(Xin,Yin,NData,Zx0,NDims,Order0,Order1)
% 
% [VarConditional]=FindSecondConditionalMoment(YCoeffH,BetaCorrH,Zx0)
% 
% 
% BetaCorrHxy
% BetaCorrH
% 
% str=input("Look at conditionalVariance");

clf;

PlotZSeriesDensity(YCoeffz(1),YCoeffz(2:Order+1),'b')


hold on

plot(IndexOut(1:IndexMax),XDensity(1:IndexMax),'r');

%PlotHermiteSeriesDensityAndRvGraph(YCoeffh(1),YCoeffh1(2:yOrder+1),'g')
title("Graph of Density of sum of Three Correlated SDE Variables conditional on starting values of the three SDE random variables.")
%title(sprintf('V00 = %.3f; kappa=%.3f, theta=%.3f,gamma=%.3f,sigma=%.3f;V10=%.3f at TT1=%.3f;TT2=%.3f;Paths=%d ',V00, kappaV, thetaV,gamma,sigma0, V10,(TT1*dt),(TT2*dt),paths));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
legend({'Analytical Density','Numerical Density1'}, ...
    'Location','northeast')

str=input("This is the first version without altering the moments of conditional density")

hold off 
 





end

.
.
.

function [BetaCorrH,YCoeffH,XCoeffH,HCorrXX,Zy,Zx] = CalculateHermiteCorrelations0NA04(Xin,Yin,NData,NDims)

%This program finds required correlations between NDims explanatory
%variables given in Xin that has NDims dimensionality. Yin is variable
%being explained and is one-diemsnional.


Order=7;
NMoments=8;

[ZI] = MomentMatchedStandardNormalDiscretizedDensity(NData);


for nn=1:NDims
[XCoeffH(1:Order+1,nn),Zx(1:NData,nn)] = CalculateHermiteSeriesFromData02(Xin(1:NData,nn),Order,ZI);

[Zx(1:NData,nn)] = CalculateZgivenXAndHSeriesCoeffs(Xin(1:NData,nn),XCoeffH(1:Order+1,nn),Order);

end


[YCoeffH,Zy] = CalculateHermiteSeriesFromData02(Yin,Order,ZI); 
[Zy] = CalculateZgivenXAndHSeriesCoeffs(Yin,YCoeffH,Order);


for nn1=1:NDims
    for nn2=1:NDims
        if(nn1==nn2)
            HCorrXX(nn1,nn2,1:Order)=1.0;
        else
            HCorrXX(nn1,nn2,1:Order) = CalculateCorrelationBivariateHermiteCH0(Zx(1:NData,nn1),Zx(1:NData,nn2),NData,Order);
        end
    end
end
for nn1=1:NDims
    HCorrXX(nn1,nn1,1:Order)=1.0;
end


for hh=1:Order
    for nn1=1:NDims
        HCorrXX
        ss2=0;
        for nn2=1:NDims
            if(nn2~=nn1)
                ss2=ss2+1;
                ss3=0;
                for nn3=1:NDims
                    if(nn3~=nn1)
                        ss3=ss3+1;
                            HCorrXXmn(ss2,ss3)=HCorrXX(nn2,nn3,hh);
                    end
                    HCorrXmXn(ss2,1)=HCorrXX(nn1,nn2,hh);
                end
            end
        end
        BetaCorrXX(1:NDims-1,nn1,hh)=HCorrXXmn(:,:)\HCorrXmXn(1:NDims-1,1);
            
    end
        
    BetaCorrXX
        
    nn1
       
        str=input("Look at important variables");
    
end
BetaCorrXX            
str=input("look at the BetacorrXX");       
for nn1=1:NDims
    for nn2=1:NDims-1
        if(nn2<nn1)
    BetaCorrXX0(nn2,nn1,1:Order)=BetaCorrXX(nn2,nn1,1:Order);
        elseif(nn2==nn1)
               % BetaCorrXX0(nn2,nn1,1:Order)=1;
            
            BetaCorrXX0(nn2+1:NDims,nn1,1:Order)=BetaCorrXX(nn2:NDims-1,nn1,1:Order);
        end
    end
    BetaCorrXX0(nn1,nn1,1:Order)=1;
end
BetaCorrXX
BetaCorrXX0
% 
% % % % str=input("Look at BetaCorrXX pairs");
% % % % for hh=1:Order
% % % %     for nn1=1:NDims
% % % %         for nn2=1:NDims-1
% % % %             BetaCorrXXn(nn1,nn2,hh)=sign(BetaCorrXXn(nn1,nn2,hh)).*(abs(BetaCorrXXn(nn1,nn2,hh))).^(1/hh);
% % % %         end
% % % %     end
% % % % end



for nn1=1:NDims
    HCorrXY(nn1,1:Order) = CalculateCorrelationBivariateHermiteCH0(Zx(1:NData,nn1),Zy(1:NData),NData,Order);
end

for hh=1:Order
    [BCorrXY(:,hh),Is(:,hh)]=sort(HCorrXY(:,hh),'descend');
end

for hh=1:Order
    CorrXX(1:NDims,1:NDims)=HCorrXX(Is(1:NDims,hh),Is(1:NDims,hh),hh);
    CorrXY(1:NDims,1)=HCorrXY(Is(1:NDims,hh),hh);
    BetaCorr0(1:NDims)=CorrXX\CorrXY;
    BetaCorrH(Is(1:NDims,hh),hh)=BetaCorr0(1:NDims);
end



HCorrXX=BetaCorrXX0;
HCorrXX
str=input("Look at HCOrrXX");

% for hh=1:Order
%     CorrXX(1:NDims,1:NDims)=HCorrXX(1:NDims,1:NDims,hh);
%     CorrXY(1:NDims,1)=HCorrXY(1:NDims,hh);
%     BetaCorrH(1:NDims,hh)=CorrXX\CorrXY;
% end




%CorrXX1st(1:NDims,1:NDims)=HCorrXX(1:NDims,1:NDims,1)


BetaCorrH(:,1)'*CorrXX

CorrXY

HCorrXX

HCorrXY

BetaCorrH

str=input("Look at BetaCorrH");
end

.
.
.

function [cn,ch] = CalculateGaussianCorrelationFromHermiteCorrelationsNDimNew(YCoeffH,XCoeffH,BetaCorrHxy,CorrHxx,NDims,Order)





for hh=1:Order
    for nn=1:NDims
        ch(nn+1,hh+1)=BetaCorrHxy(nn,hh).*YCoeffH(hh+1);%make sure indices on all hermite-correlations match.check.
        %%%% nn+1 since 1st index would go to y-tilde.
    end
end


ch
YCoeffH

str=input("Look at above---00-1-1-1-1-");

[YCoeffZ(1:Order+1)] = ConvertHermiteSeriesToZSeries01(YCoeffH(1:Order+1),Order);

for nn=1:NDims
    [XCoeffZ(1:Order+1,nn)] = ConvertHermiteSeriesToZSeries01(XCoeffH(1:Order+1,nn),Order);
end




YCoeffH

YCoeffZ

str=input("Look at the Hermite and Z-series coeffs");

% YCoeffZA=YCoeffZ;
% YCoeffZ=YCoeffZA;
% if(abs(YCoeffZ(4))>abs(YCoeffZ(2)))
%    YCoeffZ(4)=sign(YCoeffZ(4)).*abs(YCoeffZ(2));
% end
% if(abs(YCoeffZ(6))>abs(YCoeffZ(4)))
%    YCoeffZ(6)=sign(YCoeffZ(6)).*abs(YCoeffZ(4));
% end
% if(abs(YCoeffZ(8))>abs(YCoeffZ(6)))
%    YCoeffZ(6)=sign(YCoeffZ(6)).*abs(YCoeffZ(8));
% end
% 
% 
% if(abs(YCoeffZ(5))>abs(YCoeffZ(3)))
%    YCoeffZ(5)=sign(YCoeffZ(5)).*abs(YCoeffZ(3));
% end
% if(abs(YCoeffZ(7))>abs(YCoeffZ(5)))
%    YCoeffZ(5)=sign(YCoeffZ(5)).*abs(YCoeffZ(7));
% end
% 




for nn=1:NDims
[c(nn+1,1:Order+1)] = ConvertHermiteSeriesToZSeries01(ch(nn+1,1:Order+1),Order);
end


for nn=1:NDims
    rhoxy(nn+1,8)=sign(BetaCorrHxy(nn,7)).^1.*(abs(BetaCorrHxy(nn,7))).^(1/7);
    rhoxy(nn+1,7)=sign(BetaCorrHxy(nn,6)).^1.*(abs(BetaCorrHxy(nn,6))).^(1/6);
    rhoxy(nn+1,6)=sign(BetaCorrHxy(nn,5)).^1.*(abs(   (YCoeffH(6).*BetaCorrHxy(nn,5)-21*YCoeffZ(8).*rhoxy(nn+1,8).^(5))./(YCoeffH(6)-21*YCoeffZ(8)))).^(1/5);
    rhoxy(nn+1,5)=sign(BetaCorrHxy(nn,4)).^1.*(abs(  (YCoeffH(5).*BetaCorrHxy(nn,4)-15*YCoeffZ(7).*rhoxy(nn+1,7).^(4))./(YCoeffH(5)-15*YCoeffZ(7)))).^(1/4);
    rhoxy(nn+1,4)=sign(BetaCorrHxy(nn,3)).^1.*(abs(  (YCoeffH(4).*BetaCorrHxy(nn,3)-10*YCoeffZ(6).*rhoxy(nn+1,6).^(3)-105*YCoeffZ(8).*rhoxy(nn+1,8).^(3))...
        ./(YCoeffH(4)-10*YCoeffZ(6)-105*YCoeffZ(8)))).^(1/3);
    rhoxy(nn+1,3)=sign(BetaCorrHxy(nn,2)).^1.*(abs( (YCoeffH(3).*BetaCorrHxy(nn,2)-6*YCoeffZ(5).*rhoxy(nn+1,5).^(2)-45*YCoeffZ(7).*rhoxy(nn+1,7).^(2)) ...
        ./(YCoeffH(3)-6*YCoeffZ(5)-45*YCoeffZ(7)))).^(1/2);
    rhoxy(nn+1,2)=sign(BetaCorrHxy(nn,1)).^1.*(abs(  (YCoeffH(2).*BetaCorrHxy(nn,1)-3*YCoeffZ(4).*rhoxy(nn+1,4).^(1)-15*YCoeffZ(6).*rhoxy(nn+1,6).^(1) ...
        -105*YCoeffZ(8).*rhoxy(nn+1,8).^(1))...
        ./(YCoeffH(2)-3*YCoeffZ(4)-15*YCoeffZ(6)-105*YCoeffZ(8)))).^(1/1);
    
%    rhoxy(nn,2+1)=sign(4/9*BetaCorrHxy(nn,2)+1/36).*abs((4/9*BetaCorrHxy(nn,2)+1/36))^(.5);
% rhoxy(nn,3+1)=sign((36/225*BetaCorrHxy(nn,3)-3/225*BetaCorrHxy(nn,1))).*(abs((36/225*BetaCorrHxy(nn,3)-3/225*BetaCorrHxy(nn,1))))^(1/3);
% rhoxy(nn,4+1)=sign((144*4/(105)^2*BetaCorrHxy(nn,4)-6/105^2*(4*BetaCorrHxy(nn,2)+.25)-18/105^2)).*(abs(((144*4/(105)^2*BetaCorrHxy(nn,4)-6/105^2*(4*BetaCorrHxy(nn,2)+.25)-18/105^2))))^(1/4);
% rhoxy(nn,1+1)=BetaCorrHxy(nn,1);
end


% 
for nn=1:NDims
    for mm=1:NDims
        if(nn==mm)
            rhoxx(nn,mm,2)=1;
            rhoxx(nn,mm,3)=1;
            rhoxx(nn,mm,4)=1;
            rhoxx(nn,mm,5)=1;
            rhoxx(nn,mm,6)=1;
            rhoxx(nn,mm,7)=1;
            rhoxx(nn,mm,8)=1;
            
            
            
        elseif(nn<mm)
            rhoxx(nn,mm,8)=sign(CorrHxx(nn,mm,7)).*(abs(CorrHxx(nn,mm,7))).^(1/7);
            rhoxx(nn,mm,7)=sign(CorrHxx(nn,mm,6)).*(abs(CorrHxx(nn,mm,6))).^(1/6);
            rhoxx(nn,mm,6)=sign(CorrHxx(nn,mm,5)).*(abs( (XCoeffH(6,nn).*CorrHxx(nn,mm,5)-21*XCoeffZ(8,nn).*rhoxx(nn,mm,8).^(5))./(XCoeffH(6,nn)-21*XCoeffZ(8,nn)))).^(1/5);
            rhoxx(nn,mm,5)=sign(CorrHxx(nn,mm,4)).*(abs( (XCoeffH(5,nn).*CorrHxx(nn,mm,4)-15*XCoeffZ(7,nn).*rhoxx(nn,mm,7).^(4))./(XCoeffH(5,nn)-15*XCoeffZ(7,nn)))).^(1/4);
            rhoxx(nn,mm,4)=sign(CorrHxx(nn,mm,3)).*(abs( (XCoeffH(4,nn).*CorrHxx(nn,mm,3)-10*XCoeffZ(6,nn).*rhoxx(nn,mm,6).^(3)-105*XCoeffZ(8,nn).*rhoxx(nn,mm,8).^(3))...
                ./(XCoeffH(4,nn)-10*XCoeffZ(6,nn)-105*XCoeffZ(8,nn)))).^(1/3);
            rhoxx(nn,mm,3)=sign(CorrHxx(nn,mm,2)).*(abs( (XCoeffH(3,nn).*CorrHxx(nn,mm,2)-6*XCoeffZ(5,nn).*rhoxx(nn,mm,5).^(2)-45*XCoeffZ(7,nn).*rhoxx(nn,mm,7).^(2)) ...
                ./(XCoeffH(3,nn)-6*XCoeffZ(5,nn)-45*XCoeffZ(7,nn)))).^(1/2);
            rhoxx(nn,mm,2)=sign(CorrHxx(nn,mm,1)).*(abs( (XCoeffH(2,nn).*CorrHxx(nn,mm,1)-3*XCoeffZ(4,nn).*rhoxx(nn,mm,4).^(1)-15*XCoeffZ(6,nn).*rhoxx(nn,mm,6).^(1) ...
                -105*XCoeffZ(8,nn).*rhoxx(nn,mm,8).^(1))...
                ./(XCoeffH(2,nn)-3*XCoeffZ(4,nn)-15*XCoeffZ(6,nn)-105*XCoeffZ(8,nn)))).^(1/1);
        

           % rhoxx1(nn,mm,6)=sign(CorrHxx(nn,mm,5)).*(abs( (CorrHxx(nn,mm,5)-21.*rhoxx(nn,mm,8).^(5)))).^(1/5);
           % rhoxx1(nn,mm,5)=(abs( CorrHxx(nn,mm,4)-15.*rhoxx(nn,mm,7).^(4))).^(1/4);
           %rhoxx1(mm,nn,6)=rhoxx1(nn,mm,6);
           % rhoxx1(mm,nn,5)=rhoxx1(nn,mm,5);
            
           % rhoxx(mm,nn,8)=rhoxx(nn,mm,8);
           % rhoxx(mm,nn,7)=rhoxx(nn,mm,7);
rhoxx(mm,nn,6)=rhoxx(nn,mm,6);
rhoxx(mm,nn,5)=rhoxx(nn,mm,5);
rhoxx(mm,nn,4)=rhoxx(nn,mm,4);
rhoxx(mm,nn,3)=rhoxx(nn,mm,3);
rhoxx(mm,nn,2)=rhoxx(nn,mm,2);
        end
    end
end
% % 
% 
% rhoxx(:,:,6)
% rhoxx1(:,:,6)
% rhoxx(:,:,5)
% rhoxx1(:,:,5)
% 
% str=input("Look at rhoxx and rhoxx1");

% 
% % 
% % 
% % 
% % %rhoxx(1:NDims,1:NDims,2:8)=CorrHxx(1:NDims,1:NDims,1:7);
%  for hh=1:7
% % %rhoxx(1:NDims,1:NDims,hh+1)=rhoxx(1:NDims,1:NDims,2);
%  rhoxx(1:NDims,1:NDims,hh+1)=(sign(CorrHxx(1:NDims,1:NDims,hh))).^(1).*(abs(CorrHxx(1:NDims,1:NDims,hh))).^(1/hh);
%  for nn=1:NDims
%      for mm=1:NDims
%          if(nn~=mm)
%  rhoxx1(nn,mm,2+1)=-sign(4/9*CorrHxx(nn,mm,2)+1/36).^0.*abs((4/9*CorrHxx(nn,mm,2)+1/36))^(.5);
%  rhoxx1(nn,mm,3+1)=-sign((36/225*CorrHxx(nn,mm,3)-3/225*CorrHxx(nn,mm,1))).^0.*(abs((36/225*CorrHxx(nn,mm,3)-3/225*CorrHxx(nn,mm,1))))^(1/3);
%  rhoxx1(nn,mm,4+1)=-sign((144*4/(105)^2*CorrHxx(nn,mm,4)-6/105^2*(4*CorrHxx(nn,mm,2)+.25)-18/105^2)).^0.*(abs(((144*4/(105)^2*CorrHxx(nn,mm,4)-6/105^2*(4*CorrHxx(nn,mm,2)+.25)-18/105^2))))^(1/4);
%  %rhoxx1(nn,mm,2+1)=sign(4/9*CorrHxx(nn,mm,2)+1/36).^(2).*abs((4/9*CorrHxx(nn,mm,2)+1/36))^(.5);
%  %rhoxx1(nn,mm,3+1)=sign((36/225*CorrHxx(nn,mm,3)-3/225*CorrHxx(nn,mm,1))).^(3).*(abs((36/225*CorrHxx(nn,mm,3)-3/225*CorrHxx(nn,mm,1))))^(1/3);
%  %rhoxx1(nn,mm,4+1)=sign((144*4/(105)^2*CorrHxx(nn,mm,4)-6/105^2*(4*CorrHxx(nn,mm,2)+.25)-18/105^2)).^(4).*(abs(((144*4/(105)^2*CorrHxx(nn,mm,4)-6/105^2*(4*CorrHxx(nn,mm,2)+.25)-18/105^2))))^(1/4);
%  
%          else
%              rhoxx1(nn,mm,2+1)=1.0;
%               rhoxx1(nn,mm,3+1)=1.0;
%                rhoxx1(nn,mm,4+1)=1.0;
%          end
%      end
%  end
%  end
%  
%  rhoxx(1:NDims,1:NDims,3)=rhoxx1(1:NDims,1:NDims,3);
%  rhoxx(1:NDims,1:NDims,4:5)=rhoxx1(1:NDims,1:NDims,4:5);
%  
 for hh=1:7
 %rhoxx(1:NDims,1:NDims,hh+1)=sign(CorrHxx(1:NDims,1:NDims,hh)).*abs(sign(CorrHxx(1:NDims,1:NDims,hh))).^(1/hh);
 end
 %rhoxx
 %rhoxx1
% 
% 
 str=input("Look at rhoxx");

% 
% CorrHxx0(1:NDims,1:NDims,2:Order+1)=CorrHxx(1:NDims,1:NDims,1:Order);
% 
% for nn=1:NDims
%     for mm=1:NDims
% 
%         if(nn==mm)
%             rhoxx(nn,mm,2)=1;
%             rhoxx(nn,mm,3)=1;
%             rhoxx(nn,mm,4)=1;
%             rhoxx(nn,mm,5)=1;
%             rhoxx(nn,mm,6)=1;
%             rhoxx(nn,mm,7)=1;
%             rhoxx(nn,mm,8)=1;
%             
%         else
% [rrhoxx(nn+1,mm+1,1:Order+1)] = ConvertHermiteSeriesToZSeries01(CorrHxx0(nn,mm,1:Order+1),Order);
% end
% end
% end
% 
% for hh=1:7
% %rhoxx(1:NDims,1:NDims,hh+1)=rhoxx(1:NDims,1:NDims,2);
% rhoxx(1:NDims,1:NDims,hh+1)=sign(rrhoxx(1:NDims,1:NDims,hh+1)).*(abs(rrhoxx(1:NDims,1:NDims,hh+1))).^(1/hh);
% end







%rhoxx(:,:)=0.0;

%(YCoeffH(6).*BetaCorrHxy(nn,5)-21*YCoeffZ(8).*rhoxy(nn+1,8).^(5/7))
%(abs((abs(YCoeffH(6)).*BetaCorrHxy(nn,5)-21*abs(YCoeffZ(8)).*rhoxy(nn+1,8).^(5/7))./abs(YCoeffH(6)-21*(YCoeffZ(8))))).^1/5
%(YCoeffH(6)-21*YCoeffZ(8))
%(YCoeffZ(6)-21*YCoeffZ(8))

%(abs(abs(YCoeffH(6)).*BetaCorrHxy(nn,5)-21*abs(YCoeffZ(8)).*rhoxy(nn+1,8).^(5/7))).^(1/5)



%YCoeffH
%YCoeffZ
%rhoxy

%str=input("Look at rhoxy");



cn=rhoxy;





%for nn=1:NDims
%    rho_z(nn+1,7)=c(nn+1,8)./YCoeffZ(8);
%    cn(nn+1,8)=(abs(rho_z(nn+1,7))).^(1/7);
%    cn(nn+1,8)=cn(nn+1,8).*sign(rho_z(nn+1,7));
%end

cnSqrSum=0;
for nn=1:NDims
    %ccn0(nn+1,8)=cn(nn+1,8).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,8).*cn(nn+1,8);
    for mm=1:NDims
        if(nn~=mm)
            %ccn0(nn+1,8)=ccn0(nn+1,8)+2*cn(nn+1,8).*cn(mm+1,8).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,8).*cn(mm+1,8).*rhoxx(nn,mm,8);
        end
    end
    %ccn(nn+1,8)=sqrt(1-ccn0(nn+1,8));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,8)=cn(2:NDims+1,8)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end

cn(1,8)=sqrt(1-cnSqrSum);
%ccn(1,8)=sqrt(1-cn(1,8).^2);
%Cn(1,8)=cn(1,8).*YCoeffZ(8);

cn
%ccn
str=input("Compare the Eighth Element From NDim---2");


% ccn00(2,8)=cn(2,8).^2+2*cn(2,8).*cn(3,8).*CorrXX1st(1,2)+2*cn(2,8).*cn(4,8).*CorrXX1st(1,3)+ ...
%     cn(3,8).^2+2*cn(3,8).*cn(4,8).*CorrXX1st(2,3)+cn(4,8).^2
% 
% ccn00(2,8)=sqrt(1-(cn(2,8).^2+2*cn(2,8).*cn(3,8).*CorrXX1st(1,2)+2*cn(2,8).*cn(4,8).*CorrXX1st(1,3)))
% 
% cnSqrSum
% 
% ccn0(2,8)
% 
% ccn(1,8)
% 
% ccn(2,8)
% 
% str=input("Look at the variables Inside");
% 







cnSqrSum=0;
%for nn=1:NDims
%    rho_z(nn+1,6)=c(nn+1,7)./YCoeffZ(7);
%    cn(nn+1,7)=(abs(rho_z(nn+1,6))).^(1/6);
%    cn(nn+1,7)=cn(nn+1,7).*sign(rho_z(nn+1,6));
%    %ccn(nn+1,7)=sqrt(1-cn(nn+1,7).^2);
%    %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
%end

cnSqrSum=0;
for nn=1:NDims
    %ccn0(nn+1,7)=cn(nn+1,7).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,7).*cn(nn+1,7);
    for mm=1:NDims
        if(nn~=mm)
         %   ccn0(nn+1,7)=ccn0(nn+1,7)+2*cn(nn+1,7).*cn(mm+1,7).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,7).*cn(mm+1,7).*rhoxx(nn,mm,7);
        end
    end
    %ccn(nn+1,7)=sqrt(1-ccn0(nn+1,7));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,7)=cn(2:NDims+1,7)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end

cn(1,7)=sqrt(1-cnSqrSum);
%ccn(1,7)=sqrt(1-cn(1,7).^2);
%Cn(1,7)=cn(1,7).*YCoeffZ(7);


cn
%ccn
str=input("Compare the Seventh Element From NDim---2");


%cnSqrSum=0;
%for nn=1:NDims
%    CC(3,6)=YCoeffZ(8).*(21*cn(nn+1,8).^5.*ccn(nn+1,8).^2);
%    rho_z(nn+1,5)=(c(nn+1,6)-CC(3,6))./YCoeffZ(6);
%    if(abs(rho_z(nn+1,5))>.995)
%        rho_z(nn+1,5)=sign(rho_z(nn+1,5)).*.995;
%    end
%    cn(nn+1,6)=(abs(rho_z(nn+1,5))).^(1/5);
%    cn(nn+1,6)=cn(nn+1,6).*sign(rho_z(nn+1,5));
%%%    nn
%%%    c(nn+1,6)
%%%    CC(3,6)
%%%    YCoeffZ(8)
%%%    YCoeffZ(6)
%%%    rho_z(nn+1,5)
%%%    cn(nn+1,6)
%%%    
%%%    str=input("Look inside the correlations module loop");
    
    
    
    
    
%    ccn(nn+1,6)=sqrt(1-cn(nn+1,6).^2);
%    %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
%end

cnSqrSum=0;
for nn=1:NDims
%    ccn0(nn+1,6)=cn(nn+1,6).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,6).*cn(nn+1,6);
    for mm=1:NDims
        if(nn~=mm)
          %  ccn0(nn+1,6)=ccn0(nn+1,6)+2*cn(nn+1,6).*cn(mm+1,6).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,6).*cn(mm+1,6).*rhoxx(nn,mm,6);
        end
    end
    %ccn(nn+1,6)=sqrt(1-ccn0(nn+1,6));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,6)=cn(2:NDims+1,6)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end



cn(1,6)=sqrt(1-cnSqrSum);
%ccn(1,6)=sqrt(1-cn(1,6).^2);
%Cn(1,6)=cn(1,6).*YCoeffZ(6);

cn
%ccn
str=input("Compare the Sixth Element From NDim---2");



% cnSqrSum=0;
% for nn=1:NDims
%     CC(3,5)=YCoeffZ(7).*(15*cn(nn+1,7).^4.*ccn(nn+1,7).^2);
%     rho_z(nn+1,4)=(c(nn+1,5)-CC(3,5))./YCoeffZ(5);
%     cn(nn+1,5)=(abs(rho_z(nn+1,4))).^(1/4);
%     cn(nn+1,5)=cn(nn+1,5).*sign(rho_z(nn+1,4));
%     %ccn(nn+1,5)=sqrt(1-cn(nn+1,5).^2);
%     %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
% end

cnSqrSum=0;
for nn=1:NDims
    %ccn0(nn+1,5)=cn(nn+1,5).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,5).*cn(nn+1,5);
    for mm=1:NDims
        if(nn~=mm)
           % ccn0(nn+1,5)=ccn0(nn+1,5)+2*cn(nn+1,5).*cn(mm+1,5).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,5).*cn(mm+1,5).*rhoxx(nn,mm,5);
        end
    end
    %ccn(nn+1,5)=sqrt(1-ccn0(nn+1,5));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,5)=cn(2:NDims+1,5)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end


cn(1,5)=sqrt(1-cnSqrSum);
%ccn(1,5)=sqrt(1-cn(1,5).^2);
%Cn(1,5)=cn(1,5).*YCoeffZ(5);


cn
%ccn
str=input("Compare the Fifth Element From NDim---2");




%CC(3,4)=YCoeffZ(6).*(10*cn(6).^3.*dn(6).^2*zeta^2);
%CC(5,4)=YCoeffZ(8).*(35*cn(8).^3.*dn(8).^4.*zeta^4);
% cnSqrSum=0;
% for nn=1:NDims
%     CC(5,4)=YCoeffZ(8).*(35*cn(nn+1,8).^3.*ccn(nn+1,8).^4);
%     CC(3,4)=YCoeffZ(6).*(10*cn(nn+1,6).^3.*ccn(nn+1,6).^2);
%     rho_z(nn+1,3)=(c(nn+1,4)-3*CC(5,4)-CC(3,4))./YCoeffZ(4);
%     cn(nn+1,4)=(abs(rho_z(nn+1,3))).^(1/3);
%     cn(nn+1,4)=cn(nn+1,4).*sign(rho_z(nn+1,3));
%     %ccn(nn+1,4)=sqrt(1-cn(nn+1,4).^2);
%     %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
% end

cnSqrSum=0;
for nn=1:NDims
%    ccn0(nn+1,4)=cn(nn+1,4).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,4).*cn(nn+1,4);
    for mm=1:NDims
        if(nn~=mm)
          %  ccn0(nn+1,4)=ccn0(nn+1,4)+2*cn(nn+1,4).*cn(mm+1,4).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,4).*cn(mm+1,4).*rhoxx(nn,mm,4);
        end
    end
    %ccn(nn+1,4)=sqrt(1-ccn0(nn+1,4));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,4)=cn(2:NDims+1,4)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end

cn(1,4)=sqrt(1-cnSqrSum);
%ccn(1,4)=sqrt(1-cn(1,4).^2);
%Cn(1,4)=cn(1,4).*YCoeffZ(4);

cn
%ccn
str=input("Compare the Fourth Element From NDim---2");



%CC(3,3)=YCoeffZ(5).*(6* cn(5).^2.* dn(5).^2);
%CC(5,3)=YCoeffZ(7).*(15*cn(7).^2.*dn(7).^4);
% cnSqrSum=0;
% for nn=1:NDims
%     CC(3,3)=YCoeffZ(5).*(6*cn(nn+1,5).^2.*ccn(nn+1,5).^2);
%     CC(5,3)=YCoeffZ(7).*(15*cn(nn+1,7).^2.*ccn(nn+1,7).^4);
%     rho_z(nn+1,2)=(c(nn+1,3)-3*CC(5,3)-CC(3,3))./YCoeffZ(3);
%     cn(nn+1,3)=(abs(rho_z(nn+1,2))).^(1/2);
%     cn(nn+1,3)=cn(nn+1,3).*sign(rho_z(nn+1,2));
%     %ccn(nn+1,3)=sqrt(1-cn(nn+1,3).^2);
%     %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
% end

cnSqrSum=0;
for nn=1:NDims
    %ccn0(nn+1,3)=cn(nn+1,3).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,3).*cn(nn+1,3);
    for mm=1:NDims
        if(nn~=mm)
           % ccn0(nn+1,3)=ccn0(nn+1,3)+2*cn(nn+1,3).*cn(mm+1,3).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,3).*cn(mm+1,3).*rhoxx(nn,mm,3);
        end
    end
    %ccn(nn+1,3)=sqrt(1-ccn0(nn+1,3));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,3)=cn(2:NDims+1,3)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end

cn(1,3)=sqrt(1-cnSqrSum);
%ccn(1,3)=sqrt(1-cn(1,3).^2);
%Cn(1,3)=cn(1,3).*YCoeffZ(3);

cn
%ccn
str=input("Compare the Third Element From NDim---2");



%CC(7,2)=YCoeffZ(8).*(7*cn(8).*dn(8).^6*zeta^6);
%CC(5,2)=YCoeffZ(6).*(5*cn(6).*dn(6).^4.*zeta^4);
%CC(3,2)=YCoeffZ(4).* (3.* cn(4).* dn(4).^2.* zeta^2);
% cnSqrSum=0;
% for nn=1:NDims
%     CC(4,2)=YCoeffZ(4).*(3*cn(nn+1,4).*ccn(nn+1,4).^2);
%     CC(6,2)=YCoeffZ(6).*(5*cn(nn+1,5).*ccn(nn+1,6).^4);
%     CC(8,2)=YCoeffZ(8).*(7*cn(nn+1,8).*ccn(nn+1,8).^6);
%     rho_z(nn+1,1)=(c(nn+1,2)-15*CC(8,2)-3*CC(6,2)-CC(4,2))./YCoeffZ(2);
%     cn(nn+1,2)=(abs(rho_z(nn+1,1)));
%     %ccn(nn+1,2)=sqrt(1-cn(nn+1,2).^2);
%     %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
% end

cnSqrSum=0;
for nn=1:NDims
%    ccn0(nn+1,2)=cn(nn+1,2).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,2).*cn(nn+1,2);
    for mm=1:NDims
        if(nn~=mm)
        %    ccn0(nn+1,2)=ccn0(nn+1,2)+2*cn(nn+1,2).*cn(mm+1,2).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,2).*cn(mm+1,2).*rhoxx(nn,mm,2);
        end
    end
%    ccn(nn+1,2)=sqrt(1-ccn0(nn+1,2));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,2)=cn(2:NDims+1,2)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end
cn(1,2)=sqrt(1-cnSqrSum);
%ccn(1,2)=sqrt(1-cn(1,2).^2);
%Cn(1,2)=cn(1,2).*YCoeffZ(2);

cn
%ccn
str=input("Compare the Second Element From NDim---2");


% 
% %CC(3,1)=YCoeffZ(2).*(ccn(nn+1,2).^2);
% %CC(5,1)=YCoeffZ(4).*(ccn(nn+1,4).^4);
% %CC(7,1)=YCoeffZ(6).*(ccn(nn+1,6).^6);
% CC(3,1)=YCoeffZ(2).*(cn(1,2).^2);
% CC(5,1)=YCoeffZ(4).*(cn(1,4).^4);
% CC(7,1)=YCoeffZ(6).*(cn(1,6).^6);
% 
% 
% 
% 
% 
% cn(1,1)=(YCoeffZ(1)-CC(3,1)-3*CC(5,1)-15*CC(7,1))./YCoeffZ(1);
% 
% 
% cn
% ccn
% str=input("Look at coefficients cn----PleASE lOOK");
% 
% 
 dn(1:Order+1)=cn(1,1:Order+1);
% bn(1:Order+1)=ccn(1,1:Order+1);
% 
% 
% 
% 

end

.
.
.

function [rhoxy,rhoxx] = CalculateGaussianCorrelationFromHermiteCorrelations001(YCoeffH,XCoeffH,BetaCorrHxy,CorrHxx,NDims,Order)





for hh=1:Order
    for nn=1:NDims
        ch(nn+1,hh+1)=BetaCorrHxy(nn,hh).*YCoeffH(hh+1);%make sure indices on all hermite-correlations match.check.
        %%%% nn+1 since 1st index would go to y-tilde.
    end
end


ch
YCoeffH

str=input("Look at above---00-1-1-1-1-");

[YCoeffZ(1:Order+1)] = ConvertHermiteSeriesToZSeries01(YCoeffH(1:Order+1),Order);

for nn=1:NDims
    [XCoeffZ(1:Order+1,nn)] = ConvertHermiteSeriesToZSeries01(XCoeffH(1:Order+1,nn),Order);
end




YCoeffH

YCoeffZ

str=input("Look at the Hermite and Z-series coeffs-----here---here");

% YCoeffZA=YCoeffZ;
% YCoeffZ=YCoeffZA;
% if(abs(YCoeffZ(4))>abs(YCoeffZ(2)))
%    YCoeffZ(4)=sign(YCoeffZ(4)).*abs(YCoeffZ(2));
% end
% if(abs(YCoeffZ(6))>abs(YCoeffZ(4)))
%    YCoeffZ(6)=sign(YCoeffZ(6)).*abs(YCoeffZ(4));
% end
% if(abs(YCoeffZ(8))>abs(YCoeffZ(6)))
%    YCoeffZ(6)=sign(YCoeffZ(6)).*abs(YCoeffZ(8));
% end
% 
% 
% if(abs(YCoeffZ(5))>abs(YCoeffZ(3)))
%    YCoeffZ(5)=sign(YCoeffZ(5)).*abs(YCoeffZ(3));
% end
% if(abs(YCoeffZ(7))>abs(YCoeffZ(5)))
%    YCoeffZ(5)=sign(YCoeffZ(5)).*abs(YCoeffZ(7));
% end
% 




for nn=1:NDims
[c(nn+1,1:Order+1)] = ConvertHermiteSeriesToZSeries01(ch(nn+1,1:Order+1),Order);
end


for nn=1:NDims
    rhoxy(nn+1,8)=sign(BetaCorrHxy(nn,7)).^0.*(abs(BetaCorrHxy(nn,7))).^(1/7);
    rhoxy(nn+1,7)=sign(BetaCorrHxy(nn,6)).^0.*(abs(BetaCorrHxy(nn,6))).^(1/6);
    rhoxy(nn+1,6)=sign(BetaCorrHxy(nn,5)).^0.*(abs(   (YCoeffH(6).*BetaCorrHxy(nn,5)-21*YCoeffZ(8).*rhoxy(nn+1,8).^(5))./(YCoeffH(6)-21*YCoeffZ(8)))).^(1/5);
    rhoxy(nn+1,5)=sign(BetaCorrHxy(nn,4)).^0.*(abs(  (YCoeffH(5).*BetaCorrHxy(nn,4)-15*YCoeffZ(7).*rhoxy(nn+1,7).^(4))./(YCoeffH(5)-15*YCoeffZ(7)))).^(1/4);
    rhoxy(nn+1,4)=sign(BetaCorrHxy(nn,3)).^0.*(abs(  (YCoeffH(4).*BetaCorrHxy(nn,3)-10*YCoeffZ(6).*rhoxy(nn+1,6).^(3)-105*YCoeffZ(8).*rhoxy(nn+1,8).^(3))...
        ./(YCoeffH(4)-10*YCoeffZ(6)-105*YCoeffZ(8)))).^(1/3);
    rhoxy(nn+1,3)=sign(BetaCorrHxy(nn,2)).^0.*(abs( (YCoeffH(3).*BetaCorrHxy(nn,2)-6*YCoeffZ(5).*rhoxy(nn+1,5).^(2)-45*YCoeffZ(7).*rhoxy(nn+1,7).^(2)) ...
        ./(YCoeffH(3)-6*YCoeffZ(5)-45*YCoeffZ(7)))).^(1/2);
    rhoxy(nn+1,2)=sign(BetaCorrHxy(nn,1)).^0.*((  (YCoeffH(2).*BetaCorrHxy(nn,1)-3*YCoeffZ(4).*rhoxy(nn+1,4).^(1)-15*YCoeffZ(6).*rhoxy(nn+1,6).^(1) ...
        -105*YCoeffZ(8).*rhoxy(nn+1,8).^(1))...
        ./(YCoeffH(2)-3*YCoeffZ(4)-15*YCoeffZ(6)-105*YCoeffZ(8)))).^(1/1);
    
%    rhoxy(nn,2+1)=sign(4/9*BetaCorrHxy(nn,2)+1/36).*abs((4/9*BetaCorrHxy(nn,2)+1/36))^(.5);
% rhoxy(nn,3+1)=sign((36/225*BetaCorrHxy(nn,3)-3/225*BetaCorrHxy(nn,1))).*(abs((36/225*BetaCorrHxy(nn,3)-3/225*BetaCorrHxy(nn,1))))^(1/3);
% rhoxy(nn,4+1)=sign((144*4/(105)^2*BetaCorrHxy(nn,4)-6/105^2*(4*BetaCorrHxy(nn,2)+.25)-18/105^2)).*(abs(((144*4/(105)^2*BetaCorrHxy(nn,4)-6/105^2*(4*BetaCorrHxy(nn,2)+.25)-18/105^2))))^(1/4);
% rhoxy(nn,1+1)=BetaCorrHxy(nn,1);
end


% 
for nn=1:NDims
    for mm=1:NDims
        if(nn==mm)
            rhoxx(nn,mm,2)=1;
            rhoxx(nn,mm,3)=1;
            rhoxx(nn,mm,4)=1;
            rhoxx(nn,mm,5)=1;
            rhoxx(nn,mm,6)=1;
            rhoxx(nn,mm,7)=1;
            rhoxx(nn,mm,8)=1;
            
            
            
        elseif(nn<mm)
            rhoxx(nn,mm,8)=sign(CorrHxx(nn,mm,7)).^0.*(abs(CorrHxx(nn,mm,7))).^(1/7);
            rhoxx(nn,mm,7)=sign(CorrHxx(nn,mm,6)).^0.*(abs(CorrHxx(nn,mm,6))).^(1/6);
            rhoxx(nn,mm,6)=sign(CorrHxx(nn,mm,5)).^0.*(abs( (XCoeffH(6,nn).*CorrHxx(nn,mm,5)-21*XCoeffZ(8,nn).*rhoxx(nn,mm,8).^(5))./(XCoeffH(6,nn)-21*XCoeffZ(8,nn)))).^(1/5);
            rhoxx(nn,mm,5)=sign(CorrHxx(nn,mm,4)).^0.*(abs( (XCoeffH(5,nn).*CorrHxx(nn,mm,4)-15*XCoeffZ(7,nn).*rhoxx(nn,mm,7).^(4))./(XCoeffH(5,nn)-15*XCoeffZ(7,nn)))).^(1/4);
            rhoxx(nn,mm,4)=sign(CorrHxx(nn,mm,3)).^0.*(abs( (XCoeffH(4,nn).*CorrHxx(nn,mm,3)-10*XCoeffZ(6,nn).*rhoxx(nn,mm,6).^(3)-105*XCoeffZ(8,nn).*rhoxx(nn,mm,8).^(3))...
                ./(XCoeffH(4,nn)-10*XCoeffZ(6,nn)-105*XCoeffZ(8,nn)))).^(1/3);
            rhoxx(nn,mm,3)=sign(CorrHxx(nn,mm,2)).^0.*(abs( (XCoeffH(3,nn).*CorrHxx(nn,mm,2)-6*XCoeffZ(5,nn).*rhoxx(nn,mm,5).^(2)-45*XCoeffZ(7,nn).*rhoxx(nn,mm,7).^(2)) ...
                ./(XCoeffH(3,nn)-6*XCoeffZ(5,nn)-45*XCoeffZ(7,nn)))).^(1/2);
            rhoxx(nn,mm,2)=sign(CorrHxx(nn,mm,1)).^0.*(( (XCoeffH(2,nn).*CorrHxx(nn,mm,1)-3*XCoeffZ(4,nn).*rhoxx(nn,mm,4).^(1)-15*XCoeffZ(6,nn).*rhoxx(nn,mm,6).^(1) ...
                -105*XCoeffZ(8,nn).*rhoxx(nn,mm,8).^(1))...
                ./(XCoeffH(2,nn)-3*XCoeffZ(4,nn)-15*XCoeffZ(6,nn)-105*XCoeffZ(8,nn)))).^(1/1);
        

           % rhoxx1(nn,mm,6)=sign(CorrHxx(nn,mm,5)).*(abs( (CorrHxx(nn,mm,5)-21.*rhoxx(nn,mm,8).^(5)))).^(1/5);
           % rhoxx1(nn,mm,5)=(abs( CorrHxx(nn,mm,4)-15.*rhoxx(nn,mm,7).^(4))).^(1/4);
           %rhoxx1(mm,nn,6)=rhoxx1(nn,mm,6);
           % rhoxx1(mm,nn,5)=rhoxx1(nn,mm,5);
            
           % rhoxx(mm,nn,8)=rhoxx(nn,mm,8);
           % rhoxx(mm,nn,7)=rhoxx(nn,mm,7);

           rhoxx(mm,nn,8)=rhoxx(nn,mm,8);
            rhoxx(mm,nn,7)=rhoxx(nn,mm,7);
           rhoxx(mm,nn,6)=rhoxx(nn,mm,6);

            rhoxx(mm,nn,5)=rhoxx(nn,mm,5);
            rhoxx(mm,nn,4)=rhoxx(nn,mm,4);
            rhoxx(mm,nn,3)=rhoxx(nn,mm,3);
            rhoxx(mm,nn,2)=rhoxx(nn,mm,2);
        end
    end
end

rhoxy

rhoxx

str=input("Look at rhoxy and rhoxx at end");
end

.
.
.
.

function [cnOptim] = OptimizeCorrFromHermitesNewton(Zy,Zx,cnn,HCorrXX,NDims,paths,Order)

%Order;Max Order of Hermites Calibration is seven.
%Zx is paths X K array for joint values of Zx arranged for "paths" number of observations
%Zy is paths X 1 array for Zy associated with each obsewrvation.
%cnn is the original input of correlation coefficients used to derive seed for Newton methjod.
%cnn is in format of the old program and its order and indexing is changed
%in this NEwton program.

ccn(1:NDims+1,1:Order)=cnn(1:NDims+1,2:Order+1); %Shift the indices back by one index.
%so that index matches with order of associated hermite polynomial.

cn(1:NDims,1:Order)=ccn(2:NDims+1,1:Order);%The first row in each hermite column earlier was 
%associated zeta but zeta is shifted from first row to (K+1)th row.

cn(NDims+1,1:Order)=ccn(1,1:Order);

% cn(1:NDims+1,1:Order)=1/(NDims+1);
% 
% for nn=1:NDims
%     cn(nn,1:Order)=cn(nn,1:Order)+sum(HCorrXX(nn,:,1));
% end


%The function below calculates the first hermites till "Order" from data of
%Zy. WHy is (paths X K) array. kth hermite associated with pth observation appears in
%WHy(p,k)
[WHy] = EvalHermiteArrayH1(Zy,Order);

%The program below calculates moments of standard normal. Odd moments are
%zero. 
[EZ] = EZZ2();

%The loop below repeats itself for optimization of correlation coefficients 
%associated with highest power in each hermite. It repeatrs over all hermites one by one. 
for nH=1:Order
    cnIn(1:NDims+1,1:nH)=0.0;
    if(nH==1)
        cnIn(1:NDims+1,1)=cn(1:NDims+1,1);
    else
        %cnOptim is already optimized value of lower powers of Gaussians
        %calculated in earlier loops of optimization over smaller hermites.
        %This is already known from previous iterations of this loop.
        cnIn(1:NDims+1,1:nH-1)=cnOptim(1:NDims+1,1:nH-1);
        %The value of cnIn below is initial guess for correlation
        %coefficients associated with optimization of current hermite
        cnIn(1:NDims+1,nH)=cn(1:NDims+1,nH);
    end
    
    %Below nH-th hermite observation values are assigned to one dimensional array for optimization target in Newton. 
    WHyn(1:paths,1)=WHy(1:paths,nH);
    
    %Total iterations of Newton root search.
    Iterations=500;
    %Below is the function that actually does the Multi-dimensional Newton.
    [cnbest,InitialObj(nH),FinalObj(nH)] = OptimizenthHermiteGaussianWithNewton(cnIn,Zx,WHyn,HCorrXX,nH,EZ,NDims,paths, Iterations);

    %Below is optimum value of correlation coefficients associated with
    %nH-th hermite after optimization. This will be used by the program as
    %an input to optimization of higher hermites
    cnOptim(1:NDims+1,nH)=cnbest(1:NDims+1,1);
    
    nH
end


InitialObj %Obj function at start of optimization of correlation coefficients for nH-th hermite
FinalObj%Obj function after optimization of correlation coefficients for nH-th hermite
str=input("Look at initial and Final Objs");

end

.
.
.

function [bnbest,InitialObj,FinalObj] = OptimizenthHermiteGaussianWithNewton(cn,Zx,WHyn,HCorrXX,nH,EZ,NDims,paths,Iterations)

%Why0 in array of nth Hermite Values

KK=NDims;

bn(1:NDims+1,1)=cn(1:NDims+1,nH);
%Below, Calculate the model contribution of lower Gaussian powers that have
%already been calculated through associated lower hermites towards current
%hermite polynomial.
[dHermite] = EvaluatePrevGaussPowersForHCorrs(nH,cn,Zx,NDims,paths,EZ);
%Remove the model contribution of lower Gaussian powers from data value of
%current hermite.
WHy0=WHyn+dHermite;

%Below calcualte variance
var=0.0;
for kk1=1:KK
    for kk2=1:KK    
        var=var+bn(kk1,1).*bn(kk2,1).*sign(HCorrXX(kk1,kk2,nH)).*(abs(HCorrXX(kk1,kk2,nH)));%.^(1/nH);%Changed.
    end
end

var=var+bn(KK+1,1).^2;

bn

var
%Below, make sure that variance is one.
bn(:,1)=bn(:,1)/sqrt(var);

%str=input("have reached 1st Stone")        


%Below calculate the combined value of K correlated variables as described
%on Wilmott post and Calculate G.
GG1(1:paths,1)=0.0;
for kk=1:KK
    GG1(1:paths,1)=GG1(1:paths,1)+bn(kk,1).*Zx(1:paths,kk);
end

Zeta0=bn(KK+1,1);
%Evaluate the nHth Gaussian power below.
[GaussPH,GaussPh] = EvaluateHermitePolyForGaussCorrs(GG1,Zeta0,EZ,nH,paths);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Starting ObjBest and calculations here
%cnbest here
bnbest(1:NDims+1,1)=bn(1:NDims+1);

ObjSqrdBestNewton=sum((GaussPH(:,1)-WHy0(:,1)).^2)/paths;

InitialObj=ObjSqrdBestNewton;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%Below, enter the Newton loop.

StallCount=0;
iter=0;
Mul=.1;
while ((iter < Iterations) &&(StallCount<=5))
    iter=iter+1;
    dGG1(1:paths,1:NDims)=0.0;
    for kk=1:KK
        dGG1(1:paths,kk)=Zx(1:paths,kk);
    end
    
    bn=bnbest;
    
    GG1(1:paths,1)=0.0;
    for kk=1:KK
        GG1(1:paths,1)=GG1(1:paths,1)+bn(kk,1).*Zx(1:paths,kk);
    end

    Zeta0=bn(KK+1,1);
    
    bbn=0;
    bbn(1:NDims,1)=bn(1:NDims,1);

    
    [GaussPH,GaussPh,DGaussPH] = EvaluateHermitePolyAndDerivsForGaussCorrs(GG1,dGG1,Zeta0,EZ,nH,NDims,paths);


    %size(DGaussPH')
    %size(DGaussPH)
    %size(GaussPH)
    %size(WHy0(:,1))
    
    %str=input("Look at sizes");
    
    %I am using a multiple of .1 in Newton equation below. this might have
    %to be further decreased for some cases. Below is the Newton
    %over-constrained iteration.
    bbn=bbn+Mul*((DGaussPH'*DGaussPH)\(DGaussPH'*(WHy0(:,1)-GaussPH)));
    
    %DGaussPH'*DGaussPH
    
    %str=input("Look at matrix product");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Midfunction ObjBest and calculations here
%cnbest here

var=0.0;
for kk1=1:KK
    for kk2=1:KK    
        var=var+bbn(kk1,1).*bbn(kk2,1).*sign(HCorrXX(kk1,kk2,nH)).*(abs(HCorrXX(kk1,kk2,nH)));%.^(1/nH);%changed
    end
end

bbn(KK+1,1)=sqrt(1-var);

var=var+bbn(KK+1,1).^2;
bn
bbn

var

bbn=bbn/sqrt(var);



GG1(1:paths,1)=0.0;
for kk=1:KK
    GG1(1:paths,1)=GG1(1:paths,1)+bbn(kk,1).*Zx(1:paths,kk);
end


Zeta0=bbn(KK+1,1);

[GaussPH,GaussPh] = EvaluateHermitePolyForGaussCorrs(GG1,Zeta0,EZ,nH,paths);

ObjSqrdNewton=sum((GaussPH(:,1)-WHy0(:,1)).^2)/paths;
ObjSqrdBestNewton
if(ObjSqrdNewton<ObjSqrdBestNewton)
    if((ObjSqrdBestNewton-ObjSqrdNewton)/ObjSqrdBestNewton<.00001)
        StallCount=StallCount+1;
        Mul=Mul*.75
        
    end
        
    ObjSqrdBestNewton=ObjSqrdNewton;
    bnbest=bbn;
    
else
    
    StallCount=StallCount+1;
    %%%%change-step or Newton parameters
end
    
ObjSqrdBestNewton
ObjSqrdNewton
iter

%str=input("Look at the solution");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5    

end

FinalObj=ObjSqrdBestNewton;

end

.
.
.

function [CoeffR,XCoeffH,ZI] = RegressYOnX(Y,Xin,NData,Dim,Order)


[ZI] = StandardNormalIdealDesity02(NData,1);
for mm=1:Dim
    %[CoeffX(mm,1:Order+1),Zx(1:NData,mm)] = CalculateHermiteSeriesFromData(Xin(1:NData,mm),Order,ZI);
    [XCoeffH(1:Order+1,mm),Zx(1:NData,mm)] = CalculateHermiteSeriesFromData02(Xin(1:NData,mm),Order,ZI);

    [Zx(1:NData,mm)] = CalculateZgivenXAndHSeriesCoeffs(Xin(1:NData,mm),XCoeffH(1:Order+1,mm),Order);
end

%Zx

%str=input("Match Zx with Zx in python")

WW(1:NData,1:Order+1)=EvalHermiteArrayH0(Zx(1:NData,1),Order);


%WW

%str=input("Show me match of first EvalHermite Array0 with python")



for mm=2:Dim
    WW(1:NData,(mm-1)*Order+2:mm*Order+1)=EvalHermiteArrayH1(Zx(1:NData,mm),Order);
    WW0(1:NData,1:Order)=EvalHermiteArrayH1(Zx(1:NData,mm),Order);
    %WW0
    %mm
    
 
    %str=input("Compare with EvalHermiteArray1")
    
    
end




Ymu(1:NData,1)=Y(1:NData);

WWW=(WW'*WW);

for mm=1:Dim*Order+1
    WWW(mm,mm)=WWW(mm,mm)+.000001;
end

CoeffR=(WWW)\(WW'*Ymu);









% 
% for nn=1:NData
%     YProjected(nn,1)=0.0;
%    for mm=1:Dim*Order+1
%    YProjected(nn,1)=YProjected(nn,1)+WW(nn,mm).*CoeffR(mm,1);    
%     
%    end
% end
% 
%  STot=sum((Ymu(1:NData,1)-sum(Ymu(1:NData,1))/NData).^2);
%  SRes=sum((Ymu(1:NData,1)-YProjected(1:NData,1)).^2);
% % 
%  RSquared=1-SRes/STot;
% % 
%  RSquared
% % 
% YRight=0.0;
% YWrong=0.0;
% YPRExp=0.0;
% YNRExp=0.0;
% YPWExp=0.0;
% YNWExp=0.0;
% for nn=1:NData
%     
%    for mm=1:Dim*Order+1
%    YProjected(nn,1)=YProjected(nn,1)+WW(nn,mm).*CoeffR(mm,1);    
%    end
%    if(Ymu(nn,1)>=0)
%        if(YProjected(nn,1)>=0)
%            YRight=YRight+1;
%            YPRExp=YPRExp+(Ymu(nn,1));
%        else
%            YWrong=YWrong+1;
%            YPWExp=YPWExp+(Ymu(nn,1));
%        end
%    end
%    if(Ymu(nn,1)<0)
%        if(YProjected(nn,1)<0)
%            YRight=YRight+1;
%            YNRExp=YNRExp+(Ymu(nn,1));
%        else
%            YWrong=YWrong+1;
%            YNWExp=YNWExp+Ymu(nn,1);
%        end
%    end
% end
% 
% YRight
% YWrong
% 
% YPRExp
% YNRExp
% 
% YPWExp
% YNWExp
% 
% abs(YPRExp)+abs(YNRExp)-abs(YPWExp)-abs(YNWExp)
% 
% str=input("Look at RSquared--0002")
% 
% 
% 
%  mdl=fitlm(WW,Ymu)    
%     
% str=input("Matlab's first regression")    
% 
% 




end

.
.
.

Amin · Today, 3:57 pm

Friends, I am posting a few graphs for some SDE scenarios. These are conditional density of added SDE at T2 given the starting values of component SDEs at time T1. This post will be followed by another post in which I will give brief notes. And in third post I will attach matlab code used to make these graphs.
.
.
.
Here are parameters of the SDE

%SDE1
thetaV(1)=1.0;
kappaV(1)=0.000;
gamma(1)=.75;
sigma0(1)=.5;

%SDE2
thetaV(2)=1.0;
kappaV(2)=0.000;
gamma(2)=.65;
sigma0(2)=.5;

%SDE3
thetaV(3)=1.0;
kappaV(3)=0.000;
gamma(3)=.95;
sigma0(3)=.45;

and here are the four conditional starting points for the SDEs in each of the following graphs.

Set1
V10(1)=1.2;
V10(2)=.65;
V10(3)=1.2;

%Set2
V10(1)=.75;
V10(2)=1.2;
V10(3)=.7;

%Set3
V10(1)=1.25;
V10(2)=1.35;
V10(3)=1.25;

%Set4
V10(1)=.85;
V10(2)=.81;%1.2;
V10(3)=.77;

I have made three sets of graphs for three different correlation scenarios.

and correlation matrix for first scenario has zero correlations
Corr1(1,1)=1;
Corr1(1,2)=0;
Corr1(1,3)=0;
Corr1(2,2)=1;
Corr1(2,1)=0;
Corr1(2,3)=0;
Corr1(3,1)=0;
Corr1(3,2)=0;
Corr1(3,3)=1;

Here are four sets of graphs(four sets of starting conditional values) with zero correlation scenario.
.
.

Correlation matrix for second scenario:

Corr1(1,1)=1;
Corr1(1,2)=.7/3;
Corr1(1,3)=.8/3;
Corr1(2,2)=1;
Corr1(2,1)=.7/3;
Corr1(2,3)=.5/3;
Corr1(3,1)=.8/3;
Corr1(3,2)=.5/3;
Corr1(3,3)=1;

.
.
.

Third correlation scenario:

Corr1(1,1)=1;
Corr1(1,2)=-.7/3;
Corr1(1,3)=.8/3;
Corr1(2,2)=1;
Corr1(2,1)=-.7/3;
Corr1(2,3)=-.5/3;
Corr1(3,1)=.8/3;
Corr1(3,2)=-.5/3;
Corr1(3,3)=1;

Amin · Today, 4:14 pm

Friends, here are some notes.

But first I would like to say that everyday, we are improving on previous work but it will still take a few weeks to make things totally perfect.

I noticed that conditional mean from regression is more precise (as compared to earlier algorithm that we have used) for low correlations scenarios, and negative correlation scenarios. For high correlations scenarios, the conditional mean from earlier algorithm is more precise than the one from regression. So basically these two types of SDEs, at least as a start, has to be dealt a bit differently. Large number of things are common but we still need some differences in algorithm to deal with high correlation SDEs and the rest of the SDEs.

I have in above graphs dealt with small to medium positive correlations, and mixed positive and negative correlations and zero correlations. I have not dealth with high positive correlations and high negative correlations with this algorithm.

The algorithm would be clear to friends from the code but I noticed that for most of the SDEs in this category (small to medium positiv-negative correlations), the conditional mean from earlier method was higher while true mean was lower. Our conditional mean from regression is very close to true mean. I noticed that error in coefficients of hermites was also proportional to mean from old method and true mean. (slightly greater than ratio of first power of both means and slightly lesser than ratio of second power of both condiitional means).

Using the above logic, I ad-hocly improved the coefficients by following line of code: (Here YProjected is conditional mean from regression while bh00 is conditional mean from old algorithm).

YCoeffz(2:Order+1)=YCoeffz(2:Order+1).*YProjected/bh00;

I noticed that using squred ratio was somewhat better but I remained conservative and used the above formula in all of the above graphs. The above equation does not hold for high positive correlations scenario since there YProjected is worse than bh00 and bh00 is very close to true conditional mean.

You have to use improved versions of the old algorithm for very high positive correlations scenarios.

There is still a lot of work that needs to be done and I am sure we would do that in the following few weeks to get all theory on solid grounds. We are making progress.

In next post, I will attach the code used to make the graphs in previous post.

Amin · Today, 4:26 pm

Please note that following program has some repetitions where some of the things are calculated again in different parts but a later version will make everything more compact. Please note that this is still a very intermediate preliminary version and many new versions will follow in coming weeks.

.
.

function []= ConditionalDensityHermitesNvariateNewton()

%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
%The present program works only when drift is zero (for pure noises)
%meaning mu1=0 and mu2=0;

Order=7;
%NDim=4;%Three assets and one SV.
%NMomentsY=Order+1;
%NMomentsX=Order+1;

%%%%%%%%%%%%%%%%%%%5
dt=.03125/2/8;
Tt=64*2*8;
T=Tt*dt;
seed0=12130649;
%seed0=94210876;
%seed0=14210876;
%seed0=53612297
seed0=97459782
rng(seed0, 'twister')
paths=50000;
TT1=100*8/2; %Transition distribution start
TT2=256*8/2; %Transition distribution end

%%%%%%%%%%%%%%%%%%%5

NDims=3;
% 
V0(1)=1.000;%.32;
V00(1)=V0(1);

%SDE1
 thetaV(1)=1.0;%.045;%1;%.04;
 kappaV(1)=0.000;%1.5;%1.5;
 gamma(1)=.75;%.75
 sigma0(1)=.5;%.35


 
%SDE2%
%thetaV(1)=.90;%.045;%1;%.04;
%kappaV(1)=1.000;%1.5;%1.5;
%gamma(1)=.75;%.75
%sigma0(1)=.85;%.35

%SDE3%
%thetaV(1)=1.350;%.045;%1;%.04;
%kappaV(1)=2.000;%1.5;%1.5;
%gamma(1)=.75;%.75
%sigma0(1)=.95;%.35


mu1(1)=kappaV(1)*thetaV(1);
mu2(1)=-kappaV(1);
beta1(1)=0;
beta2(1)=1;



V0(2)=1.000;%.32;
V00(2)=V0(2);

%SDE1
thetaV(2)=1.0;%.045;%1;%.04;
 kappaV(2)=0.000;%1.5;%1.5;
 gamma(2)=.65;%.75
 sigma0(2)=.5;%.35

%SDE2 
% thetaV(2)=1.0;%.045;%1;%.04;
% kappaV(2)=2.000;%1.5;%1.5;
% gamma(2)=.65;%.75
% sigma0(2)=.95;%.35

%SDE3
%thetaV(2)=.850;%.045;%1;%.04;
%kappaV(2)=.35000;%1.5;%1.5;
%gamma(2)=.65;%.75
%sigma0(2)=.55;%.35


mu1(2)=kappaV(2)*thetaV(2);
mu2(2)=-kappaV(2);
beta1(2)=0;
beta2(2)=1;



V0(3)=1.000;%.33;
V00(3)=V0(3);

%SDE1
thetaV(3)=1.0;%.045;%1;%.04;
 kappaV(3)=0.000;%1.5;%1.5;
 gamma(3)=.95;%.75
 sigma0(3)=.45;%.35

%SDE2 
% thetaV(3)=.750;%.045;%1;%.04;
% kappaV(3)=0.5000;%1.5;%1.5;
% gamma(3)=.95;%.75
% sigma0(3)=.45;%.35


%SDE3 
%thetaV(3)=1.150;%.045;%1;%.04;
%kappaV(3)=1.5000;%1.5;%1.5;
%gamma(3)=.95;%.75
%sigma0(3)=.75;%.35

mu1(3)=kappaV(3)*thetaV(3);
mu2(3)=-kappaV(3);
beta1(3)=0;
beta2(3)=1;



% Corr1(1,1)=1;
% Corr1(1,2)=.7;
% Corr1(1,3)=.8;
% Corr1(2,2)=1;
% Corr1(2,1)=.7;
% Corr1(2,3)=.5;
% Corr1(3,1)=.8;
% Corr1(3,2)=.5;
% Corr1(3,3)=1;


Corr1(1,1)=1;
Corr1(1,2)=.7/3;
Corr1(1,3)=.8/3;
Corr1(2,2)=1;
Corr1(2,1)=.7/3;
Corr1(2,3)=.5/3;
Corr1(3,1)=.8/3;
Corr1(3,2)=.5/3;
Corr1(3,3)=1;



Corr1(1,1)=1;
Corr1(1,2)=-.7/3;
Corr1(1,3)=.8/3;
Corr1(2,2)=1;
Corr1(2,1)=-.7/3;
Corr1(2,3)=-.5/3;
Corr1(3,1)=.8/3;
Corr1(3,2)=-.5/3;
Corr1(3,3)=1;



% Corr1(1,1)=1;
% Corr1(1,2)=-.27;
% Corr1(1,3)=-.28;
% Corr1(2,2)=1;
% Corr1(2,1)=-.27;
% Corr1(2,3)=-.15;
% Corr1(3,1)=-.28;
% Corr1(3,2)=-.15;
% Corr1(3,3)=1;
%87%Corr1(1:3,1:3)=0.0;

R=chol(Corr1);
% 
%  R(1:3,1:3)=0;
%  R(1,1)=1;
%  R(2,2)=1;
%  R(3,3)=1;


V0(4)=1.000;%.44;
V00(4)=V0(4);
thetaV(4)=1.0;%.045;%1;%.04;
kappaV(4)=1.000;%1.5;%1.5;
mu1(4)=kappaV(4)*thetaV(4);
mu2(4)=-kappaV(4);
beta1(4)=0;
beta2(4)=1;
gamma(4)=.95;%.75
sigma0(4)=.55;%.45


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
for nn=1:NDims
V(1:paths,nn)=V0(nn);
end
%V(1:paths,1:NDims)=V0(1,1:NDims);

Random1(1:paths/2,1:NDims)=0;

for tt=1:Tt
    
    Random2(1:paths/2,1:NDims)=randn(size(Random1));
    Random2(paths/2+1:paths,1:NDims)=-Random2(1:paths/2,1:NDims);

    %Random2=R'*Random2;
    random1=Random2*R(:,1);
    random2=Random2*R(:,2);
    random3=Random2*R(:,3);
    
    Random2(:,1)=random1;
    Random2(:,2)=random2;
    Random2(:,3)=random3;
    
    
    for nn=1:NDims
    
    V(1:paths,nn)=V(1:paths,nn)+ ...
        (mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn))*dt + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*Random2(1:paths,nn)*sqrt(dt) + ...
        (mu1(nn).*beta1(nn)*V(1:paths,nn).^(beta1(nn)-1) + mu2(nn).*beta2(nn).*V(1:paths,nn).^(beta2(nn)-1)).* ...
        ((mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn))*dt^2/2 + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*Random2(1:paths,nn)*(1-1/sqrt(3))*dt^1.5) + ...
        .5*(mu1(nn).*beta1(nn).*(beta1(nn)-1).*V(1:paths,nn).^(beta1(nn)-2) + mu2(nn).*beta2(nn).*(beta2(nn)-1).*V(1:paths,nn).^(beta2(nn)-2)).* ...    
        sigma0(nn)^2.*V(1:paths,nn).^(2*gamma(nn)).*dt^2/2 + ...
        sigma0(nn).*gamma(nn).*V(1:paths,nn).^(gamma(nn)-1) .* ...
        ((mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn)).*Random2(1:paths,nn).*1/sqrt(3)*dt^1.5 + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*(Random2(1:paths,nn).^2-1)*dt/2) + ...
        .5*sigma0(nn).*gamma(nn).*(gamma(nn)-1).*V(1:paths,nn).^(gamma(nn)-2) .* ...
        sigma0(nn).^2.*V(1:paths,nn).^(2*gamma(nn)) .*Random2(1:paths,nn).*1/sqrt(3)*dt^1.5;

    V(V(:,nn)<=0)=.0001;

    
    if(tt==TT1)
        Xin(1:paths,nn)=real(V(1:paths,nn));
    end
        
    if(tt==TT2)
        
        Xin2(1:paths,nn)=real(V(1:paths,nn));
    end

    end    
    
end


Yin(1:paths)=0.0;

for pp=1:paths
for mm=1:NDims
    Yin(pp)=Yin(pp)+Xin2(pp,mm);%-Xin(pp,mm);
%-Xin(pp,mm);
end
end


str=input("I have reached stone 1");


% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
% 



%save("C:\Users\Lenovo\Documents\MATLAB\MATLAB1\YXfile.mat","Yin","Xin")


%str=input("Data has been saved");




% load("C:\Users\Lenovo\Documents\MATLAB\MATLAB1\YXfile.mat","Yin","Xin")
% 
% str=input("Data has been loaded");
% 
% 
% 
% 
 Order=7;
 clf;
[BetaCorrH,YCoeffH,XCoeffH,HCorrXX,Zy,Zx] = CalculateHermiteCorrelations0NA04(Xin,Yin,paths,NDims);
YCoeffZ(1:Order+1)= ConvertHermiteSeriesToZSeries01(YCoeffH(1:Order+1),Order);

%HCorrXX(:,:)=0.0;

BetaCorrH

HCorrXX
%CorrXX1st
str=input("I have reached Stone2");


%[cn,dn,bnv,bbnv,ch] = Calculate2DZSeriesNDimCorrected(YCoeffH,BetaCorrH,HCorrXX,NDims,Order);
%cn=bnv;
[cn,ch]=CalculateGaussianCorrelationFromHermiteCorrelationsNDimNew(YCoeffH,XCoeffH,BetaCorrH,HCorrXX,NDims,Order);
% 
% 
% 
% save("C:\Users\Lenovo\Documents\MATLAB\MATLAB1\Zyxfile.mat","Zy","Zx","cn","HCorrXX")
% 
% 
% str=input("Data has been saved");

%load("C:\Users\Lenovo\Documents\MATLAB\MATLAB1\Zyxfile.mat","Zy","Zx","cn","HCorrXX")

%str=input("Data has been loaded");

%size(Zy)
%size(Zx)
%size(cn)
HCorrXX0=HCorrXX;
[HCorrXY,HCorrXX,YCoeffHNorm,XCoeffHNorm]=CalculateHermiteCorrelationCoefficientsForGaussianCorrs(Xin,Yin,paths,NDims);

HCorrXX
HCorrXX0

str=input("Look at HCorrXX and HCorrXX0");

%[rhoxy,rhoxx] = CalculateGaussianCorrelationFromHermiteCorrelations001(YCoeffH,XCoeffH,HCorrXY,HCorrXX,NDims,Order)
%[rhoxy,rhoxx] = CalculateGaussianCorrelationFromHermiteCorrelations001(YCoeffH,XCoeffH,BetaCorrH,HCorrXX,NDims,Order)
[rhoxy,rhoxx] = CalculateGaussianCorrelationFromHermiteCorrelations001(YCoeffHNorm,XCoeffH,HCorrXY,HCorrXX0,NDims,Order);

%rhoxy
%rhoxx

%str=input("Look at these correlations")

%OptimizeCorrFromHermitesIterative(Zy,Zx,cn,HCorrXX,NDims,paths,Order);

%OptimizeCorrFromHermitesIterative02(Zy,Zx,cn,HCorrXX,NDims,paths,Order);

%OptimizeCorrFromHermitesNewton02(Zy,Zx,cn,HCorrXX,NDims,paths,Order);

% HCorrXX0=HCorrXX;
% for nn1=1:NDims
%     for nn2=1:NDims
%         VariableHIn(2:Order+1)=HCorrXX(nn1,nn2,1:Order);
%         VariableZOut(1:Order+1)= ConvertHermiteSeriesToZSeries01(VariableHIn(1:Order+1),Order);
%         HCorrXXZ(nn1,nn2,1:Order)=VariableZOut(2:Order+1);
%     end
% end



%[cnOptim] = OptimizeCorrFromHermitesNewton(Zy,Zx,cn,HCorrXX,NDims,paths,Order);
[cnOptim] = OptimizeCorrFromHermitesNewton(Zy,Zx,rhoxy,rhoxx(:,:,2:8),NDims,paths,Order);

cn=cnOptim;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5555

[CoeffR,XCoeffH,ZI] = RegressYOnX(Yin,Xin,paths,NDims,Order);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Set1
V10(1)=1.2;%.65;%
V10(2)=.65;
V10(3)=1.2;
%V10(4)=1.85;



% 
% % % % %Set2
 V10(1)=.75;
 V10(2)=1.2;%1.2;
 V10(3)=.7;
% % % %V10(4)=.65
% % % % 
% % % % %Set3
  V10(1)=1.25;
  V10(2)=1.35;%1.2;
  V10(3)=1.25;
% % %   V10(4)=1.4
% % % % 
% % % %   
  V10(1)=.85;
  V10(2)=.81;%1.2;
  V10(3)=.77;
% % %  V10(4)=.7

 

%V10(1)=.65;
%V10(2)=.51;%1.2;
%V10(3)=.73;

%V10(1)=1.65;
%V10(2)=1.51;%1.2;
%V10(3111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111111)=1.63;

for nn=1:NDims
V(1:paths,nn)=V10(nn);
end
%V(1:paths,1:NDims)=V0(1,1:NDims);

Random1(1:paths/2,1:NDims)=0;

for tt=1:TT2-TT1
    
    Random2(1:paths/2,1:NDims)=randn(size(Random1));
    Random2(paths/2+1:paths,1:NDims)=-Random2(1:paths/2,1:NDims);

    random1=Random2*R(:,1);
    random2=Random2*R(:,2);
    random3=Random2*R(:,3);
    
    Random2(:,1)=random1;
    Random2(:,2)=random2;
    Random2(:,3)=random3;
    
    
    
    for nn=1:NDims
    
    V(1:paths,nn)=V(1:paths,nn)+ ...
        (mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn))*dt + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*Random2(1:paths,nn)*sqrt(dt) + ...
        (mu1(nn).*beta1(nn)*V(1:paths,nn).^(beta1(nn)-1) + mu2(nn).*beta2(nn).*V(1:paths,nn).^(beta2(nn)-1)).* ...
        ((mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn))*dt^2/2 + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*Random2(1:paths,nn)*(1-1/sqrt(3))*dt^1.5) + ...
        .5*(mu1(nn).*beta1(nn).*(beta1(nn)-1).*V(1:paths,nn).^(beta1(nn)-2) + mu2(nn).*beta2(nn).*(beta2(nn)-1).*V(1:paths,nn).^(beta2(nn)-2)).* ...    
        sigma0(nn)^2.*V(1:paths,nn).^(2*gamma(nn)).*dt^2/2 + ...
        sigma0(nn).*gamma(nn).*V(1:paths,nn).^(gamma(nn)-1) .* ...
        ((mu1(nn).*V(1:paths,nn).^beta1(nn) + mu2(nn).*V(1:paths,nn).^beta2(nn)).*Random2(1:paths,nn).*1/sqrt(3)*dt^1.5 + ...
        sigma0(nn).*V(1:paths,nn).^gamma(nn) .*(Random2(1:paths,nn).^2-1)*dt/2) + ...
        .5*sigma0(nn).*gamma(nn).*(gamma(nn)-1).*V(1:paths,nn).^(gamma(nn)-2) .* ...
        sigma0(nn).^2.*V(1:paths,nn).^(2*gamma(nn)) .*Random2(1:paths,nn).*1/sqrt(3)*dt^1.5;

    V(V(:,nn)<0)=.0001;

    
    %if(tt==TT1)
    %    XXin(1:paths,nn)=V(1:paths,nn);
    %end
        
    if(tt==TT2-TT1)
        
        XXin2(1:paths,nn)=V(1:paths,nn);
    end

    end    
    
end


Yin2(1:paths)=0.0;

for pp=1:paths
for mm=1:NDims
    Yin2(pp)=Yin2(pp)+XXin2(pp,mm);%-V10(mm);
end
end



randn()

str=input("I have reached stone 4");



NoOfBins=200;
MaxCutOff=20;

[XDensity,IndexOut,IndexMax] = MakeDensityFromSimulation_Infiniti_NEW(Yin2,paths,NoOfBins,MaxCutOff );

for qq=1:8
    
    DataMoments(qq)=sum(Yin2(:).^qq)/paths;
    
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

for nn=1:NDims
    [Zx0(nn)] = CalculateZgivenXAndHSeriesCoeffs(V10(nn),XCoeffH(1:Order+1,nn),Order);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55

WW(1,1:Order+1)=EvalHermiteArrayH0(Zx0(1),Order);
for mm=2:NDims
    WW(1,(mm-1)*Order+2:mm*Order+1)=EvalHermiteArrayH1(Zx0(mm),Order);
end    




     YProjected=0.0;
    for mm=1:NDims*Order+1
    YProjected=YProjected+WW(1,mm).*CoeffR(mm,1);    
     
    end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
Zx(1:Order)=0; 

for hh=1:Order
    for nn=1:NDims
        %Zx(hh) = Zx(hh) + bnv(nn+1,hh+1).^1.*Zx0(nn);
        Zx(hh) = Zx(hh) + cn(nn,hh).^1.*Zx0(nn);
    end
end

bh0=YCoeffH(1);

[W] = EvalHermiteArrayH0(Zx0,Order);

for hh=1:Order
    for nn=1:NDims
        bh0 = bh0 + ch(nn+1,hh+1).^1.*W(nn,hh+1);
    end
end
bh00=bh0;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


a0=YCoeffZ(1);
a1=YCoeffZ(2);
a2=YCoeffZ(3);
a3=YCoeffZ(4);
a4=YCoeffZ(5);
a5=YCoeffZ(6);
a6=YCoeffZ(7);
a7=YCoeffZ(8);


d1=cn(4,1);
d2=cn(4,2);
d3=cn(4,3);
d4=cn(4,4);
d5=cn(4,5);
d6=cn(4,6);
d7=cn(4,7);


%b0=(a0+a1.* Zx(1) +a2.* Zx(2).^2+a3.* Zx(3).^3+a4* Zx(4).^4+a5 *Zx(5).^5+a6* Zx(6).^6+a7* Zx(7).^7);
b1=(a1* d1+2* a2* Zx(2)* d2 +3* a3* Zx(3).^2* d3+4* a4* Zx(4).^3* d4+5* a5* Zx(5).^4* d5+6* a6* Zx(6).^5* d6+7* a7* Zx(7).^6* d7);
b2=(a2* d2.^2+3* a3* Zx(3)* d3^2 +6* a4* Zx(4).^2* d4^2+10* a5* Zx(5).^3* d5^2+15* a6* Zx(6).^4* d6^2+21* a7* Zx(7).^5* d7^2);
b3=(a3* d3^3+4* a4 *Zx(4)* d4^3 +10* a5* Zx(5).^2* d5^3+20* a6* Zx(6).^3* d6^3+35* a7* Zx(7).^4* d7^3);
b4=(a4* d4^4+5* a5* Zx(5)* d5^4 +15* a6 *Zx(6).^2* d6^4+35* a7* Zx(7).^3 *d7^4);
b5=(a5* d5^5+6* a6* Zx(6)* d6^5 +21* a7 *Zx(7).^2* d7^5);
b6=(a6* d6^6+7* a7* Zx(7)* d7^6 );
b7=(a7* d7^7);


%bh0=b0+b2+3*b4+15*b6;
bh1=b1+3*b3+15*b5+105*b7;
bh2=b2+6*b4+45*b6;
bh3=b3+10*b5+105*b7;
bh4=b4+15*b6;
bh5=b5+21*b7;
bh6=b6;
bh7=b7;


bh(1)=YProjected;
bh(2)=bh1;
bh(3)=bh2;
bh(4)=bh3;
bh(5)=bh4;
bh(6)=bh5;
bh(7)=bh6;
bh(8)=bh7;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5

YCoeffz(1:Order+1)= ConvertHermiteSeriesToZSeries01(bh(1:Order+1),Order);


[YMomentsModel]=CalculateMomentsOfZSeries(YCoeffz(1),YCoeffz(2:Order+1),Order,Order+1);
 
YMomentsModel
 
DataMoments
 
str=input("Look at comparison of moments");

[cMu1,YCentralMoments]=ConvertRawMomentsToCentralMoments(YMomentsModel,Order+1);
[cMu2,DataCentralMoments]=ConvertRawMomentsToCentralMoments(DataMoments,Order+1);


YCentralMoments

DataCentralMoments

bh00
YProjected

YCoeffz(2:Order+1)=YCoeffz(2:Order+1).*YProjected/bh00;
str=input("Look at central moments");

% Order0=7;
% Order1=13;
% NData=paths;
% [CondMomentYgivX,BetaCorrHxy] = CalculateConditionalMomentYgivXNDims(Xin,Yin,NData,Zx0,NDims,Order0,Order1)
% 
% [VarConditional]=FindSecondConditionalMoment(YCoeffH,BetaCorrH,Zx0)
% 
% 
% BetaCorrHxy
% BetaCorrH
% 
% str=input("Look at conditionalVariance");

clf;

PlotZSeriesDensity(YCoeffz(1),YCoeffz(2:Order+1),'b')


hold on

plot(IndexOut(1:IndexMax),XDensity(1:IndexMax),'r');

%PlotHermiteSeriesDensityAndRvGraph(YCoeffh(1),YCoeffh1(2:yOrder+1),'g')
title("Graph of Density of sum of Three Correlated SDE Variables conditional on starting values of the three SDE random variables.")
%title(sprintf('V00 = %.3f; kappa=%.3f, theta=%.3f,gamma=%.3f,sigma=%.3f;V10=%.3f at TT1=%.3f;TT2=%.3f;Paths=%d ',V00, kappaV, thetaV,gamma,sigma0, V10,(TT1*dt),(TT2*dt),paths));%,sprintf('theta= %f', theta), sprintf('kappa = %f', kappa),sprintf('sigma = %f', sigma0),sprintf('T = %f', T));
legend({'Analytical Density','Numerical Density1'}, ...
    'Location','northeast')

str=input("This is the first version without altering the moments of conditional density")

hold off 
 





end

.
.
.
.

function [BetaCorrH,YCoeffH,XCoeffH,HCorrXX,Zy,Zx] = CalculateHermiteCorrelations0NA04(Xin,Yin,NData,NDims)

%This program finds required correlations between NDims explanatory
%variables given in Xin that has NDims dimensionality. Yin is variable
%being explained and is one-diemsnional.


Order=7;
NMoments=8;

[ZI] = MomentMatchedStandardNormalDiscretizedDensity(NData);


for nn=1:NDims
[XCoeffH(1:Order+1,nn),Zx(1:NData,nn)] = CalculateHermiteSeriesFromData02(Xin(1:NData,nn),Order,ZI);

[Zx(1:NData,nn)] = CalculateZgivenXAndHSeriesCoeffs(Xin(1:NData,nn),XCoeffH(1:Order+1,nn),Order);

end


[YCoeffH,Zy] = CalculateHermiteSeriesFromData02(Yin,Order,ZI); 
[Zy] = CalculateZgivenXAndHSeriesCoeffs(Yin,YCoeffH,Order);


for nn1=1:NDims
    for nn2=1:NDims
        if(nn1==nn2)
            HCorrXX(nn1,nn2,1:Order)=1.0;
        else
            HCorrXX(nn1,nn2,1:Order) = CalculateCorrelationBivariateHermiteCH0(Zx(1:NData,nn1),Zx(1:NData,nn2),NData,Order);
        end
    end
end
for nn1=1:NDims
    HCorrXX(nn1,nn1,1:Order)=1.0;
end


for hh=1:Order
    for nn1=1:NDims
        HCorrXX
        ss2=0;
        for nn2=1:NDims
            if(nn2~=nn1)
                ss2=ss2+1;
                ss3=0;
                for nn3=1:NDims
                    if(nn3~=nn1)
                        ss3=ss3+1;
                            HCorrXXmn(ss2,ss3)=HCorrXX(nn2,nn3,hh);
                    end
                    HCorrXmXn(ss2,1)=HCorrXX(nn1,nn2,hh);
                end
            end
        end
        BetaCorrXX(1:NDims-1,nn1,hh)=HCorrXXmn(:,:)\HCorrXmXn(1:NDims-1,1);
            
    end
        
    BetaCorrXX
        
    nn1
       
        str=input("Look at important variables");
    
end
BetaCorrXX            
str=input("look at the BetacorrXX");       
for nn1=1:NDims
    for nn2=1:NDims-1
        if(nn2<nn1)
    BetaCorrXX0(nn2,nn1,1:Order)=BetaCorrXX(nn2,nn1,1:Order);
        elseif(nn2==nn1)
               % BetaCorrXX0(nn2,nn1,1:Order)=1;
            
            BetaCorrXX0(nn2+1:NDims,nn1,1:Order)=BetaCorrXX(nn2:NDims-1,nn1,1:Order);
        end
    end
    BetaCorrXX0(nn1,nn1,1:Order)=1;
end
BetaCorrXX
BetaCorrXX0
% 
% % % % str=input("Look at BetaCorrXX pairs");
% % % % for hh=1:Order
% % % %     for nn1=1:NDims
% % % %         for nn2=1:NDims-1
% % % %             BetaCorrXXn(nn1,nn2,hh)=sign(BetaCorrXXn(nn1,nn2,hh)).*(abs(BetaCorrXXn(nn1,nn2,hh))).^(1/hh);
% % % %         end
% % % %     end
% % % % end



for nn1=1:NDims
    HCorrXY(nn1,1:Order) = CalculateCorrelationBivariateHermiteCH0(Zx(1:NData,nn1),Zy(1:NData),NData,Order);
end

for hh=1:Order
    [BCorrXY(:,hh),Is(:,hh)]=sort(HCorrXY(:,hh),'descend');
end

for hh=1:Order
    CorrXX(1:NDims,1:NDims)=HCorrXX(Is(1:NDims,hh),Is(1:NDims,hh),hh);
    CorrXY(1:NDims,1)=HCorrXY(Is(1:NDims,hh),hh);
    BetaCorr0(1:NDims)=CorrXX\CorrXY;
    BetaCorrH(Is(1:NDims,hh),hh)=BetaCorr0(1:NDims);
end



HCorrXX=BetaCorrXX0;
HCorrXX
str=input("Look at HCOrrXX");

% for hh=1:Order
%     CorrXX(1:NDims,1:NDims)=HCorrXX(1:NDims,1:NDims,hh);
%     CorrXY(1:NDims,1)=HCorrXY(1:NDims,hh);
%     BetaCorrH(1:NDims,hh)=CorrXX\CorrXY;
% end




%CorrXX1st(1:NDims,1:NDims)=HCorrXX(1:NDims,1:NDims,1)


BetaCorrH(:,1)'*CorrXX

CorrXY

HCorrXX

HCorrXY

BetaCorrH

str=input("Look at BetaCorrH");
end

.
.
.
.

function [cn,ch] = CalculateGaussianCorrelationFromHermiteCorrelationsNDimNew(YCoeffH,XCoeffH,BetaCorrHxy,CorrHxx,NDims,Order)





for hh=1:Order
    for nn=1:NDims
        ch(nn+1,hh+1)=BetaCorrHxy(nn,hh).*YCoeffH(hh+1);%make sure indices on all hermite-correlations match.check.
        %%%% nn+1 since 1st index would go to y-tilde.
    end
end


ch
YCoeffH

str=input("Look at above---00-1-1-1-1-");

[YCoeffZ(1:Order+1)] = ConvertHermiteSeriesToZSeries01(YCoeffH(1:Order+1),Order);

for nn=1:NDims
    [XCoeffZ(1:Order+1,nn)] = ConvertHermiteSeriesToZSeries01(XCoeffH(1:Order+1,nn),Order);
end




YCoeffH

YCoeffZ

str=input("Look at the Hermite and Z-series coeffs");

% YCoeffZA=YCoeffZ;
% YCoeffZ=YCoeffZA;
% if(abs(YCoeffZ(4))>abs(YCoeffZ(2)))
%    YCoeffZ(4)=sign(YCoeffZ(4)).*abs(YCoeffZ(2));
% end
% if(abs(YCoeffZ(6))>abs(YCoeffZ(4)))
%    YCoeffZ(6)=sign(YCoeffZ(6)).*abs(YCoeffZ(4));
% end
% if(abs(YCoeffZ(8))>abs(YCoeffZ(6)))
%    YCoeffZ(6)=sign(YCoeffZ(6)).*abs(YCoeffZ(8));
% end
% 
% 
% if(abs(YCoeffZ(5))>abs(YCoeffZ(3)))
%    YCoeffZ(5)=sign(YCoeffZ(5)).*abs(YCoeffZ(3));
% end
% if(abs(YCoeffZ(7))>abs(YCoeffZ(5)))
%    YCoeffZ(5)=sign(YCoeffZ(5)).*abs(YCoeffZ(7));
% end
% 




for nn=1:NDims
[c(nn+1,1:Order+1)] = ConvertHermiteSeriesToZSeries01(ch(nn+1,1:Order+1),Order);
end


for nn=1:NDims
    rhoxy(nn+1,8)=sign(BetaCorrHxy(nn,7)).^1.*(abs(BetaCorrHxy(nn,7))).^(1/7);
    rhoxy(nn+1,7)=sign(BetaCorrHxy(nn,6)).^1.*(abs(BetaCorrHxy(nn,6))).^(1/6);
    rhoxy(nn+1,6)=sign(BetaCorrHxy(nn,5)).^1.*(abs(   (YCoeffH(6).*BetaCorrHxy(nn,5)-21*YCoeffZ(8).*rhoxy(nn+1,8).^(5))./(YCoeffH(6)-21*YCoeffZ(8)))).^(1/5);
    rhoxy(nn+1,5)=sign(BetaCorrHxy(nn,4)).^1.*(abs(  (YCoeffH(5).*BetaCorrHxy(nn,4)-15*YCoeffZ(7).*rhoxy(nn+1,7).^(4))./(YCoeffH(5)-15*YCoeffZ(7)))).^(1/4);
    rhoxy(nn+1,4)=sign(BetaCorrHxy(nn,3)).^1.*(abs(  (YCoeffH(4).*BetaCorrHxy(nn,3)-10*YCoeffZ(6).*rhoxy(nn+1,6).^(3)-105*YCoeffZ(8).*rhoxy(nn+1,8).^(3))...
        ./(YCoeffH(4)-10*YCoeffZ(6)-105*YCoeffZ(8)))).^(1/3);
    rhoxy(nn+1,3)=sign(BetaCorrHxy(nn,2)).^1.*(abs( (YCoeffH(3).*BetaCorrHxy(nn,2)-6*YCoeffZ(5).*rhoxy(nn+1,5).^(2)-45*YCoeffZ(7).*rhoxy(nn+1,7).^(2)) ...
        ./(YCoeffH(3)-6*YCoeffZ(5)-45*YCoeffZ(7)))).^(1/2);
    rhoxy(nn+1,2)=sign(BetaCorrHxy(nn,1)).^1.*(abs(  (YCoeffH(2).*BetaCorrHxy(nn,1)-3*YCoeffZ(4).*rhoxy(nn+1,4).^(1)-15*YCoeffZ(6).*rhoxy(nn+1,6).^(1) ...
        -105*YCoeffZ(8).*rhoxy(nn+1,8).^(1))...
        ./(YCoeffH(2)-3*YCoeffZ(4)-15*YCoeffZ(6)-105*YCoeffZ(8)))).^(1/1);
    
%    rhoxy(nn,2+1)=sign(4/9*BetaCorrHxy(nn,2)+1/36).*abs((4/9*BetaCorrHxy(nn,2)+1/36))^(.5);
% rhoxy(nn,3+1)=sign((36/225*BetaCorrHxy(nn,3)-3/225*BetaCorrHxy(nn,1))).*(abs((36/225*BetaCorrHxy(nn,3)-3/225*BetaCorrHxy(nn,1))))^(1/3);
% rhoxy(nn,4+1)=sign((144*4/(105)^2*BetaCorrHxy(nn,4)-6/105^2*(4*BetaCorrHxy(nn,2)+.25)-18/105^2)).*(abs(((144*4/(105)^2*BetaCorrHxy(nn,4)-6/105^2*(4*BetaCorrHxy(nn,2)+.25)-18/105^2))))^(1/4);
% rhoxy(nn,1+1)=BetaCorrHxy(nn,1);
end


% 
for nn=1:NDims
    for mm=1:NDims
        if(nn==mm)
            rhoxx(nn,mm,2)=1;
            rhoxx(nn,mm,3)=1;
            rhoxx(nn,mm,4)=1;
            rhoxx(nn,mm,5)=1;
            rhoxx(nn,mm,6)=1;
            rhoxx(nn,mm,7)=1;
            rhoxx(nn,mm,8)=1;
            
            
            
        elseif(nn<mm)
            rhoxx(nn,mm,8)=sign(CorrHxx(nn,mm,7)).*(abs(CorrHxx(nn,mm,7))).^(1/7);
            rhoxx(nn,mm,7)=sign(CorrHxx(nn,mm,6)).*(abs(CorrHxx(nn,mm,6))).^(1/6);
            rhoxx(nn,mm,6)=sign(CorrHxx(nn,mm,5)).*(abs( (XCoeffH(6,nn).*CorrHxx(nn,mm,5)-21*XCoeffZ(8,nn).*rhoxx(nn,mm,8).^(5))./(XCoeffH(6,nn)-21*XCoeffZ(8,nn)))).^(1/5);
            rhoxx(nn,mm,5)=sign(CorrHxx(nn,mm,4)).*(abs( (XCoeffH(5,nn).*CorrHxx(nn,mm,4)-15*XCoeffZ(7,nn).*rhoxx(nn,mm,7).^(4))./(XCoeffH(5,nn)-15*XCoeffZ(7,nn)))).^(1/4);
            rhoxx(nn,mm,4)=sign(CorrHxx(nn,mm,3)).*(abs( (XCoeffH(4,nn).*CorrHxx(nn,mm,3)-10*XCoeffZ(6,nn).*rhoxx(nn,mm,6).^(3)-105*XCoeffZ(8,nn).*rhoxx(nn,mm,8).^(3))...
                ./(XCoeffH(4,nn)-10*XCoeffZ(6,nn)-105*XCoeffZ(8,nn)))).^(1/3);
            rhoxx(nn,mm,3)=sign(CorrHxx(nn,mm,2)).*(abs( (XCoeffH(3,nn).*CorrHxx(nn,mm,2)-6*XCoeffZ(5,nn).*rhoxx(nn,mm,5).^(2)-45*XCoeffZ(7,nn).*rhoxx(nn,mm,7).^(2)) ...
                ./(XCoeffH(3,nn)-6*XCoeffZ(5,nn)-45*XCoeffZ(7,nn)))).^(1/2);
            rhoxx(nn,mm,2)=sign(CorrHxx(nn,mm,1)).*(abs( (XCoeffH(2,nn).*CorrHxx(nn,mm,1)-3*XCoeffZ(4,nn).*rhoxx(nn,mm,4).^(1)-15*XCoeffZ(6,nn).*rhoxx(nn,mm,6).^(1) ...
                -105*XCoeffZ(8,nn).*rhoxx(nn,mm,8).^(1))...
                ./(XCoeffH(2,nn)-3*XCoeffZ(4,nn)-15*XCoeffZ(6,nn)-105*XCoeffZ(8,nn)))).^(1/1);
        

           % rhoxx1(nn,mm,6)=sign(CorrHxx(nn,mm,5)).*(abs( (CorrHxx(nn,mm,5)-21.*rhoxx(nn,mm,8).^(5)))).^(1/5);
           % rhoxx1(nn,mm,5)=(abs( CorrHxx(nn,mm,4)-15.*rhoxx(nn,mm,7).^(4))).^(1/4);
           %rhoxx1(mm,nn,6)=rhoxx1(nn,mm,6);
           % rhoxx1(mm,nn,5)=rhoxx1(nn,mm,5);
            
           % rhoxx(mm,nn,8)=rhoxx(nn,mm,8);
           % rhoxx(mm,nn,7)=rhoxx(nn,mm,7);
rhoxx(mm,nn,6)=rhoxx(nn,mm,6);
rhoxx(mm,nn,5)=rhoxx(nn,mm,5);
rhoxx(mm,nn,4)=rhoxx(nn,mm,4);
rhoxx(mm,nn,3)=rhoxx(nn,mm,3);
rhoxx(mm,nn,2)=rhoxx(nn,mm,2);
        end
    end
end
% % 
% 
% rhoxx(:,:,6)
% rhoxx1(:,:,6)
% rhoxx(:,:,5)
% rhoxx1(:,:,5)
% 
% str=input("Look at rhoxx and rhoxx1");

% 
% % 
% % 
% % 
% % %rhoxx(1:NDims,1:NDims,2:8)=CorrHxx(1:NDims,1:NDims,1:7);
%  for hh=1:7
% % %rhoxx(1:NDims,1:NDims,hh+1)=rhoxx(1:NDims,1:NDims,2);
%  rhoxx(1:NDims,1:NDims,hh+1)=(sign(CorrHxx(1:NDims,1:NDims,hh))).^(1).*(abs(CorrHxx(1:NDims,1:NDims,hh))).^(1/hh);
%  for nn=1:NDims
%      for mm=1:NDims
%          if(nn~=mm)
%  rhoxx1(nn,mm,2+1)=-sign(4/9*CorrHxx(nn,mm,2)+1/36).^0.*abs((4/9*CorrHxx(nn,mm,2)+1/36))^(.5);
%  rhoxx1(nn,mm,3+1)=-sign((36/225*CorrHxx(nn,mm,3)-3/225*CorrHxx(nn,mm,1))).^0.*(abs((36/225*CorrHxx(nn,mm,3)-3/225*CorrHxx(nn,mm,1))))^(1/3);
%  rhoxx1(nn,mm,4+1)=-sign((144*4/(105)^2*CorrHxx(nn,mm,4)-6/105^2*(4*CorrHxx(nn,mm,2)+.25)-18/105^2)).^0.*(abs(((144*4/(105)^2*CorrHxx(nn,mm,4)-6/105^2*(4*CorrHxx(nn,mm,2)+.25)-18/105^2))))^(1/4);
%  %rhoxx1(nn,mm,2+1)=sign(4/9*CorrHxx(nn,mm,2)+1/36).^(2).*abs((4/9*CorrHxx(nn,mm,2)+1/36))^(.5);
%  %rhoxx1(nn,mm,3+1)=sign((36/225*CorrHxx(nn,mm,3)-3/225*CorrHxx(nn,mm,1))).^(3).*(abs((36/225*CorrHxx(nn,mm,3)-3/225*CorrHxx(nn,mm,1))))^(1/3);
%  %rhoxx1(nn,mm,4+1)=sign((144*4/(105)^2*CorrHxx(nn,mm,4)-6/105^2*(4*CorrHxx(nn,mm,2)+.25)-18/105^2)).^(4).*(abs(((144*4/(105)^2*CorrHxx(nn,mm,4)-6/105^2*(4*CorrHxx(nn,mm,2)+.25)-18/105^2))))^(1/4);
%  
%          else
%              rhoxx1(nn,mm,2+1)=1.0;
%               rhoxx1(nn,mm,3+1)=1.0;
%                rhoxx1(nn,mm,4+1)=1.0;
%          end
%      end
%  end
%  end
%  
%  rhoxx(1:NDims,1:NDims,3)=rhoxx1(1:NDims,1:NDims,3);
%  rhoxx(1:NDims,1:NDims,4:5)=rhoxx1(1:NDims,1:NDims,4:5);
%  
 for hh=1:7
 %rhoxx(1:NDims,1:NDims,hh+1)=sign(CorrHxx(1:NDims,1:NDims,hh)).*abs(sign(CorrHxx(1:NDims,1:NDims,hh))).^(1/hh);
 end
 %rhoxx
 %rhoxx1
% 
% 
 str=input("Look at rhoxx");

% 
% CorrHxx0(1:NDims,1:NDims,2:Order+1)=CorrHxx(1:NDims,1:NDims,1:Order);
% 
% for nn=1:NDims
%     for mm=1:NDims
% 
%         if(nn==mm)
%             rhoxx(nn,mm,2)=1;
%             rhoxx(nn,mm,3)=1;
%             rhoxx(nn,mm,4)=1;
%             rhoxx(nn,mm,5)=1;
%             rhoxx(nn,mm,6)=1;
%             rhoxx(nn,mm,7)=1;
%             rhoxx(nn,mm,8)=1;
%             
%         else
% [rrhoxx(nn+1,mm+1,1:Order+1)] = ConvertHermiteSeriesToZSeries01(CorrHxx0(nn,mm,1:Order+1),Order);
% end
% end
% end
% 
% for hh=1:7
% %rhoxx(1:NDims,1:NDims,hh+1)=rhoxx(1:NDims,1:NDims,2);
% rhoxx(1:NDims,1:NDims,hh+1)=sign(rrhoxx(1:NDims,1:NDims,hh+1)).*(abs(rrhoxx(1:NDims,1:NDims,hh+1))).^(1/hh);
% end







%rhoxx(:,:)=0.0;

%(YCoeffH(6).*BetaCorrHxy(nn,5)-21*YCoeffZ(8).*rhoxy(nn+1,8).^(5/7))
%(abs((abs(YCoeffH(6)).*BetaCorrHxy(nn,5)-21*abs(YCoeffZ(8)).*rhoxy(nn+1,8).^(5/7))./abs(YCoeffH(6)-21*(YCoeffZ(8))))).^1/5
%(YCoeffH(6)-21*YCoeffZ(8))
%(YCoeffZ(6)-21*YCoeffZ(8))

%(abs(abs(YCoeffH(6)).*BetaCorrHxy(nn,5)-21*abs(YCoeffZ(8)).*rhoxy(nn+1,8).^(5/7))).^(1/5)



%YCoeffH
%YCoeffZ
%rhoxy

%str=input("Look at rhoxy");



cn=rhoxy;





%for nn=1:NDims
%    rho_z(nn+1,7)=c(nn+1,8)./YCoeffZ(8);
%    cn(nn+1,8)=(abs(rho_z(nn+1,7))).^(1/7);
%    cn(nn+1,8)=cn(nn+1,8).*sign(rho_z(nn+1,7));
%end

cnSqrSum=0;
for nn=1:NDims
    %ccn0(nn+1,8)=cn(nn+1,8).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,8).*cn(nn+1,8);
    for mm=1:NDims
        if(nn~=mm)
            %ccn0(nn+1,8)=ccn0(nn+1,8)+2*cn(nn+1,8).*cn(mm+1,8).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,8).*cn(mm+1,8).*rhoxx(nn,mm,8);
        end
    end
    %ccn(nn+1,8)=sqrt(1-ccn0(nn+1,8));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,8)=cn(2:NDims+1,8)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end

cn(1,8)=sqrt(1-cnSqrSum);
%ccn(1,8)=sqrt(1-cn(1,8).^2);
%Cn(1,8)=cn(1,8).*YCoeffZ(8);

cn
%ccn
str=input("Compare the Eighth Element From NDim---2");


% ccn00(2,8)=cn(2,8).^2+2*cn(2,8).*cn(3,8).*CorrXX1st(1,2)+2*cn(2,8).*cn(4,8).*CorrXX1st(1,3)+ ...
%     cn(3,8).^2+2*cn(3,8).*cn(4,8).*CorrXX1st(2,3)+cn(4,8).^2
% 
% ccn00(2,8)=sqrt(1-(cn(2,8).^2+2*cn(2,8).*cn(3,8).*CorrXX1st(1,2)+2*cn(2,8).*cn(4,8).*CorrXX1st(1,3)))
% 
% cnSqrSum
% 
% ccn0(2,8)
% 
% ccn(1,8)
% 
% ccn(2,8)
% 
% str=input("Look at the variables Inside");
% 







cnSqrSum=0;
%for nn=1:NDims
%    rho_z(nn+1,6)=c(nn+1,7)./YCoeffZ(7);
%    cn(nn+1,7)=(abs(rho_z(nn+1,6))).^(1/6);
%    cn(nn+1,7)=cn(nn+1,7).*sign(rho_z(nn+1,6));
%    %ccn(nn+1,7)=sqrt(1-cn(nn+1,7).^2);
%    %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
%end

cnSqrSum=0;
for nn=1:NDims
    %ccn0(nn+1,7)=cn(nn+1,7).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,7).*cn(nn+1,7);
    for mm=1:NDims
        if(nn~=mm)
         %   ccn0(nn+1,7)=ccn0(nn+1,7)+2*cn(nn+1,7).*cn(mm+1,7).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,7).*cn(mm+1,7).*rhoxx(nn,mm,7);
        end
    end
    %ccn(nn+1,7)=sqrt(1-ccn0(nn+1,7));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,7)=cn(2:NDims+1,7)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end

cn(1,7)=sqrt(1-cnSqrSum);
%ccn(1,7)=sqrt(1-cn(1,7).^2);
%Cn(1,7)=cn(1,7).*YCoeffZ(7);


cn
%ccn
str=input("Compare the Seventh Element From NDim---2");


%cnSqrSum=0;
%for nn=1:NDims
%    CC(3,6)=YCoeffZ(8).*(21*cn(nn+1,8).^5.*ccn(nn+1,8).^2);
%    rho_z(nn+1,5)=(c(nn+1,6)-CC(3,6))./YCoeffZ(6);
%    if(abs(rho_z(nn+1,5))>.995)
%        rho_z(nn+1,5)=sign(rho_z(nn+1,5)).*.995;
%    end
%    cn(nn+1,6)=(abs(rho_z(nn+1,5))).^(1/5);
%    cn(nn+1,6)=cn(nn+1,6).*sign(rho_z(nn+1,5));
%%%    nn
%%%    c(nn+1,6)
%%%    CC(3,6)
%%%    YCoeffZ(8)
%%%    YCoeffZ(6)
%%%    rho_z(nn+1,5)
%%%    cn(nn+1,6)
%%%    
%%%    str=input("Look inside the correlations module loop");
    
    
    
    
    
%    ccn(nn+1,6)=sqrt(1-cn(nn+1,6).^2);
%    %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
%end

cnSqrSum=0;
for nn=1:NDims
%    ccn0(nn+1,6)=cn(nn+1,6).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,6).*cn(nn+1,6);
    for mm=1:NDims
        if(nn~=mm)
          %  ccn0(nn+1,6)=ccn0(nn+1,6)+2*cn(nn+1,6).*cn(mm+1,6).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,6).*cn(mm+1,6).*rhoxx(nn,mm,6);
        end
    end
    %ccn(nn+1,6)=sqrt(1-ccn0(nn+1,6));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,6)=cn(2:NDims+1,6)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end



cn(1,6)=sqrt(1-cnSqrSum);
%ccn(1,6)=sqrt(1-cn(1,6).^2);
%Cn(1,6)=cn(1,6).*YCoeffZ(6);

cn
%ccn
str=input("Compare the Sixth Element From NDim---2");



% cnSqrSum=0;
% for nn=1:NDims
%     CC(3,5)=YCoeffZ(7).*(15*cn(nn+1,7).^4.*ccn(nn+1,7).^2);
%     rho_z(nn+1,4)=(c(nn+1,5)-CC(3,5))./YCoeffZ(5);
%     cn(nn+1,5)=(abs(rho_z(nn+1,4))).^(1/4);
%     cn(nn+1,5)=cn(nn+1,5).*sign(rho_z(nn+1,4));
%     %ccn(nn+1,5)=sqrt(1-cn(nn+1,5).^2);
%     %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
% end

cnSqrSum=0;
for nn=1:NDims
    %ccn0(nn+1,5)=cn(nn+1,5).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,5).*cn(nn+1,5);
    for mm=1:NDims
        if(nn~=mm)
           % ccn0(nn+1,5)=ccn0(nn+1,5)+2*cn(nn+1,5).*cn(mm+1,5).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,5).*cn(mm+1,5).*rhoxx(nn,mm,5);
        end
    end
    %ccn(nn+1,5)=sqrt(1-ccn0(nn+1,5));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,5)=cn(2:NDims+1,5)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end


cn(1,5)=sqrt(1-cnSqrSum);
%ccn(1,5)=sqrt(1-cn(1,5).^2);
%Cn(1,5)=cn(1,5).*YCoeffZ(5);


cn
%ccn
str=input("Compare the Fifth Element From NDim---2");




%CC(3,4)=YCoeffZ(6).*(10*cn(6).^3.*dn(6).^2*zeta^2);
%CC(5,4)=YCoeffZ(8).*(35*cn(8).^3.*dn(8).^4.*zeta^4);
% cnSqrSum=0;
% for nn=1:NDims
%     CC(5,4)=YCoeffZ(8).*(35*cn(nn+1,8).^3.*ccn(nn+1,8).^4);
%     CC(3,4)=YCoeffZ(6).*(10*cn(nn+1,6).^3.*ccn(nn+1,6).^2);
%     rho_z(nn+1,3)=(c(nn+1,4)-3*CC(5,4)-CC(3,4))./YCoeffZ(4);
%     cn(nn+1,4)=(abs(rho_z(nn+1,3))).^(1/3);
%     cn(nn+1,4)=cn(nn+1,4).*sign(rho_z(nn+1,3));
%     %ccn(nn+1,4)=sqrt(1-cn(nn+1,4).^2);
%     %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
% end

cnSqrSum=0;
for nn=1:NDims
%    ccn0(nn+1,4)=cn(nn+1,4).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,4).*cn(nn+1,4);
    for mm=1:NDims
        if(nn~=mm)
          %  ccn0(nn+1,4)=ccn0(nn+1,4)+2*cn(nn+1,4).*cn(mm+1,4).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,4).*cn(mm+1,4).*rhoxx(nn,mm,4);
        end
    end
    %ccn(nn+1,4)=sqrt(1-ccn0(nn+1,4));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,4)=cn(2:NDims+1,4)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end

cn(1,4)=sqrt(1-cnSqrSum);
%ccn(1,4)=sqrt(1-cn(1,4).^2);
%Cn(1,4)=cn(1,4).*YCoeffZ(4);

cn
%ccn
str=input("Compare the Fourth Element From NDim---2");



%CC(3,3)=YCoeffZ(5).*(6* cn(5).^2.* dn(5).^2);
%CC(5,3)=YCoeffZ(7).*(15*cn(7).^2.*dn(7).^4);
% cnSqrSum=0;
% for nn=1:NDims
%     CC(3,3)=YCoeffZ(5).*(6*cn(nn+1,5).^2.*ccn(nn+1,5).^2);
%     CC(5,3)=YCoeffZ(7).*(15*cn(nn+1,7).^2.*ccn(nn+1,7).^4);
%     rho_z(nn+1,2)=(c(nn+1,3)-3*CC(5,3)-CC(3,3))./YCoeffZ(3);
%     cn(nn+1,3)=(abs(rho_z(nn+1,2))).^(1/2);
%     cn(nn+1,3)=cn(nn+1,3).*sign(rho_z(nn+1,2));
%     %ccn(nn+1,3)=sqrt(1-cn(nn+1,3).^2);
%     %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
% end

cnSqrSum=0;
for nn=1:NDims
    %ccn0(nn+1,3)=cn(nn+1,3).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,3).*cn(nn+1,3);
    for mm=1:NDims
        if(nn~=mm)
           % ccn0(nn+1,3)=ccn0(nn+1,3)+2*cn(nn+1,3).*cn(mm+1,3).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,3).*cn(mm+1,3).*rhoxx(nn,mm,3);
        end
    end
    %ccn(nn+1,3)=sqrt(1-ccn0(nn+1,3));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,3)=cn(2:NDims+1,3)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end

cn(1,3)=sqrt(1-cnSqrSum);
%ccn(1,3)=sqrt(1-cn(1,3).^2);
%Cn(1,3)=cn(1,3).*YCoeffZ(3);

cn
%ccn
str=input("Compare the Third Element From NDim---2");



%CC(7,2)=YCoeffZ(8).*(7*cn(8).*dn(8).^6*zeta^6);
%CC(5,2)=YCoeffZ(6).*(5*cn(6).*dn(6).^4.*zeta^4);
%CC(3,2)=YCoeffZ(4).* (3.* cn(4).* dn(4).^2.* zeta^2);
% cnSqrSum=0;
% for nn=1:NDims
%     CC(4,2)=YCoeffZ(4).*(3*cn(nn+1,4).*ccn(nn+1,4).^2);
%     CC(6,2)=YCoeffZ(6).*(5*cn(nn+1,5).*ccn(nn+1,6).^4);
%     CC(8,2)=YCoeffZ(8).*(7*cn(nn+1,8).*ccn(nn+1,8).^6);
%     rho_z(nn+1,1)=(c(nn+1,2)-15*CC(8,2)-3*CC(6,2)-CC(4,2))./YCoeffZ(2);
%     cn(nn+1,2)=(abs(rho_z(nn+1,1)));
%     %ccn(nn+1,2)=sqrt(1-cn(nn+1,2).^2);
%     %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
% end

cnSqrSum=0;
for nn=1:NDims
%    ccn0(nn+1,2)=cn(nn+1,2).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,2).*cn(nn+1,2);
    for mm=1:NDims
        if(nn~=mm)
        %    ccn0(nn+1,2)=ccn0(nn+1,2)+2*cn(nn+1,2).*cn(mm+1,2).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,2).*cn(mm+1,2).*rhoxx(nn,mm,2);
        end
    end
%    ccn(nn+1,2)=sqrt(1-ccn0(nn+1,2));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,2)=cn(2:NDims+1,2)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end
cn(1,2)=sqrt(1-cnSqrSum);
%ccn(1,2)=sqrt(1-cn(1,2).^2);
%Cn(1,2)=cn(1,2).*YCoeffZ(2);

cn
%ccn
str=input("Compare the Second Element From NDim---2");


% 
% %CC(3,1)=YCoeffZ(2).*(ccn(nn+1,2).^2);
% %CC(5,1)=YCoeffZ(4).*(ccn(nn+1,4).^4);
% %CC(7,1)=YCoeffZ(6).*(ccn(nn+1,6).^6);
% CC(3,1)=YCoeffZ(2).*(cn(1,2).^2);
% CC(5,1)=YCoeffZ(4).*(cn(1,4).^4);
% CC(7,1)=YCoeffZ(6).*(cn(1,6).^6);
% 
% 
% 
% 
% 
% cn(1,1)=(YCoeffZ(1)-CC(3,1)-3*CC(5,1)-15*CC(7,1))./YCoeffZ(1);
% 
% 
% cn
% ccn
% str=input("Look at coefficients cn----PleASE lOOK");
% 
% 
 dn(1:Order+1)=cn(1,1:Order+1);
% bn(1:Order+1)=ccn(1,1:Order+1);
% 
% 
% 
% 

end

.
.
.

function [HCorrXY,HCorrXX,YCoeffHNorm,XCoeffHNorm] = CalculateHermiteCorrelationCoefficientsForGaussianCorrs(Xin,Yin,NData,NDims)

Order=7;
NMoments=8;

[ZI] = MomentMatchedStandardNormalDiscretizedDensity(NData);


for nn=1:NDims
[XCoeffH(1:Order+1,nn),Zx(1:NData,nn)] = CalculateHermiteSeriesFromData02(Xin(1:NData,nn),Order,ZI);

[Zx(1:NData,nn)] = CalculateZgivenXAndHSeriesCoeffs(Xin(1:NData,nn),XCoeffH(1:Order+1,nn),Order);

end


[YCoeffH,Zy] = CalculateHermiteSeriesFromData02(Yin,Order,ZI); 
[Zy] = CalculateZgivenXAndHSeriesCoeffs(Yin,YCoeffH,Order);


for nn1=1:NDims
    XXin(1:NData,nn1)=Xin(1:NData,nn1)./(sum(Xin(1:NData,nn1).^2)/NData-(sum(Xin(1:NData,nn1))/NData).^2);
    ss2=0;
    for nn2=1:NDims
        if(nn2~=nn1)
            ss2=ss2+1;
            XX(1:NData,ss2)=Xin(1:NData,nn2);%/(sum(Xin(1:NData,nn2).^2)/NData-(sum(Xin(1:NData,nn2))/NData).^2);
        end
    end
    [CoeffR(1,:),XCoeffH,ZI] = RegressYOnX(XXin(1:NData,nn1),XX,NData,NDims-1,Order);
    CoeffRX(:,nn1)=CoeffR';
end

%for mm=1:NDims*Order+1
%    YProjected=YProjected+WW(1,mm).*CoeffR(mm,1);        
%end

for nn=1:NDims
    ss=1;
    for mm=1:NDims-1
        if(mm<nn)
            for hh=1:Order
                ss=ss+1;
                HCorrXX(mm,nn,hh)=CoeffRX(ss,nn);
            end
        elseif(mm==nn)
            for hh=1:Order
                ss=ss+1;
                HCorrXX(mm,nn,hh)=1.0;
                HCorrXX(mm+1,nn,hh)=CoeffRX(ss,nn);
            end
        elseif(mm>nn)
            for hh=1:Order
                ss=ss+1;
                HCorrXX(mm+1,nn,hh)=CoeffRX(ss,nn);
            end
        end
    end
end
HCorrXX(NDims,NDims,1:Order)=1;

for nn=1:NDims
[XCoeffHNorm(1:Order+1,nn),Zx(1:NData,nn)] = CalculateHermiteSeriesFromData02(XXin(1:NData,nn),Order,ZI);
end


YYin(1:NData)=Yin(1:NData)./(sum(Yin(1:NData).^2)/NData-(sum(Yin(1:NData))/NData).^2);
[CoeffRY(1,:),XCoeffH,ZI] = RegressYOnX(YYin(1:NData),XXin,NData,NDims,Order);            


ss=1;
for nn=1:NDims
    for hh=1:Order
        ss=ss+1;
        HCorrXY(nn,hh)=CoeffRY(1,ss);
    end
end
        
CoeffRY
HCorrXY


[YCoeffHNorm(1:Order+1),Zy(1:NData)] = CalculateHermiteSeriesFromData02(YYin(1:NData),Order,ZI);


%CoeffRX

%HCorrXX
        
str=input("Look at comparison of values"); 


            
end

.
.
.
.

function [cn,rhoxx] = CalculateGaussianCorrelationFromHermiteCorrelations001(YCoeffH,XCoeffH,BetaCorrHxy,CorrHxx,NDims,Order)





for hh=1:Order
    for nn=1:NDims
        ch(nn+1,hh+1)=BetaCorrHxy(nn,hh).*YCoeffH(hh+1);%make sure indices on all hermite-correlations match.check.
        %%%% nn+1 since 1st index would go to y-tilde.
    end
end


ch
YCoeffH

str=input("Look at above---00-1-1-1-1-");

[YCoeffZ(1:Order+1)] = ConvertHermiteSeriesToZSeries01(YCoeffH(1:Order+1),Order);

for nn=1:NDims
    [XCoeffZ(1:Order+1,nn)] = ConvertHermiteSeriesToZSeries01(XCoeffH(1:Order+1,nn),Order);
end




YCoeffH

YCoeffZ

str=input("Look at the Hermite and Z-series coeffs-----here---here");

% YCoeffZA=YCoeffZ;
% YCoeffZ=YCoeffZA;
% if(abs(YCoeffZ(4))>abs(YCoeffZ(2)))
%    YCoeffZ(4)=sign(YCoeffZ(4)).*abs(YCoeffZ(2));
% end
% if(abs(YCoeffZ(6))>abs(YCoeffZ(4)))
%    YCoeffZ(6)=sign(YCoeffZ(6)).*abs(YCoeffZ(4));
% end
% if(abs(YCoeffZ(8))>abs(YCoeffZ(6)))
%    YCoeffZ(6)=sign(YCoeffZ(6)).*abs(YCoeffZ(8));
% end
% 
% 
% if(abs(YCoeffZ(5))>abs(YCoeffZ(3)))
%    YCoeffZ(5)=sign(YCoeffZ(5)).*abs(YCoeffZ(3));
% end
% if(abs(YCoeffZ(7))>abs(YCoeffZ(5)))
%    YCoeffZ(5)=sign(YCoeffZ(5)).*abs(YCoeffZ(7));
% end
% 




for nn=1:NDims
[c(nn+1,1:Order+1)] = ConvertHermiteSeriesToZSeries01(ch(nn+1,1:Order+1),Order);
end


for nn=1:NDims
    rhoxy(nn+1,8)=sign(BetaCorrHxy(nn,7)).^1.*(abs(BetaCorrHxy(nn,7))).^(1/7);
    rhoxy(nn+1,7)=sign(BetaCorrHxy(nn,6)).^1.*(abs(BetaCorrHxy(nn,6))).^(1/6);
    rhoxy(nn+1,6)=sign(BetaCorrHxy(nn,5)).^1.*(abs(   (YCoeffH(6).*BetaCorrHxy(nn,5)-21*YCoeffZ(8).*rhoxy(nn+1,8).^(5))./(YCoeffH(6)-21*YCoeffZ(8)))).^(1/5);
    rhoxy(nn+1,5)=sign(BetaCorrHxy(nn,4)).^1.*(abs(  (YCoeffH(5).*BetaCorrHxy(nn,4)-15*YCoeffZ(7).*rhoxy(nn+1,7).^(4))./(YCoeffH(5)-15*YCoeffZ(7)))).^(1/4);
    rhoxy(nn+1,4)=sign(BetaCorrHxy(nn,3)).^1.*(abs(  (YCoeffH(4).*BetaCorrHxy(nn,3)-10*YCoeffZ(6).*rhoxy(nn+1,6).^(3)-105*YCoeffZ(8).*rhoxy(nn+1,8).^(3))...
        ./(YCoeffH(4)-10*YCoeffZ(6)-105*YCoeffZ(8)))).^(1/3);
    rhoxy(nn+1,3)=sign(BetaCorrHxy(nn,2)).^1.*(abs( (YCoeffH(3).*BetaCorrHxy(nn,2)-6*YCoeffZ(5).*rhoxy(nn+1,5).^(2)-45*YCoeffZ(7).*rhoxy(nn+1,7).^(2)) ...
        ./(YCoeffH(3)-6*YCoeffZ(5)-45*YCoeffZ(7)))).^(1/2);
    rhoxy(nn+1,2)=sign(BetaCorrHxy(nn,1)).^0.*(((  (YCoeffH(2).*BetaCorrHxy(nn,1)-3*YCoeffZ(4).*rhoxy(nn+1,4).^(1)-15*YCoeffZ(6).*rhoxy(nn+1,6).^(1) ...
        -105*YCoeffZ(8).*rhoxy(nn+1,8).^(1))...
        ./(YCoeffH(2)-3*YCoeffZ(4)-15*YCoeffZ(6)-105*YCoeffZ(8)))).^(1/1));
    
%    rhoxy(nn,2+1)=sign(4/9*BetaCorrHxy(nn,2)+1/36).*abs((4/9*BetaCorrHxy(nn,2)+1/36))^(.5);
% rhoxy(nn,3+1)=sign((36/225*BetaCorrHxy(nn,3)-3/225*BetaCorrHxy(nn,1))).*(abs((36/225*BetaCorrHxy(nn,3)-3/225*BetaCorrHxy(nn,1))))^(1/3);
% rhoxy(nn,4+1)=sign((144*4/(105)^2*BetaCorrHxy(nn,4)-6/105^2*(4*BetaCorrHxy(nn,2)+.25)-18/105^2)).*(abs(((144*4/(105)^2*BetaCorrHxy(nn,4)-6/105^2*(4*BetaCorrHxy(nn,2)+.25)-18/105^2))))^(1/4);
% rhoxy(nn,1+1)=BetaCorrHxy(nn,1);
end


% 
for nn=1:NDims
    for mm=1:NDims
        if(nn==mm)
            rhoxx(nn,mm,2)=1;
            rhoxx(nn,mm,3)=1;
            rhoxx(nn,mm,4)=1;
            rhoxx(nn,mm,5)=1;
            rhoxx(nn,mm,6)=1;
            rhoxx(nn,mm,7)=1;
            rhoxx(nn,mm,8)=1;
            
            
            
        else    %if(nn<mm)
            rhoxx(nn,mm,8)=sign(CorrHxx(nn,mm,7)).^1.*(abs(CorrHxx(nn,mm,7))).^(1/7);
            rhoxx(nn,mm,7)=sign(CorrHxx(nn,mm,6)).^1.*(abs(CorrHxx(nn,mm,6))).^(1/6);
            rhoxx(nn,mm,6)=sign(CorrHxx(nn,mm,5)).^1.*(abs( (XCoeffH(6,nn).*CorrHxx(nn,mm,5)-21*XCoeffZ(8,nn).*rhoxx(nn,mm,8).^(5))./(XCoeffH(6,nn)-21*XCoeffZ(8,nn)))).^(1/5);
            rhoxx(nn,mm,5)=sign(CorrHxx(nn,mm,4)).^1.*(abs( (XCoeffH(5,nn).*CorrHxx(nn,mm,4)-15*XCoeffZ(7,nn).*rhoxx(nn,mm,7).^(4))./(XCoeffH(5,nn)-15*XCoeffZ(7,nn)))).^(1/4);
            rhoxx(nn,mm,4)=sign(CorrHxx(nn,mm,3)).^1.*(abs( (XCoeffH(4,nn).*CorrHxx(nn,mm,3)-10*XCoeffZ(6,nn).*rhoxx(nn,mm,6).^(3)-105*XCoeffZ(8,nn).*rhoxx(nn,mm,8).^(3))...
                ./(XCoeffH(4,nn)-10*XCoeffZ(6,nn)-105*XCoeffZ(8,nn)))).^(1/3);
            rhoxx(nn,mm,3)=sign(CorrHxx(nn,mm,2)).^1.*(abs( (XCoeffH(3,nn).*CorrHxx(nn,mm,2)-6*XCoeffZ(5,nn).*rhoxx(nn,mm,5).^(2)-45*XCoeffZ(7,nn).*rhoxx(nn,mm,7).^(2)) ...
                ./(XCoeffH(3,nn)-6*XCoeffZ(5,nn)-45*XCoeffZ(7,nn)))).^(1/2);
            rhoxx(nn,mm,2)=sign(CorrHxx(nn,mm,1)).^0.*((( (XCoeffH(2,nn).*CorrHxx(nn,mm,1)-3*XCoeffZ(4,nn).*rhoxx(nn,mm,4).^(1)-15*XCoeffZ(6,nn).*rhoxx(nn,mm,6).^(1) ...
                -105*XCoeffZ(8,nn).*rhoxx(nn,mm,8).^(1))...
                ./(XCoeffH(2,nn)-3*XCoeffZ(4,nn)-15*XCoeffZ(6,nn)-105*XCoeffZ(8,nn)))).^(1/1));
        

           % rhoxx1(nn,mm,6)=sign(CorrHxx(nn,mm,5)).*(abs( (CorrHxx(nn,mm,5)-21.*rhoxx(nn,mm,8).^(5)))).^(1/5);
           % rhoxx1(nn,mm,5)=(abs( CorrHxx(nn,mm,4)-15.*rhoxx(nn,mm,7).^(4))).^(1/4);
           %rhoxx1(mm,nn,6)=rhoxx1(nn,mm,6);
           % rhoxx1(mm,nn,5)=rhoxx1(nn,mm,5);
            
           % rhoxx(mm,nn,8)=rhoxx(nn,mm,8);
           % rhoxx(mm,nn,7)=rhoxx(nn,mm,7);

%            rhoxx(mm,nn,8)=rhoxx(nn,mm,8);
%             rhoxx(mm,nn,7)=rhoxx(nn,mm,7);
%            rhoxx(mm,nn,6)=rhoxx(nn,mm,6);
% 
%             rhoxx(mm,nn,5)=rhoxx(nn,mm,5);
%             rhoxx(mm,nn,4)=rhoxx(nn,mm,4);
%             rhoxx(mm,nn,3)=rhoxx(nn,mm,3);
%             rhoxx(mm,nn,2)=rhoxx(nn,mm,2);
        end
    end
end



for nn=1:NDims
    for mm=1:NDims
        rhoxx0(mm,nn,2:Order+1)=(rhoxx(mm,nn,2:Order+1)+rhoxx(nn,mm,2:Order+1))/2;
    end
end

rhoxy

rhoxx0

str=input("Look at rhoxy and rhoxx at end");

rhoxx=rhoxx0;

cn=rhoxy;





%for nn=1:NDims
%    rho_z(nn+1,7)=c(nn+1,8)./YCoeffZ(8);
%    cn(nn+1,8)=(abs(rho_z(nn+1,7))).^(1/7);
%    cn(nn+1,8)=cn(nn+1,8).*sign(rho_z(nn+1,7));
%end

cnSqrSum=0;
for nn=1:NDims
    %ccn0(nn+1,8)=cn(nn+1,8).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,8).*cn(nn+1,8);
    for mm=1:NDims
        if(nn~=mm)
            %ccn0(nn+1,8)=ccn0(nn+1,8)+2*cn(nn+1,8).*cn(mm+1,8).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,8).*cn(mm+1,8).*rhoxx(nn,mm,8);
        end
    end
    %ccn(nn+1,8)=sqrt(1-ccn0(nn+1,8));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,8)=cn(2:NDims+1,8)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end

cn(1,8)=sqrt(1-cnSqrSum);
%ccn(1,8)=sqrt(1-cn(1,8).^2);
%Cn(1,8)=cn(1,8).*YCoeffZ(8);

cn
%ccn
str=input("Compare the Eighth Element From NDim---2");


% ccn00(2,8)=cn(2,8).^2+2*cn(2,8).*cn(3,8).*CorrXX1st(1,2)+2*cn(2,8).*cn(4,8).*CorrXX1st(1,3)+ ...
%     cn(3,8).^2+2*cn(3,8).*cn(4,8).*CorrXX1st(2,3)+cn(4,8).^2
% 
% ccn00(2,8)=sqrt(1-(cn(2,8).^2+2*cn(2,8).*cn(3,8).*CorrXX1st(1,2)+2*cn(2,8).*cn(4,8).*CorrXX1st(1,3)))
% 
% cnSqrSum
% 
% ccn0(2,8)
% 
% ccn(1,8)
% 
% ccn(2,8)
% 
% str=input("Look at the variables Inside");
% 







cnSqrSum=0;
%for nn=1:NDims
%    rho_z(nn+1,6)=c(nn+1,7)./YCoeffZ(7);
%    cn(nn+1,7)=(abs(rho_z(nn+1,6))).^(1/6);
%    cn(nn+1,7)=cn(nn+1,7).*sign(rho_z(nn+1,6));
%    %ccn(nn+1,7)=sqrt(1-cn(nn+1,7).^2);
%    %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
%end

cnSqrSum=0;
for nn=1:NDims
    %ccn0(nn+1,7)=cn(nn+1,7).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,7).*cn(nn+1,7);
    for mm=1:NDims
        if(nn~=mm)
         %   ccn0(nn+1,7)=ccn0(nn+1,7)+2*cn(nn+1,7).*cn(mm+1,7).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,7).*cn(mm+1,7).*rhoxx(nn,mm,7);
        end
    end
    %ccn(nn+1,7)=sqrt(1-ccn0(nn+1,7));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,7)=cn(2:NDims+1,7)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end

cn(1,7)=sqrt(1-cnSqrSum);
%ccn(1,7)=sqrt(1-cn(1,7).^2);
%Cn(1,7)=cn(1,7).*YCoeffZ(7);


cn
%ccn
str=input("Compare the Seventh Element From NDim---2");


%cnSqrSum=0;
%for nn=1:NDims
%    CC(3,6)=YCoeffZ(8).*(21*cn(nn+1,8).^5.*ccn(nn+1,8).^2);
%    rho_z(nn+1,5)=(c(nn+1,6)-CC(3,6))./YCoeffZ(6);
%    if(abs(rho_z(nn+1,5))>.995)
%        rho_z(nn+1,5)=sign(rho_z(nn+1,5)).*.995;
%    end
%    cn(nn+1,6)=(abs(rho_z(nn+1,5))).^(1/5);
%    cn(nn+1,6)=cn(nn+1,6).*sign(rho_z(nn+1,5));
%%%    nn
%%%    c(nn+1,6)
%%%    CC(3,6)
%%%    YCoeffZ(8)
%%%    YCoeffZ(6)
%%%    rho_z(nn+1,5)
%%%    cn(nn+1,6)
%%%    
%%%    str=input("Look inside the correlations module loop");
    
    
    
    
    
%    ccn(nn+1,6)=sqrt(1-cn(nn+1,6).^2);
%    %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
%end

cnSqrSum=0;
for nn=1:NDims
%    ccn0(nn+1,6)=cn(nn+1,6).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,6).*cn(nn+1,6);
    for mm=1:NDims
        if(nn~=mm)
          %  ccn0(nn+1,6)=ccn0(nn+1,6)+2*cn(nn+1,6).*cn(mm+1,6).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,6).*cn(mm+1,6).*rhoxx(nn,mm,6);
        end
    end
    %ccn(nn+1,6)=sqrt(1-ccn0(nn+1,6));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,6)=cn(2:NDims+1,6)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end



cn(1,6)=sqrt(1-cnSqrSum);
%ccn(1,6)=sqrt(1-cn(1,6).^2);
%Cn(1,6)=cn(1,6).*YCoeffZ(6);

cn
%ccn
str=input("Compare the Sixth Element From NDim---2");



% cnSqrSum=0;
% for nn=1:NDims
%     CC(3,5)=YCoeffZ(7).*(15*cn(nn+1,7).^4.*ccn(nn+1,7).^2);
%     rho_z(nn+1,4)=(c(nn+1,5)-CC(3,5))./YCoeffZ(5);
%     cn(nn+1,5)=(abs(rho_z(nn+1,4))).^(1/4);
%     cn(nn+1,5)=cn(nn+1,5).*sign(rho_z(nn+1,4));
%     %ccn(nn+1,5)=sqrt(1-cn(nn+1,5).^2);
%     %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
% end

cnSqrSum=0;
for nn=1:NDims
    %ccn0(nn+1,5)=cn(nn+1,5).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,5).*cn(nn+1,5);
    for mm=1:NDims
        if(nn~=mm)
           % ccn0(nn+1,5)=ccn0(nn+1,5)+2*cn(nn+1,5).*cn(mm+1,5).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,5).*cn(mm+1,5).*rhoxx(nn,mm,5);
        end
    end
    %ccn(nn+1,5)=sqrt(1-ccn0(nn+1,5));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,5)=cn(2:NDims+1,5)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end


cn(1,5)=sqrt(1-cnSqrSum);
%ccn(1,5)=sqrt(1-cn(1,5).^2);
%Cn(1,5)=cn(1,5).*YCoeffZ(5);


cn
%ccn
str=input("Compare the Fifth Element From NDim---2");




%CC(3,4)=YCoeffZ(6).*(10*cn(6).^3.*dn(6).^2*zeta^2);
%CC(5,4)=YCoeffZ(8).*(35*cn(8).^3.*dn(8).^4.*zeta^4);
% cnSqrSum=0;
% for nn=1:NDims
%     CC(5,4)=YCoeffZ(8).*(35*cn(nn+1,8).^3.*ccn(nn+1,8).^4);
%     CC(3,4)=YCoeffZ(6).*(10*cn(nn+1,6).^3.*ccn(nn+1,6).^2);
%     rho_z(nn+1,3)=(c(nn+1,4)-3*CC(5,4)-CC(3,4))./YCoeffZ(4);
%     cn(nn+1,4)=(abs(rho_z(nn+1,3))).^(1/3);
%     cn(nn+1,4)=cn(nn+1,4).*sign(rho_z(nn+1,3));
%     %ccn(nn+1,4)=sqrt(1-cn(nn+1,4).^2);
%     %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
% end

cnSqrSum=0;
for nn=1:NDims
%    ccn0(nn+1,4)=cn(nn+1,4).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,4).*cn(nn+1,4);
    for mm=1:NDims
        if(nn~=mm)
          %  ccn0(nn+1,4)=ccn0(nn+1,4)+2*cn(nn+1,4).*cn(mm+1,4).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,4).*cn(mm+1,4).*rhoxx(nn,mm,4);
        end
    end
    %ccn(nn+1,4)=sqrt(1-ccn0(nn+1,4));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,4)=cn(2:NDims+1,4)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end

cn(1,4)=sqrt(1-cnSqrSum);
%ccn(1,4)=sqrt(1-cn(1,4).^2);
%Cn(1,4)=cn(1,4).*YCoeffZ(4);

cn
%ccn
str=input("Compare the Fourth Element From NDim---2");



%CC(3,3)=YCoeffZ(5).*(6* cn(5).^2.* dn(5).^2);
%CC(5,3)=YCoeffZ(7).*(15*cn(7).^2.*dn(7).^4);
% cnSqrSum=0;
% for nn=1:NDims
%     CC(3,3)=YCoeffZ(5).*(6*cn(nn+1,5).^2.*ccn(nn+1,5).^2);
%     CC(5,3)=YCoeffZ(7).*(15*cn(nn+1,7).^2.*ccn(nn+1,7).^4);
%     rho_z(nn+1,2)=(c(nn+1,3)-3*CC(5,3)-CC(3,3))./YCoeffZ(3);
%     cn(nn+1,3)=(abs(rho_z(nn+1,2))).^(1/2);
%     cn(nn+1,3)=cn(nn+1,3).*sign(rho_z(nn+1,2));
%     %ccn(nn+1,3)=sqrt(1-cn(nn+1,3).^2);
%     %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
% end

cnSqrSum=0;
for nn=1:NDims
    %ccn0(nn+1,3)=cn(nn+1,3).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,3).*cn(nn+1,3);
    for mm=1:NDims
        if(nn~=mm)
           % ccn0(nn+1,3)=ccn0(nn+1,3)+2*cn(nn+1,3).*cn(mm+1,3).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,3).*cn(mm+1,3).*rhoxx(nn,mm,3);
        end
    end
    %ccn(nn+1,3)=sqrt(1-ccn0(nn+1,3));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,3)=cn(2:NDims+1,3)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end

cn(1,3)=sqrt(1-cnSqrSum);
%ccn(1,3)=sqrt(1-cn(1,3).^2);
%Cn(1,3)=cn(1,3).*YCoeffZ(3);

cn
%ccn
str=input("Compare the Third Element From NDim---2");



%CC(7,2)=YCoeffZ(8).*(7*cn(8).*dn(8).^6*zeta^6);
%CC(5,2)=YCoeffZ(6).*(5*cn(6).*dn(6).^4.*zeta^4);
%CC(3,2)=YCoeffZ(4).* (3.* cn(4).* dn(4).^2.* zeta^2);
% cnSqrSum=0;
% for nn=1:NDims
%     CC(4,2)=YCoeffZ(4).*(3*cn(nn+1,4).*ccn(nn+1,4).^2);
%     CC(6,2)=YCoeffZ(6).*(5*cn(nn+1,5).*ccn(nn+1,6).^4);
%     CC(8,2)=YCoeffZ(8).*(7*cn(nn+1,8).*ccn(nn+1,8).^6);
%     rho_z(nn+1,1)=(c(nn+1,2)-15*CC(8,2)-3*CC(6,2)-CC(4,2))./YCoeffZ(2);
%     cn(nn+1,2)=(abs(rho_z(nn+1,1)));
%     %ccn(nn+1,2)=sqrt(1-cn(nn+1,2).^2);
%     %cnSqrSum=cnSqrSum+cn(7,nn+1).^2;
% end

cnSqrSum=0;
for nn=1:NDims
%    ccn0(nn+1,2)=cn(nn+1,2).^2;
    cnSqrSum=cnSqrSum+cn(nn+1,2).*cn(nn+1,2);
    for mm=1:NDims
        if(nn~=mm)
        %    ccn0(nn+1,2)=ccn0(nn+1,2)+2*cn(nn+1,2).*cn(mm+1,2).*CorrXX1st(nn,mm);
            cnSqrSum=cnSqrSum+cn(nn+1,2).*cn(mm+1,2).*rhoxx(nn,mm,2);
        end
    end
%    ccn(nn+1,2)=sqrt(1-ccn0(nn+1,2));
    
end

if(cnSqrSum>.995)
    cn(2:NDims+1,2)=cn(2:NDims+1,2)*sqrt(.995)/sqrt(cnSqrSum);
    cnSqrSum=.995;
end

if(cnSqrSum<.01)
    cn(2:NDims+1,2)=cn(2:NDims+1,2)*sqrt(.01)/sqrt(abs(cnSqrSum))/sign(cnSqrSum);
    cnSqrSum=.01;
end


cn(1,2)=sqrt(1-cnSqrSum);
%ccn(1,2)=sqrt(1-cn(1,2).^2);
%Cn(1,2)=cn(1,2).*YCoeffZ(2);

cn
%ccn
str=input("Compare the Second Element From NDim---2");


% 
% %CC(3,1)=YCoeffZ(2).*(ccn(nn+1,2).^2);
% %CC(5,1)=YCoeffZ(4).*(ccn(nn+1,4).^4);
% %CC(7,1)=YCoeffZ(6).*(ccn(nn+1,6).^6);
% CC(3,1)=YCoeffZ(2).*(cn(1,2).^2);
% CC(5,1)=YCoeffZ(4).*(cn(1,4).^4);
% CC(7,1)=YCoeffZ(6).*(cn(1,6).^6);
% 
% 
% 
% 
% 
% cn(1,1)=(YCoeffZ(1)-CC(3,1)-3*CC(5,1)-15*CC(7,1))./YCoeffZ(1);
% 
% 
% cn
% ccn
% str=input("Look at coefficients cn----PleASE lOOK");
% 
% 
 dn(1:Order+1)=cn(1,1:Order+1);
% bn(1:Order+1)=ccn(1,1:Order+1);
% 
% 
% 
% 

end

.
.
.
.

function [cnOptim] = OptimizeCorrFromHermitesNewton(Zy,Zx,cnn,HCorrXX,NDims,paths,Order)

%Order;Max Order of Hermites Calibration is seven.
%Zx is paths X K array for joint values of Zx arranged for "paths" number of observations
%Zy is paths X 1 array for Zy associated with each obsewrvation.
%cnn is the original input of correlation coefficients used to derive seed for Newton methjod.
%cnn is in format of the old program and its order and indexing is changed
%in this NEwton program.

ccn(1:NDims+1,1:Order)=cnn(1:NDims+1,2:Order+1); %Shift the indices back by one index.
%so that index matches with order of associated hermite polynomial.

cn(1:NDims,1:Order)=ccn(2:NDims+1,1:Order);%The first row in each hermite column earlier was 
%associated zeta but zeta is shifted from first row to (K+1)th row.

cn(NDims+1,1:Order)=ccn(1,1:Order);

% cn(1:NDims+1,1:Order)=1/(NDims+1);
% 
% for nn=1:NDims
%     cn(nn,1:Order)=cn(nn,1:Order)+sum(HCorrXX(nn,:,1));
% end


%The function below calculates the first hermites till "Order" from data of
%Zy. WHy is (paths X K) array. kth hermite associated with pth observation appears in
%WHy(p,k)
[WHy] = EvalHermiteArrayH1(Zy,Order);

%The program below calculates moments of standard normal. Odd moments are
%zero. 
[EZ] = EZZ2();

%The loop below repeats itself for optimization of correlation coefficients 
%associated with highest power in each hermite. It repeatrs over all hermites one by one. 
for nH=1:Order
    cnIn(1:NDims+1,1:nH)=0.0;
    if(nH==1)
        cnIn(1:NDims+1,1)=cn(1:NDims+1,1);
    else
        %cnOptim is already optimized value of lower powers of Gaussians
        %calculated in earlier loops of optimization over smaller hermites.
        %This is already known from previous iterations of this loop.
        cnIn(1:NDims+1,1:nH-1)=cnOptim(1:NDims+1,1:nH-1);
        %The value of cnIn below is initial guess for correlation
        %coefficients associated with optimization of current hermite
        cnIn(1:NDims+1,nH)=cn(1:NDims+1,nH);
    end
    
    %Below nH-th hermite observation values are assigned to one dimensional array for optimization target in Newton. 
    WHyn(1:paths,1)=WHy(1:paths,nH);
    
    %Total iterations of Newton root search.
    Iterations=500;
    %Below is the function that actually does the Multi-dimensional Newton.
    [cnbest,InitialObj(nH),FinalObj(nH)] = OptimizenthHermiteGaussianWithNewton(cnIn,Zx,WHyn,HCorrXX,nH,EZ,NDims,paths, Iterations);

    %Below is optimum value of correlation coefficients associated with
    %nH-th hermite after optimization. This will be used by the program as
    %an input to optimization of higher hermites
    cnOptim(1:NDims+1,nH)=cnbest(1:NDims+1,1);
    
    nH
end

cnOptim
InitialObj %Obj function at start of optimization of correlation coefficients for nH-th hermite
FinalObj%Obj function after optimization of correlation coefficients for nH-th hermite
str=input("Look at initial and Final Objs");

end

.
.
.
.

function [bnbest,InitialObj,FinalObj] = OptimizenthHermiteGaussianWithNewton(cn,Zx,WHyn,HCorrXX,nH,EZ,NDims,paths,Iterations)

%Why0 in array of nth Hermite Values

KK=NDims;

bn(1:NDims+1,1)=cn(1:NDims+1,nH);
%Below, Calculate the model contribution of lower Gaussian powers that have
%already been calculated through associated lower hermites towards current
%hermite polynomial.
[dHermite] = EvaluatePrevGaussPowersForHCorrs(nH,cn,Zx,NDims,paths,EZ);
%Remove the model contribution of lower Gaussian powers from data value of
%current hermite.
WHy0=WHyn+dHermite;

%Below calcualte variance
var=0.0;
for kk1=1:KK
    for kk2=1:KK    
        var=var+bn(kk1,1).*bn(kk2,1).*sign(HCorrXX(kk1,kk2,nH)).*(abs(HCorrXX(kk1,kk2,nH)));%.^(1/nH);%Changed.
    end
end

var=var+bn(KK+1,1).^2;

bn

var
%Below, make sure that variance is one.
bn(:,1)=bn(:,1)/sqrt(var);

str=input("have reached 1st Stone")        


%Below calculate the combined value of K correlated variables as described
%on Wilmott post and Calculate G.
GG1(1:paths,1)=0.0;
for kk=1:KK
    GG1(1:paths,1)=GG1(1:paths,1)+bn(kk,1).*Zx(1:paths,kk);
end

Zeta0=bn(KK+1,1);
%Evaluate the nHth Gaussian power below.
[GaussPH,GaussPh] = EvaluateHermitePolyForGaussCorrs(GG1,Zeta0,EZ,nH,paths);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Starting ObjBest and calculations here
%cnbest here
bnbest(1:NDims+1,1)=bn(1:NDims+1);

ObjSqrdBestNewton=sum((GaussPH(:,1)-WHy0(:,1)).^2)/paths;

InitialObj=ObjSqrdBestNewton;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%Below, enter the Newton loop.

StallCount=0;
iter=0;
Mul=.1;
while ((iter < Iterations) &&(StallCount<=5))
    iter=iter+1;
    dGG1(1:paths,1:NDims)=0.0;
    for kk=1:KK
        dGG1(1:paths,kk)=Zx(1:paths,kk);
    end
    
    bn=bnbest;
    
    GG1(1:paths,1)=0.0;
    for kk=1:KK
        GG1(1:paths,1)=GG1(1:paths,1)+bn(kk,1).*Zx(1:paths,kk);
    end

    Zeta0=bn(KK+1,1);
    
    bbn=0;
    bbn(1:NDims,1)=bn(1:NDims,1);

    
    [GaussPH,GaussPh,DGaussPH] = EvaluateHermitePolyAndDerivsForGaussCorrs(GG1,dGG1,Zeta0,EZ,nH,NDims,paths);


    %size(DGaussPH')
    %size(DGaussPH)
    %size(GaussPH)
    %size(WHy0(:,1))
    
    %str=input("Look at sizes");
    
    %I am using a multiple of .1 in Newton equation below. this might have
    %to be further decreased for some cases. Below is the Newton
    %over-constrained iteration.
    bbn=bbn+Mul*((DGaussPH'*DGaussPH)\(DGaussPH'*(WHy0(:,1)-GaussPH)));
    
    %DGaussPH'*DGaussPH
    
    %str=input("Look at matrix product");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Midfunction ObjBest and calculations here
%cnbest here

var=0.0;
for kk1=1:KK
    for kk2=1:KK    
        var=var+bbn(kk1,1).*bbn(kk2,1).*sign(HCorrXX(kk1,kk2,nH)).*(abs(HCorrXX(kk1,kk2,nH)));%.^(1/nH);%changed
    end
end

if(var>.9925)
    bbn(:,1)=bbn(:,1)/sqrt(var)*sqrt(.9925);
    var=.9925;
end
    

bbn(KK+1,1)=sqrt(1-var);

var=var+bbn(KK+1,1).^2;
bn
bbn

var

bbn=bbn/sqrt(var);



GG1(1:paths,1)=0.0;
for kk=1:KK
    GG1(1:paths,1)=GG1(1:paths,1)+bbn(kk,1).*Zx(1:paths,kk);
end


Zeta0=bbn(KK+1,1);

[GaussPH,GaussPh] = EvaluateHermitePolyForGaussCorrs(GG1,Zeta0,EZ,nH,paths);

ObjSqrdNewton=sum((GaussPH(:,1)-WHy0(:,1)).^2)/paths;
ObjSqrdBestNewton
if(ObjSqrdNewton<ObjSqrdBestNewton)
    if((ObjSqrdBestNewton-ObjSqrdNewton)/ObjSqrdBestNewton<.00001)
        StallCount=StallCount+1;
        Mul=Mul*.75
        
    end
        
    ObjSqrdBestNewton=ObjSqrdNewton;
    bnbest=bbn;
    
else
    
    StallCount=StallCount+1;
    %%%%change-step or Newton parameters
end
    
ObjSqrdBestNewton
ObjSqrdNewton
bnbest
iter

%str=input("Look at the solution");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5    

end

FinalObj=ObjSqrdBestNewton;

end

.
.
.

Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation

Re: Breakthrough in the theory of stochastic differential equations and their simulation