Friends, here is the new preliminary version of the program. I will post a new version with Newton root search added to it and other comprehensive additions in another two days. But here we are for now. Good thing is that density always remains stable and tail is always good.
.
.
function []= ConditionalDensityHermitesBivariateNewtonItersRegSqr02()
%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
Order=5;
%NDim=4;%Three assets and one SV.
NMomentsY=Order+1;
NMomentsX=Order+1;
w2D0(1:NMomentsY,1:NMomentsX)=1;
%%%%%%%%%%%%%%%%%%%5
V0=1.00;%.32;
V00=V0;
thetaV=1.050;%.045;%1;%.04;
kappaV=.320;%1.5;%1.5;
mu1=kappaV*thetaV;
mu2=-kappaV;
beta1=0;
beta2=1;
%gamma=.5;%.950;
%sigma0=.55;%.45;%
gamma=.55;%.950;
sigma0=.65;%.45;%
dt=.03125/2;
Tt=64*2;
T=Tt*dt;
seed0=52130649;
seed0=94210876;
rng(seed0, 'twister')
paths=10000;
V(1:paths,1)=V0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
TT1=30; %Transition distribution start
TT2=124; %Transition distribution end
Random2(1:paths/2,1)=0;
for tt=1:Tt
Random2(1:paths/2)=randn(paths/2,1);
Random2(paths/2+1:paths)=-Random2(1:paths/2);
V(1:paths,1)=V(1:paths,1)+ ...
(mu1.*V(1:paths,1).^beta1 + mu2.*V(1:paths,1).^beta2)*dt + ...
sigma0*V(1:paths,1).^gamma .*Random2(1:paths,1)*sqrt(dt) + ...
(mu1.*beta1*V(1:paths,1).^(beta1-1) + mu2.*beta2.*V(1:paths,1).^(beta2-1)).* ...
((mu1.*V(1:paths,1).^beta1 + mu2.*V(1:paths,1).^beta2)*dt^2/2 + ...
sigma0*V(1:paths,1).^gamma .*Random2(1:paths,1)*(1-1/sqrt(3))*dt^1.5) + ...
.5*(mu1.*beta1.*(beta1-1).*V(1:paths,1).^(beta1-2) + mu2.*beta2.*(beta2-1).*V(1:paths,1).^(beta2-2)).* ...
sigma0^2.*V(1:paths,1).^(2*gamma).*dt^2/2 + ...
sigma0*gamma*V(1:paths,1).^(gamma-1) .* ...
((mu1.*V(1:paths,1).^beta1 + mu2.*V(1:paths,1).^beta2).*Random2(1:paths,1).*1/sqrt(3)*dt^1.5 + ...
sigma0.*V(1:paths,1).^gamma .*(Random2(1:paths,1).^2-1)*dt/2) + ...
.5*sigma0*gamma*(gamma-1).*V(1:paths,1).^(gamma-2) .* ...
sigma0^2.*V(1:paths,1).^(2*gamma) .*Random2(1:paths,1).*1/sqrt(3)*dt^1.5;
V(V<0)=.00001;
if(tt==TT1)
Xin(1:paths,1)=V(1:paths,1);
end
if(tt==TT2)
Yin(1:paths,1)=V(1:paths,1);
end
end
str=input("I have reached stone 1");
Order=5;
%
Ndata=paths;
%
[ZI] = MomentMatchedStandardNormalDiscretizedDensity(Ndata);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%55
[XCoeffH,Zx] = CalculateHermiteSeriesFromData01(Xin,Order,ZI);
[XCoeffZ(1:Order+1)] = ConvertHermiteSeriesToZSeries01(XCoeffH(1:Order+1),Order);
for qq1=1:Order+1
XMoments(qq1)=sum(Xin(1:paths).^qq1)/paths;
end
[XCoeffZ(1),XCoeffZ(2:6)] = CalculateZSeriesDensityFromRawMomentsM6(XMoments,XCoeffZ(2:6));
for nn=1:paths
[Zx(nn)]=CalculateZgivenXAndZSeriesBisection2C5Improved(Xin(nn),XCoeffZ(1),XCoeffZ(2:6));
end
Zx0(1:Ndata,1)=Zx(1:Ndata);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[YCoeffH,Zy] = CalculateHermiteSeriesFromData01(Yin,Order,ZI);
[YCoeffZ(1:Order+1)] = ConvertHermiteSeriesToZSeries01(YCoeffH(1:Order+1),Order);
for pp1=1:Order+1
YMoments(pp1)=sum(Yin(1:paths).^pp1)/paths;
end
[YCoeffZ(1),YCoeffZ(2:6)] = CalculateZSeriesDensityFromRawMomentsM6(YMoments,YCoeffZ(2:6));
[YCoeffH0] = ConvertZSeriesToHermiteSeriesNew(YCoeffZ,Order);
for nn=1:paths
[Zy(nn)]=CalculateZgivenXAndZSeriesBisection2C5Improved(Yin(nn),YCoeffZ(1),YCoeffZ(2:6));
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
HH=5;
[CorrH] = CalculateCorrelationBivariateHermiteCH0(Zx,Zy,paths,HH)
CorrH
str=input("Look at hermites correlation");
Coeffyx2(1)=YCoeffH0(1);
Coeffyx2(2:6)=CorrH(1:5).*YCoeffH0(2:6);
Ydecorr(1:Ndata,1)=Yin(1:Ndata,1)-Coeffyx2(1) ...
-Coeffyx2(2).*Zx0(1:Ndata,1) ...
-Coeffyx2(3).*(Zx0(1:Ndata,1).^2-1) ...
-Coeffyx2(4).*(Zx0(1:Ndata,1).^3-3*Zx0(1:Ndata,1)) ...
-Coeffyx2(5).*(Zx0(1:Ndata,1).^4-6*Zx0(1:Ndata,1).^2+3) ...
-Coeffyx2(6).*(Zx0(1:Ndata,1).^5-10*Zx0(1:Ndata,1).^3+15*Zx0(1:Ndata,1));% ...
%Coeffyx2
%PlotHermiteSeriesDensityAndRvGraph(Coeffyx2(1),Coeffyx2(2:Order+1),'b')
%str=input("Look at initial correlated hermite series plot---2");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[YdCoeffH,Zyd] = CalculateHermiteSeriesFromData01(Ydecorr,Order,ZI);
[YdCoeffZ(1:Order+1)] = ConvertHermiteSeriesToZSeries01(YdCoeffH(1:Order+1),Order);
for pp1=1:Order+1
YdecorrMoments(pp1,1)=sum(Ydecorr(1:paths).^pp1)/paths;
end
[YdCoeffZ(1),YdCoeffZ(2:6)] = CalculateZSeriesDensityFromRawMomentsM6(YdecorrMoments,YdCoeffZ(2:6));
[YdCoeffH] = ConvertZSeriesToHermiteSeriesNew(YdCoeffZ,Order);
for nn=1:paths
[Zyd(nn)]=CalculateZgivenXAndZSeriesBisection2C5Improved(Ydecorr(nn),YdCoeffZ(1),YdCoeffZ(2:6));
end
Zyd0(1:paths,1)=Zyd(1:paths);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
YCoeffHH(1:6,1:6)=0.0;
YCoeffHH(2:6,1)=YdCoeffH(1,2:6);%*XCoeffRatio(1);
YCoeffHH
str=input("Look at 2D Hermite-series---5");
%%%%%%%%%%%%%%%%%%%%%%%%%%[YCoeffHH,ObjFunc] = HSeriesCoeffsFromCrossMomentsIterative2D2nd(YCoeffHH,XCoeffZ,MomentsYX0,YCoeffZ,MaxIter,w2D)
for pp1=1:Order+1
for qq1=1:Order+1
MomentsYX0((pp1-1)*(Order+1)+qq1,1)=sum(Ydecorr(1:paths,1).^pp1.*Xin(1:paths,1).^qq1)/paths;
%MomentsYX(pp1,qq1)=sum((1:paths,1).^pp1.*Zx0(1:paths,1).^qq1)/paths;
end
end
% MomentsYX
% str=input("Look at cross-moments matrix");
MaxIter=1500;
%In the function below, we do iterative optimization on all 2D Hermite
%coefficients excluding the first row. This function is applied when we
%have become too close to the true values by optimization through previous
%function.
[YCoeffHH,ObjFunc] = HSeriesCoeffsFromCrossMomentsIterative2D2nd04(YCoeffHH,XCoeffZ,MomentsYX0,YdCoeffZ,YdCoeffH,MaxIter,w2D0);
YCoeffHH
[YCoeffZZ] = Convert2DHermitesInto2DSeriesNew(YCoeffHH,Order,Order);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[YdecorrMomentsM] = CalculateMomentsOf2DZSeriesY(YCoeffZZ,Order,Order,Order+1);
YdecorrMomentsM
YdecorrMoments
str=input("Look at comparison of moments")
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[YCoeffHH] = ConvertZSeriesToHermiteSeries2D(YCoeffZZ,Order+1,Order+1);
YCoeffHH
str=input("Look at output oof 2D Hermite-Series");
YCoeffHH(1,1:Order+1)=YCoeffHH(1,1:Order+1)*0+Coeffyx2(1,1:Order+1);
[YCoeffZZ] = Convert2DHermitesInto2DSeriesNew(YCoeffHH,Order,Order);
[Y1DMoments]=CalculateMomentsOfZSeries(YCoeffZ(1),YCoeffZ(2:6),5,6)
OrderY=Order;
OrderX=Order;
NMoments=6;
[YMoments] = CalculateMomentsOf2DZSeriesY(YCoeffZZ,OrderY,OrderX,NMoments);
%
%
Y1DMoments
YMoments
str=input("Look at comparison of second moments/moments");
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%In this block, we do conditional Monte Carlo of SDE which is the true
%numerical 1D conditional density of Y conditional on X having a given
%value. This value of X on which we are conditionin is defined by V10 below.
%We can alter V10 below and it will take out a 1D slice from the 2D density
%of Y conditional on X taking a specific value (here V10).
%Zx=2.0;
%V10=XCoeffZ(1)+XCoeffZ(2).*Zx+XCoeffZ(3).*Zx^2+XCoeffZ(4).*Zx^3+XCoeffZ(5).*Zx^4+XCoeffZ(6).*Zx^5;
%V10
%str=input("Look at V10");
V10=1.750;
V(1:paths)=V10;
Random2(1:paths,1)=0;
for tt=TT1+1:TT2
Random2(1:paths/2)=randn(paths/2,1);
Random2(paths/2+1:paths)=-Random2(1:paths/2);
V(1:paths,1)=V(1:paths,1)+ ...
(mu1.*V(1:paths,1).^beta1 + mu2.*V(1:paths,1).^beta2)*dt + ...
sigma0*V(1:paths,1).^gamma .*Random2(1:paths,1)*sqrt(dt) + ...
(mu1.*beta1*V(1:paths,1).^(beta1-1) + mu2.*beta2.*V(1:paths,1).^(beta2-1)).* ...
((mu1.*V(1:paths,1).^beta1 + mu2.*V(1:paths,1).^beta2)*dt^2/2 + ...
sigma0*V(1:paths,1).^gamma .*Random2(1:paths,1)*(1-1/sqrt(3))*dt^1.5) + ...
.5*(mu1.*beta1.*(beta1-1).*V(1:paths,1).^(beta1-2) + mu2.*beta2.*(beta2-1).*V(1:paths,1).^(beta2-2)).* ...
sigma0^2.*V(1:paths,1).^(2*gamma).*dt^2/2 + ...
sigma0*gamma*V(1:paths,1).^(gamma-1) .* ...
((mu1.*V(1:paths,1).^beta1 + mu2.*V(1:paths,1).^beta2).*Random2(1:paths,1).*1/sqrt(3)*dt^1.5 + ...
sigma0.*V(1:paths,1).^gamma .*(Random2(1:paths,1).^2-1)*dt/2) + ...
.5*sigma0*gamma*(gamma-1).*V(1:paths,1).^(gamma-2) .* ...
sigma0^2.*V(1:paths,1).^(2*gamma) .*Random2(1:paths,1).*1/sqrt(3)*dt^1.5;
V(V<0)=.00001;
end
%Below
NoOfBins=200;
MaxCutOff=10;
[XDensity,IndexOut,IndexMax] = MakeDensityFromSimulation_Infiniti_NEW(V,paths,NoOfBins,MaxCutOff );
clf;
plot(IndexOut(1:IndexMax),XDensity(1:IndexMax),'r');
hold on
for qq=1:6
DataMoments(qq)=sum(V(:).^qq)/paths;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%Below, we calculate the value of Zx0 corresponding to X0 (here V10). This
%comes from inverting the Z-series of X so that X=V10;
xOrder=Order;
yOrder=Order;
[Zx0] = CalculateZgivenXAndHSeriesCoeffs(V10,XCoeffH,xOrder);
%Below, we collapse two dimensional pdf of Y|X to one dimensional pdf given
%Zx i.e. this Y|X=X(Zx)
%Here we take out a 1D slice from the 2D density
%of Y conditional on X taking a specific value (here V10).
He(1)=1;
He(2)=Zx0;
He(3)=Zx0^2-1;
He(4)=Zx0^3-3*Zx0;
He(5)=Zx0^4-6*Zx0.^2+3;
He(6)=Zx0^5-10*Zx0.^3+15*Zx0;
He(7)=Zx0^6-15*Zx0.^4+45*Zx0.^2-15;
He(8)=Zx0^7-21*Zx0.^5+105*Zx0.^3-105*Zx0;
YCoeffh(1:Order+1)=0;
for hh=1:yOrder+1
for hh2=1:xOrder+1
YCoeffh(hh)=YCoeffh(hh)+YCoeffHH(hh,hh2).*He(hh2);
end
end
[YCoeffz(1:Order+1)] = ConvertHermiteSeriesToZSeries01(YCoeffh(1:Order+1),Order);
[YMomentsModel]=CalculateMomentsOfZSeries(YCoeffz(1),YCoeffz(2:6),5,6);
YMomentsModel
DataMoments
str=input("Look at comparison of moments");
PlotHermiteSeriesDensityAndRvGraph(YCoeffh(1),YCoeffh(2:yOrder+1),'b')
%hold on
%PlotHermiteSeriesDensityAndRvGraph(YCoeffh(1),YCoeffh1(2:yOrder+1),'g')
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({'Numerical Density','Analytic Density1'}, ...
'Location','northeast')
str=input("This is the first version without altering the moments of conditional density")
[YMoments0]=CalculateMomentsOfZSeries(YCoeffz(1),YCoeffz(2:6),5,6)
end
.
.
.
function [YCoeffHH,ObjFunc] = HSeriesCoeffsFromCrossMomentsIterative2D2nd04(YCoeffHH,XCoeffZ,MomentsYX0,YdCoeffZ,YdCoeffH,MaxIter,w2D)
SeriesOrder=6;
NMoments=6;
Order=5;
VarY(2:Order+1)=0.0;
for nn=2:Order+1
VarY(nn)=YdCoeffH(nn).^2*factorial(nn-1);
end
%for nn=2:Order+1
% VarY(nn)=YCoeffHH(nn,1).^2+YCoeffHH(nn,2).^2+YCoeffHH(nn,2).^2*2+YCoeffHH(nn,4).^2*6+YCoeffHH(nn,5).^2*24+YCoeffHH(nn,6).^2*120;
% if(VarY(nn)>1.15^2*YdCoeffH(nn).^2)
% YCoeffHH(nn,1:Order+1)=YCoeffHH(nn,1:Order+1).*sqrt(YdCoeffH(nn).^2./VarY(nn))*1.15;
% end
%end
[YdMoments]=CalculateMomentsOfZSeries(YdCoeffZ(1),YdCoeffZ(2:6),5,6);
% %Below, we calculate first and 2nd moment of 2D Hermite series
% [MeanHH,M2HH] = CalculateMeanAndVarOf2DHSeriesY(YCoeffHH,Order,Order);
% %VarHH is model variance .
% VarHH=M2HH-MeanHH.^2;
% %VarIn is input variance
% VarIn=YMoments(2)-YMoments(1).^2;
% %Variance correction is applied below to all rows but not to first row of hermite series.
% %YCoeffHH(2:6,1:6)=YCoeffHH(2:6,1:6).*sqrt(VarIn/VarHH);
% %YCoeffHH(1,2:6)=YCoeffHH(1,2:6).*sqrt(VarIn/VarHH);
% %We convert final 2D Hermite series to 2D Z-series
[YCoeffZZ] = Convert2DHermitesInto2DSeriesNew(YCoeffHH,Order,Order);
for nn=1:SeriesOrder
for mm=1:SeriesOrder
DeltaX(nn,mm)=.00125/(factorial(nn+1))/factorial(mm+1);%/2^((nn-1)/2)/2^((mm-1)/2);
DeltaXMax(nn,mm)=DeltaX(nn,mm)*10;
end
end
ObjFunc1=0;
ObjFunc2=0;
%%%Weights on objective function that shows fit to moments
w2D(NMoments,NMoments)=1;
for nn=NMoments-1:-1:1
for mm=NMoments-1:-1:1
w2D(nn,mm)=w2D(nn+1,mm+1)*(nn+1).^0*(mm+1).^0;
end
end
[AMat] = CalculateCrossMomentsofTwo2Dand1DZSeries(YCoeffZZ,XCoeffZ,Order,Order);
%[BMat] = CalculateMomentsOf2DZSeriesY(YCoeffZZ,Order,Order,Order+1);
for pp1=1:Order+1
for qq1=1:Order+1
Moments2D(pp1,qq1)=AMat((pp1-1)*(Order+1)+qq1,1);
Moments2D0(pp1,qq1)=MomentsYX0((pp1-1)*(Order+1)+qq1,1);
end
% Moments1D(pp1)=AMat((Order+1)*(Order+1)+pp1,1);
end
[ObjFuncBest] = CalculateObjMoment2D1D(Moments2D0,Moments2D,w2D,Order+1,Order+1);
%ObjFuncBest=ObjFuncBest+sum((BMat-YdMoments).^2);
%[ObjFuncBest] = CalculateObjMoment(sMu,Moments,w,mOrder);
ObjFuncBest
str=input("Look objFuncBest");
%cBest=c;
ImproveFlag(1:SeriesOrder,1:SeriesOrder)=1;
ImproveFlagPrev(1:SeriesOrder,1:SeriesOrder)=0;
ImproveCount(1:SeriesOrder,1:SeriesOrder)=0;
SeriesOrderM=6;
TooSmallStep=0;
ObjFuncBestBefore=ObjFuncBest;
kk=0;
while ((kk<MaxIter)&&(ObjFuncBest>.001)&&(TooSmallStep<100))
kk=kk+1;
if(ObjFuncBestBefore-ObjFuncBest<.000001)
TooSmallStep=TooSmallStep+1;
else
TooSmallStep=0;
end
ObjFuncBestBefore=ObjFuncBest;
if(kk>400)
SeriesOrderM=6;
end
if(kk>700)
SeriesOrderM=6;
end
%We only iteratively alter c(2:7), first and second moments are
%retrieved automatically in our set up
for nn=2:SeriesOrder
for mm=1:SeriesOrderM
if((nn==1))
% ;
else
ImproveFlagPrev(nn,mm)=ImproveFlag(nn,mm);
YCoeffHHnew=YCoeffHH;
% [YCoeffZZnew] = Convert2DHermitesInto2DSeriesNew(YCoeffHHnew,Order,Order);
%
% [AMat] = CalculateCrossMomentsofTwo2Dand1DZSeries(YCoeffZZnew,XCoeffZ,Order,Order);
% for pp1=1:Order+1
% for qq1=1:Order+1
% Moments2D(pp1,qq1)=AMat((pp1-1)*(Order+1)+qq1,1);
% Moments2D0(pp1,qq1)=MomentsYX0((pp1-1)*(Order+1)+qq1,1);
% end
% end
CheckObjFlag=0;
if(Moments2D0(nn,mm)>Moments2D(nn,mm))
CheckObjFlag=1;
YCoeffHHnew(nn,mm)=YCoeffHHnew(nn,mm)+DeltaX(nn,mm);
YSD(nn)=sqrt(YCoeffHHnew(nn,1).^2+YCoeffHHnew(nn,2).^2+YCoeffHHnew(nn,3).^2*2+YCoeffHHnew(nn,4).^2*6+YCoeffHHnew(nn,5).^2*24+YCoeffHHnew(nn,6).^2*120);
if(kk>-50)
if(YSD(nn)*1.0>YdCoeffH(nn))
YCoeffHHnew(nn,:)=YCoeffHHnew(nn,:).*YdCoeffH(nn)./YSD(nn)/1.0;
end
end %
% if(nn>1)
% VarY(nn)=YCoeffHHnew(nn,1).^2+YCoeffHHnew(nn,2).^2+YCoeffHHnew(nn,2).^2*2+YCoeffHHnew(nn,4).^2*6+YCoeffHHnew(nn,5).^2*24+YCoeffHHnew(nn,6).^2*120;
% if(VarY(nn)>1.15^2*YdCoeffH(nn).^2)
% YCoeffHHnew(nn,1:Order+1)=YCoeffHHnew(nn,1:Order+1).*sqrt(YdCoeffH(nn).^2./VarY(nn))*1.15;
% end
% end
%Below, we calculate first and 2nd moment of 2D Hermite series
% [MeanHH,M2HH] = CalculateMeanAndVarOf2DHSeriesY(YCoeffHHnew,Order,Order);
%VarHH is model variance .
% VarHH=M2HH-MeanHH.^2;
%VarIn is input variance
% VarIn=YMoments(2)-YMoments(1).^2;
%Variance correction is applied below to all rows but not to first row of hermite series.
% YCoeffHHnew(2:6,1:6)=YCoeffHHnew(2:6,1:6).*sqrt(VarIn/VarHH);
% YCoeffHHnew(1,2:6)=YCoeffHHnew(1,2:6).*sqrt(VarIn/VarHH);
%We convert final 2D Hermite series to 2D Z-series
[YCoeffZZnew] = Convert2DHermitesInto2DSeriesNew(YCoeffHHnew,Order,Order);
[AMat] = CalculateCrossMomentsofTwo2Dand1DZSeries(YCoeffZZnew,XCoeffZ,Order,Order);
for pp1=1:Order+1
for qq1=1:Order+1
Moments2D(pp1,qq1)=AMat((pp1-1)*(Order+1)+qq1,1);
Moments2D0(pp1,qq1)=MomentsYX0((pp1-1)*(Order+1)+qq1,1);
end
end
end
[ObjFunc1] = CalculateObjMoment2D1D(Moments2D0,Moments2D,w2D,Order+1,Order+1);
%[BMat] = CalculateMomentsOf2DZSeriesY(YCoeffZZnew,Order,Order,Order+1);
%ObjFunc1=ObjFunc1+sum((BMat-YdMoments).^2);
if((ObjFunc1<ObjFuncBest)&&(CheckObjFlag==1))
ObjFuncBest=ObjFunc1;
YCoeffHH=YCoeffHHnew;
ImproveFlag(nn,mm)=1;
else
YCoeffHHnew=YCoeffHH;
% [YCoeffZZnew] = Convert2DHermitesInto2DSeriesNew(YCoeffHHnew,Order,Order);
%
%
% [AMat] = CalculateCrossMomentsofTwo2Dand1DZSeries(YCoeffZZnew,XCoeffZ,Order,Order);
% for pp1=1:Order+1
% for qq1=1:Order+1
% Moments2D(pp1,qq1)=AMat((pp1-1)*(Order+1)+qq1,1);
% Moments2D0(pp1,qq1)=MomentsYX0((pp1-1)*(Order+1)+qq1,1);
% end
% end
CheckObjFlag=0;
if(Moments2D0(nn,mm)<Moments2D(nn,mm))
CheckObjFlag=1;
YCoeffHHnew(nn,mm)=YCoeffHHnew(nn,mm)-DeltaX(nn,mm);
%YCoeffHHnew(nn,1)=sqrt(YdCoeffH(nn).^2-(YCoeffHHnew(nn,2).^2+YCoeffHHnew(nn,3).^2*2+YCoeffHHnew(nn,4).^2*6+YCoeffHHnew(nn,5).^2*24*1.414+YCoeffHHnew(nn,6).^2*120*2));
YSD(nn)=sqrt(YCoeffHHnew(nn,1).^2+YCoeffHHnew(nn,2).^2+YCoeffHHnew(nn,3).^2*2+YCoeffHHnew(nn,4).^2*6+YCoeffHHnew(nn,5).^2*24+YCoeffHHnew(nn,6).^2*120);
if(kk>-50)
if(YSD(nn)*1.0>YdCoeffH(nn))
YCoeffHHnew(nn,:)=YCoeffHHnew(nn,:).*YdCoeffH(nn)./YSD(nn)/1.0;
end
end
%
% if(nn>1)
% VarY(nn)=YCoeffHHnew(nn,1).^2+YCoeffHHnew(nn,2).^2+YCoeffHHnew(nn,2).^2*2+YCoeffHHnew(nn,4).^2*6+YCoeffHHnew(nn,5).^2*24+YCoeffHHnew(nn,6).^2*120;
% if(VarY(nn)>1.15^2*YdCoeffH(nn).^2)
% YCoeffHHnew(nn,1:Order+1)=YCoeffHHnew(nn,1:Order+1).*sqrt(YdCoeffH(nn).^2./VarY(nn))*1.15;
% end
% end
%Below, we calculate first and 2nd moment of 2D Hermite series
% [MeanHH,M2HH] = CalculateMeanAndVarOf2DHSeriesY(YCoeffHHnew,Order,Order);
%VarHH is model variance .
% VarHH=M2HH-MeanHH.^2;
%VarIn is input variance
% VarIn=YMoments(2)-YMoments(1).^2;
%Variance correction is applied below to all rows but not to first row of hermite series.
% YCoeffHHnew(2:6,1:6)=YCoeffHHnew(2:6,1:6).*sqrt(VarIn/VarHH);
% YCoeffHHnew(1,2:6)=YCoeffHHnew(1,2:6).*sqrt(VarIn/VarHH);
%We convert final 2D Hermite series to 2D Z-series
[YCoeffZZnew] = Convert2DHermitesInto2DSeriesNew(YCoeffHHnew,Order,Order);
[AMat] = CalculateCrossMomentsofTwo2Dand1DZSeries(YCoeffZZnew,XCoeffZ,Order,Order);
for pp1=1:Order+1
for qq1=1:Order+1
Moments2D(pp1,qq1)=AMat((pp1-1)*(Order+1)+qq1,1);
Moments2D0(pp1,qq1)=MomentsYX0((pp1-1)*(Order+1)+qq1,1);
end
end
end
[ObjFunc2] = CalculateObjMoment2D1D(Moments2D0,Moments2D,w2D,Order+1,Order+1);
% [BMat] = CalculateMomentsOf2DZSeriesY(YCoeffZZnew,Order,Order,Order+1);
% ObjFunc2=ObjFunc2+sum((BMat-YdMoments).^2);
if((ObjFunc2<ObjFuncBest)&&(CheckObjFlag==1))
ObjFuncBest=ObjFunc2;
YCoeffHH=YCoeffHHnew;
ImproveFlag(nn,mm)=-1;
else
ImproveFlag(nn,mm)=0;
end
end
end
end
end
if(rem(kk,1)==0)
kk
ObjFunc1
ObjFunc2
ObjFuncBest
DeltaX
ImproveCount
end
ImproveFlag0=0;
for nn=1:SeriesOrder
for mm=1:SeriesOrder
if(ImproveFlag(nn,mm)==0)
DeltaX(nn,mm)=DeltaX(nn,mm)*.88;
%ImproveFlag0=1;
end
if((ImproveFlag(nn,mm)==1)||(ImproveFlag(nn,mm)==-1))
DeltaX(nn,mm)=DeltaX(nn,mm)*1.1;
%ImproveFlag0=1;
end
if(DeltaX(nn,mm)>DeltaXMax(nn,mm))
DeltaX(nn,mm)=DeltaXMax(nn,mm);
end
end
end
end
ObjFunc=ObjFuncBest;
end
.
.
.
Here are the last parts of the program that would appear on screen.
Look at output oof 2D Hermite-Series
Y1DMoments =
Columns 1 through 3
1.022321303494130 1.529211477201552 3.014113407951637
Columns 4 through 6
7.394867276878588 21.734999830747970 74.292100365702012
Y1DMoments =
Columns 1 through 3
1.022321303494130 1.529211477201552 3.014113407951637
Columns 4 through 6
7.394867276878588 21.734999830747970 74.292100365702012
YMoments =
Columns 1 through 3
1.022321303494130 1.524394473942560 2.959964048421259
Columns 4 through 6
7.212836423341501 21.643914852332809 79.437115792148774
Look at comparison of second moments/moments
IndexMax =
201
YMomentsModel =
1.0e+02 *
Columns 1 through 3
0.014584519586695 0.027165562397651 0.061791190407146
Columns 4 through 6
0.168383895264301 0.545006028422809 2.093815236555474
DataMoments =
1.0e+02 *
Columns 1 through 3
0.014846685485368 0.028652042693887 0.067506320477130
Columns 4 through 6
0.187343327031092 0.600194804917567 2.203990882949031
Here is the graph you will see.
.
.
.
