Breakthrough in the theory of stochastic differential equations and their simulation

Amin · Today, 7:26 am

EDIT: Please use the program given in post 2854. This(given in this post) is older version that has some errors. That(given in post 2854) is the new good version. Notes below are valid for that program in post 2854.

Friends, please use this version (still very raw) of the function code for good results. It works well for very high volatilities for cev .55, .65 or even .75. Mostly until 1D density of Ydecorr remains stable.

Start with zeros on first row of YCoeffHH(1,:) of Ydecorr (Decorrelated part of Y). (This first row will be later replaced with correlated part of Y from Zx)
Start with one dimensional Hermite-Series (YdCoeffH(2:6)) on first column of 2D Hermite-Series YCoeffHH(2:6,1) of Ydecorr (Decorrelated part of Y).
Start with zeros in the rest of the body of YCoeffHH of Ydecorr.

Please do not use any 2D least squares regression that we were previously using. Just start iterations from the initial guess as described in above few lines.

YdCoeffZ,YdCoeffH are arguments from one dimensional Z-series and Hermite-Series of Ydecorr.
MomentsYX0 are cross moments between Ydecorr and X.

Even if you have another good working program, be sure to include the following lines from this function. These lines mean that 2D hermite series variance(actually SD) of any row is not allowed to increase the 1D Variance of that hermite polynomial of Y from YdCoeffH. It is only allowed to increase by a very small factor of .005 as in 1.005.

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(YSD(nn)>YdCoeffH(nn)*1.005)
YCoeffHHnew(nn,:)=YCoeffHHnew(nn,:).*YdCoeffH(nn)./YSD(nn)*1.005;
end

.
.
.

function [YCoeffHH,ObjFunc] = HSeriesCoeffsFromCrossMomentsIterative2D2nd(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)=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)>YdCoeffH(nn)^2)
    YCoeffHH(nn,1:Order+1)=YCoeffHH(nn,1:Order+1).*sqrt(YdCoeffH(nn).^2./VarY(nn));
    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*2/(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;
kk=0;
while ((kk<MaxIter)&&(ObjFuncBest>.001))
kk=kk+1;


    %We only iteratively alter c(2:7), first and second moments are
    %retrieved automatically in our set up
    for nn=2:6%SeriesOrder
        for mm=1:6%SeriesOrder
            if((nn==1))
            %    ;
            else
            
            ImproveFlagPrev(nn,mm)=ImproveFlag(nn,mm);


            YCoeffHHnew=YCoeffHH;
            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(YSD(nn)>YdCoeffH(nn)*1.005)
                YCoeffHHnew(nn,:)=YCoeffHHnew(nn,:).*YdCoeffH(nn)./YSD(nn)*1.005;
            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
            [ObjFunc1] = CalculateObjMoment2D1D(Moments2D0,Moments2D,w2D,Order+1,Order+1);
            [BMat] = CalculateMomentsOf2DZSeriesY(YCoeffZZnew,Order,Order,Order+1);
            %ObjFunc1=ObjFunc1+sum((BMat-YdMoments).^2);
                 [DValue1] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,-3.0,Order);
                [DValue2] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,3.0,Order);
                [DValue3] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,-2.5,Order);
                [DValue4] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,2.5,Order);
                [DValue5] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,-2.0,Order);
                [DValue6] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,2.0,Order);
                
                DFlag1=(DValue1>0)&&(DValue2>0)&&(DValue3>0)&&(DValue4>0)&&(DValue5>0)&&(DValue6>0);
                [DValue1] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,-3.0,Order);
                [DValue2] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,3.0,Order);
                [DValue3] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,-2.5,Order);
                [DValue4] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,2.5,Order);
                [DValue5] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,-2.0,Order);
                [DValue6] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,2.0,Order);
                
                DFlag2=(DValue1>0)&&(DValue2>0)&&(DValue3>0)&&(DValue4>0)&&(DValue5>0)&&(DValue6>0);
 %DValue1=1;
 %DValue2=1;
 
 %DFlag1=1;
 %DFlag2=1;
                if((ObjFunc1<ObjFuncBest)&&(DFlag1==1)&&(DFlag2==1))       
            %if(ObjFunc1<ObjFuncBest)
            
                ObjFuncBest=ObjFunc1;
            
                YCoeffHH=YCoeffHHnew;
                ImproveFlag(nn,mm)=1;
            
            else
            
                YCoeffHHnew=YCoeffHH;
                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(YSD(nn)>YdCoeffH(nn)*1.005)
                YCoeffHHnew(nn,:)=YCoeffHHnew(nn,:).*YdCoeffH(nn)./YSD(nn)*1.005;
            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
                [ObjFunc2] = CalculateObjMoment2D1D(Moments2D0,Moments2D,w2D,Order+1,Order+1);
                
                [BMat] = CalculateMomentsOf2DZSeriesY(YCoeffZZnew,Order,Order,Order+1);
             %   ObjFunc2=ObjFunc2+sum((BMat-YdMoments).^2);
                
                
                                [DValue1] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,-3.0,Order);
                [DValue2] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,3.0,Order);
                [DValue3] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,-2.5,Order);
                [DValue4] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,2.5,Order);
                [DValue5] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,-2.0,Order);
                [DValue6] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,2.0,Order);
                
                DFlag1=(DValue1>0)&&(DValue2>0)&&(DValue3>0)&&(DValue4>0)&&(DValue5>0)&&(DValue6>0);
                [DValue1] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,-3.0,Order);
                [DValue2] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,3.0,Order);
                [DValue3] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,-2.5,Order);
                [DValue4] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,2.5,Order);
                [DValue5] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,-2.0,Order);
                [DValue6] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,2.0,Order);
                
                DFlag2=(DValue1>0)&&(DValue2>0)&&(DValue3>0)&&(DValue4>0)&&(DValue5>0)&&(DValue6>0);
 %DValue1=1;
 %DValue2=1;
 
 %DFlag1=1;
 %DFlag2=1;
                if((ObjFunc2<ObjFuncBest)&&(DFlag1==1)&&(DFlag2==1))
               % if(ObjFunc2<ObjFuncBest)

                    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

Amin · Today, 8:17 am

Friends, I would love to write a very comprehensive program as I had earlier suggested, but mind control agents started altering my programs again and that is why I could not work on it day before yesterday and yesterday. I want to try doing that again today and want to see if mind control agents are still adamant on mischief. I will let friends know in the evening how it goes.

Amin · Today, 8:51 am

I am not sure how this program works, but it was working superbly three days ago. Here it is
.
.
.

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;
 
 w2D5(1:NMomentsY,1:NMomentsX)=1;
 w2D6(1:NMomentsY,1:NMomentsX)=1;
for nn=NMomentsY-1:-1:1
    for mm=NMomentsX-1:-1:1
        
        w2D5(nn,mm)=w2D5(nn+1,mm+1)*(nn+1).^5*(mm+1).^5.0;
        w2D6(nn,mm)=w2D6(nn+1,mm+1)*(nn+1).^6*(mm+1).^6.0;
    end
end
 
 
 w2D0(1:NMomentsY,1:NMomentsX)=1;
 w2Dh(1:NMomentsY,1:NMomentsX)=1;
 w2D1(1:NMomentsY,1:NMomentsX)=1;
 w2D2(1:NMomentsY,1:NMomentsX)=1;
 w2D3(1:NMomentsY,1:NMomentsX)=1;
 w2D4(1:NMomentsY,1:NMomentsX)=1;
 
for nn=NMomentsY-1:-1:1
    for mm=NMomentsX-1:-1:1
        w2Dh(nn,mm)=w2Dh(nn+1,mm+1)*(nn+1).^.5*(mm+1).^.5;
        w2D1(nn,mm)=w2D1(nn+1,mm+1)*(nn+1)*(mm+1);
        w2D2(nn,mm)=w2D2(nn+1,mm+1)*(nn+1).^2*(mm+1).^2.0;
        w2D3(nn,mm)=w2D3(nn+1,mm+1)*(nn+1).^3*(mm+1).^3.0;
        w2D4(nn,mm)=w2D4(nn+1,mm+1)*(nn+1).^4*(mm+1).^4.0;
    end
end

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

V0=1.00;%.32;
V00=V0;
thetaV=1.050;%.045;%1;%.04;
kappaV=.32000;%1.5;%1.5;
mu1=kappaV*thetaV;
mu2=-kappaV;
beta1=0;
beta2=1;
%gamma=.5;%.950;
%sigma0=.55;%.45;%
gamma=.65;%.950;
sigma0=.55;%.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=60; %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));
%[YCoeffZ(1),YCoeffZ(2:6)] = CalculateZSeriesDensityFromRawMoments02(YMoments);


[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(7).*(Zx0(1:Ndata,1).^6-15*Zx0(1:Ndata,1).^4+45*Zx0(1:Ndata,1).^2-15) ...
    %-Coeffyx2(8).*(Zx0(1:Ndata,1).^7-21*Zx0(1:Ndata,1).^5+105*Zx0(1:Ndata,1).^3-105*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));
%[YdCoeffZ(1),YdCoeffZ(2:6)] = CalculateZSeriesDensityFromRawMoments02(YdecorrMoments);

[YdCoeffH] = ConvertZSeriesToHermiteSeriesNew(YdCoeffZ,Order);  

clf;

PlotHermiteSeriesDensityAndRvGraph(YdCoeffH(1),YdCoeffH(2:Order+1),'b')

str=input("Look at YdCoeffH density");



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;
          %MomentsYX0((pp1-1)*(Order+1)+qq1,1)=sum(Ydecorr(1:paths,1).^pp1.*Zx0(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=150; 
%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.

%ZCoeffZ(1:6)=0;
%ZCoeffZ(2)=1;

%[YCoeffHH,ObjFunc] = HSeriesCoeffsFromCrossMomentsIterative2D2nd(YCoeffHH,XCoeffZ,MomentsYX0,YdCoeffZ,YdCoeffH,MaxIter,w2D2);

%[YCoeffHH,ObjFunc] = HSeriesCoeffsFromCrossMomentsIterative2D2nd04(YCoeffHH,XCoeffZ,MomentsYX0,YdCoeffZ,YdCoeffH,MaxIter,w2D2);
%YCoeffHH

TargetCrossMoments=MomentsYX0;
OrderY=Order;
OrderX=Order;
%[YCoeffHH,ObjFunc] = HSeriesCoeffsFromCrossMomentsIterative2D2nd04(YCoeffHH,XCoeffZ,TargetCrossMoments,YdCoeffH,OrderY,OrderX,MaxIter);
MaxIter=300;
[YCoeffHH,ObjFunc] = HSeriesCoeffsFromCrossMomentsIterative2D2ndNew(YCoeffHH,XCoeffZ,MomentsYX0,YdCoeffZ,YdCoeffH,MaxIter,w2D2);


%load("C:\Users\Lenovo\Documents\MATLAB\MATLAB1\pqfile02A.mat","YCoeffHH")%45000


YCoeffHH
str=input("Variable has been saved")


%YCoeffHH(6,:)=YCoeffHH(6,:);


[YCoeffZZ] = Convert2DHermitesInto2DSeriesNew(YCoeffHH,Order,Order);


%load("C:\Users\Lenovo\Documents\MATLAB\MATLAB1\pqfile02.mat","YCoeffZZ")%45000


MaxIter=2500; 



%[YCoeffZZ,AMat] = DoNewtonIterationsOn2DZSeriesSubset(YCoeffZZ,XCoeffZ,Order,MomentsYX0,YdecorrMoments,YdCoeffH,MaxIter);
%[YCoeffZZ] = Convert2DHermitesInto2DSeriesNew(YCoeffHH,Order,Order);
%LG[YCoeffZZ,AMat] = DoNewtonIterationsOn2DZSeries(YCoeffZZ,XCoeffZ,Order,MomentsYX0,YdecorrMoments,MaxIter);
%[YCoeffZZ,AMat] = DoNewtonIterationsOn2DZSeriesWithData(YCoeffZZ,XCoeffZ,Order,MomentsYX0,YdecorrMoments,MaxIter);%,Zyd,Zx,Ndata);
%[YCoeffZZ,AMat] = DoNewtonIterationsOn2DZSeriesWithData(YCoeffZZ,ZCoeff,Order,MomentsYX0,MomentsY,MaxIter,Zyd,Zx,Ndata);

[YdecorrMomentsM] = CalculateMomentsOf2DZSeriesY(YCoeffZZ,Order,Order,Order+1);

YdecorrMomentsM
YdecorrMoments






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%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

.
.
.
If 1D density of Ydecorr given on line 201 of above function is well-behaved (with proper tails) with five hermite polynomials (six moments), the optimization should go well with five hermite polynomials. If this density is not well-behaved with five polynomials, the optimization would be hopeless with five polynomials. I checked that for most cases, where the density that is not well-behaved with five polynomials, becomes well-behaved with nice tails when we calculate Ydecorr with seven hermite polynomials(For that You have to re-do all the earlier calculations of 1D Densities, correlations and correlated part with seven hermite polynomials(eight moments). I will automate this in new version)

Amin · Today, 9:03 am

Friends, very sorry, please do not use the function given at the start of this page in post 2851. Rather you need to use its new version as given here. Notes there are the same and are valid for this program.
.
.

function [YCoeffHH,ObjFunc] = HSeriesCoeffsFromCrossMomentsIterative2D2ndNew(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)=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)>YdCoeffH(nn)^2)
    YCoeffHH(nn,1:Order+1)=YCoeffHH(nn,1:Order+1).*sqrt(YdCoeffH(nn).^2./VarY(nn));
    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*2/(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;
kk=0;
while ((kk<MaxIter)&&(ObjFuncBest>.001))
kk=kk+1;

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:6%SeriesOrder
        for mm=1:6%SeriesOrderM
            if((nn==1))
            %    ;
            else
            
            ImproveFlagPrev(nn,mm)=ImproveFlag(nn,mm);


            YCoeffHHnew=YCoeffHH;
            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(YSD(nn)>YdCoeffH(nn))
                YCoeffHHnew(nn,:)=YCoeffHHnew(nn,:).*YdCoeffH(nn)./YSD(nn);
            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
            [ObjFunc1] = CalculateObjMoment2D1D(Moments2D0,Moments2D,w2D,Order+1,Order+1);
            [BMat] = CalculateMomentsOf2DZSeriesY(YCoeffZZnew,Order,Order,Order+1);
            %ObjFunc1=ObjFunc1+sum((BMat-YdMoments).^2);
                 [DValue1] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,-3.0,Order);
                [DValue2] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,3.0,Order);
                [DValue3] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,-2.5,Order);
                [DValue4] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,2.5,Order);
                [DValue5] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,-2.0,Order);
                [DValue6] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,2.0,Order);
                
                DFlag1=(DValue1>0)&&(DValue2>0)&&(DValue3>0)&&(DValue4>0)&&(DValue5>0)&&(DValue6>0);
                [DValue1] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,-3.0,Order);
                [DValue2] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,3.0,Order);
                [DValue3] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,-2.5,Order);
                [DValue4] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,2.5,Order);
                [DValue5] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,-2.0,Order);
                [DValue6] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,2.0,Order);
                
                DFlag2=(DValue1>0)&&(DValue2>0)&&(DValue3>0)&&(DValue4>0)&&(DValue5>0)&&(DValue6>0);
 %DValue1=1;
 %DValue2=1;
 
 %DFlag1=1;
 %DFlag2=1;
                if((ObjFunc1<ObjFuncBest)&&(DFlag1==1)&&(DFlag2==1))       
            %if(ObjFunc1<ObjFuncBest)
            
                ObjFuncBest=ObjFunc1;
            
                YCoeffHH=YCoeffHHnew;
                ImproveFlag(nn,mm)=1;
            
            else
            
                YCoeffHHnew=YCoeffHH;
                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(YSD(nn)>YdCoeffH(nn))
                YCoeffHHnew(nn,:)=YCoeffHHnew(nn,:).*YdCoeffH(nn)./YSD(nn);
            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
                [ObjFunc2] = CalculateObjMoment2D1D(Moments2D0,Moments2D,w2D,Order+1,Order+1);
                
                [BMat] = CalculateMomentsOf2DZSeriesY(YCoeffZZnew,Order,Order,Order+1);
             %   ObjFunc2=ObjFunc2+sum((BMat-YdMoments).^2);
                
                
                                [DValue1] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,-3.0,Order);
                [DValue2] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,3.0,Order);
                [DValue3] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,-2.5,Order);
                [DValue4] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,2.5,Order);
                [DValue5] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,-2.0,Order);
                [DValue6] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,-3.0,2.0,Order);
                
                DFlag1=(DValue1>0)&&(DValue2>0)&&(DValue3>0)&&(DValue4>0)&&(DValue5>0)&&(DValue6>0);
                [DValue1] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,-3.0,Order);
                [DValue2] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,3.0,Order);
                [DValue3] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,-2.5,Order);
                [DValue4] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,2.5,Order);
                [DValue5] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,-2.0,Order);
                [DValue6] = EvaluateDerivativeOf2DZSeries(YCoeffZZnew,3.0,2.0,Order);
                
                DFlag2=(DValue1>0)&&(DValue2>0)&&(DValue3>0)&&(DValue4>0)&&(DValue5>0)&&(DValue6>0);
 %DValue1=1;
 %DValue2=1;
 
 %DFlag1=1;
 %DFlag2=1;
                if((ObjFunc2<ObjFuncBest)&&(DFlag1==1)&&(DFlag2==1))
               % if(ObjFunc2<ObjFuncBest)

                    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

Amin · Today, 9:10 am

Friends, you need to comment the calculation of BMat in the above program everywhere since it is calculated but not used at all. If you want to run the program as such, you will also need the following function.
.
.

function [YMoments] = CalculateMomentsOf2DZSeriesY(YCoeffZZ,OrderY,OrderX,NMoments)

[EZ] = EZZ2();



EX=0;
for nn1=1:OrderY
    for nn2=1:OrderX
        EX=EX+YCoeffZZ(nn1,nn2).*EZ(nn1).*EZ(nn2);
    end
end

MeanZZ=EX;
YMoments(1)=EX;



YpCoeff(1:(OrderY+1)*(OrderX+1),1:(OrderY+1)*(OrderX+1),1:NMoments)=0.0;

YpCoeff(1:OrderY+1,1:OrderX+1,1)=YCoeffZZ(1:OrderY+1,1:OrderX+1);

SeriesOrder10=OrderY+1;
SeriesOrder01=OrderX+1;
SeriesOrder20=OrderY+1;
SeriesOrder02=OrderX+1;

for pp=2:NMoments

  SeriesOrder30=SeriesOrder10+SeriesOrder20-1;
  SeriesOrder03=SeriesOrder01+SeriesOrder02-1;
% %  
% %  %Below calculate the product: 2D product of Z-series [Y|X](Zy,Zx)*X(Zx)
  [YpCoeff(1:SeriesOrder30,1:SeriesOrder03,pp),EYp(pp)] =SeriesProduct2DAlgorithmicNew(YpCoeff(1:SeriesOrder10,1:SeriesOrder01,pp-1),YCoeffZZ(1:OrderX+1,1:OrderY+1,1),SeriesOrder10,SeriesOrder01,SeriesOrder20,SeriesOrder02);
    
  YMoments(pp)=EYp(pp);
  
  SeriesOrder10=SeriesOrder30;
SeriesOrder01=SeriesOrder03;
SeriesOrder20=OrderY+1;
SeriesOrder02=OrderX+1;
  
end




end

.
.
.

Amin · Today, 9:21 am

Friends, what I can do is to write a comprehensive function that mimics the above program. It worked extremely well 2-3 days ago but then had problems. And I hope that at the end of the weekend a better, nicer and easy to use version with graphs that would be available that you could also use for eight moments calibration.

Amin · Today, 9:51 am

Friends, I cannot work my programs properly now. Bad Jewish scumbags have taken my programs and will continue to use them and have replaced them with improper code that runs when I use my program. This must stop.

Friends, I tried the above function three days ago and it was working well with gamma=.75 and vol=.65 and zero mean-reversion. For past two days, it continued to work for these values but gave me very skewed graphs and today the iterations have stopped to converge at all (for much less volatility numbers) and I cannot even get skewed graphs that I was getting during past two days.

Amin · Today, 10:05 am

Friends, I really believe that in a civilized society of twenty first century there is absolutely no place for connivance and machinations of a group of rich jews led by jewish billionaire Godfather of mind control who have retarded, at whim, thousands of talented intelligent people in America and abroad that they did not like. These cruel people are inhuman animals even though they have learnt to adjust extremely well within American society.

These corrupt crooks are very powerful and they have already complete control over American army and a large number of politicians. Due to this huge influence and power, it is very difficult to do anything to undo the cruel mind control network led by jewish billionaire and his rich jewish cohort.

I have suggested several times previously to good American people if nothing else works, please just ask some patriotic American investigative journalists in major newspapers to run a comprehensive story about mind control network of jewish billionaire and the nexus between rich jews and American defense and about tens of thousands of people that have been retarded over past thirty years.

I want to invite patriotic American journalists from CNN, Washington Post, NewYork Times, Guardian or some other large, reputed and responsible US or European media outlets to run a thorough story on machinations of evil jews led by jewish billionaire and the nexus of bad jews with American defense and their cruel mind control network that has retarded tens of thousands of very intelligent Americans and foreigners all over the world.

Since good people may not have the power to forcefully end the evil practices of the powerful nexus between evil Jewish billionaire and US army crooks, only way to end the cruel and subhuman practices of the nexus seems to be exposing them openly in US public and the public of other countries all across the world.

I think it is important for us to make sure there is a better future of American children and children all across the world and we leave them a legacy in which fairness and justice prevails as opposed to machinations, trickery and deceit that has been unleashed in the society by the jewish billionaire and his rich jewish cohort.

Amin · 55 minutes ago

A Relevant post copied from before:

Friends, last night after writing my post, I took a lot of protective measures and tried to be serious about breaking their control. And then I went to sleep without doing any further work. Today, I was a bit better in the morning but when I went out to buy food, it was a challenge and I could see at several places that mind control agents reached when I was trying to buy something.

What I am writing is my honest opinion even though a lot of friends might differ with it to some degree. I think it is very important that all of us try to have rule of law, and good governance for the benefit of our country and society. We promote freedom of speech and freedom of press in our societies so that people can criticize bad or wrong actions of other people and make each other aware of their observations of wrongdoings by some people especially the ones in the authority in the society while, obviously, every body has to try to investigate the truth on their own without blindly believing anyone. The purpose of such criticism and related discussions is to try to root out bad practices from the society and promote good ethics and practices.

Another thing I want to suggest is that true patriotism is not about hiding bad practices taking place in our country from others, instead it is about exposing such bad practices by unethical people in the public with the final goal that these bad practices are rooted out from our societies.

Now I will come to the point that I have been suggesting for years that there is a jewish billionaire and his rich cohort who are behind indiscriminate mind control of a large number of intelligent American people and others . And a large number of people in American defense have become a mercenerary army that play in the hands of these rich jews. These rich jews hire quite a few of senior officers from Pentagon and army after their retirement in their large enterprises and these people are paid millions of dollars each and they enjoy very high social status. One of the most important purpose of these retired officers is to convince their juniors who are still in official service and have reached positions of authority to do illegal, extrajudicial and extraconstitutional mind control control of Americans and others suggested by the rich jewish cohort. This thing is so well known in some jewish circles that many jews continuously report thousands of people among intelligent students and many others who fall out with these jews for any reason to mind control by American defense. You do not always have to be extra intelligent, sometimes falling out with a jew who knows how to approach the rich jewish cohort will seal the fate of the victim by mind control forever. As always, I would like to mention that there are large number of jews who are opposed to this practice and have absolutely nothing to do with mind control of innocent people. Please let me specifically say again that there are many very rich jews and billionaires who have nothing to do with mind control. But I definitely own that I am suggesting this rich jewish billionaire and his cohort of rich jews are behind most of the mind control and this is a common knowledge in many jewish circles and many people among these circles know who to approach for mind control of anyone they do not like.

But what is even more alarming is that all the bad actors behind mind control including rich jewish cohort and mercenary criminals in American defense, when exposed in public, have absolutely no remorse about their crimes and they become even more adamant to use all machinatory tactics so that they+ can continue illegal and extrajudicial mind control of innocent and intelligent people. They have lost all sense of good and evil and have no qualms about continuing the evil practices at all. It seems many of these bad actors have no conception that retarding innocent human beings is a wrong thing and an evil practice and they think that mind control of innocent people in the society they do not like is indeed a good thing.

I know that many good Americans want to end the bad practice of mind control of innocent people from their society without any retribution against the bad actors behind such evil practices. The bad actors, on the other hand, have absolutely no remorse about their evil actions and are ready to go to any limit to thwart the attempts of good people to end mind control and they always try to undermine all good people and their actions towards amicably ending the evil of mind control from American society. These bad actors are not willing to gracefully acknowledge and be thankful that good American people do not want to hurt them beyond ending the inhuman practice of mind control from American society and instead these bad actors are adamantly using all sort of machinatory and bellicose tactics so they can continue and increase the evil of inhuman mind control of intelligent and innocent human beings.

I really want to suggest if all amicable attempts by good Americans to end mind control fail, it does not mean that they have to yield to inhuman practices in our society by continuing to play nice. As there is very huge evidence about existence of extra-judicial, extra-constitutional and cruel mind control and resulting retarding of very large number of innocent people by American defense and a large number of stories are circulating in academia about people who were targeted by American defense on behest of jewish cohort, I really believe that it is a great opportunity for some good journalists in American and international media to write in detail about the hugely widespread evil of extrajudicial, extraconstitutional, inhuman and cruel practice of mind control and also about nexus between some bad jews and corrupt people in American defense that has been the dominating force behind this inhuman evil in our societies. I want to request good American people that if all else fails, please motivate some brave and patriotic journalists to write in complete detail about the evil of mind control. If it is not possible to write about this inhuman practice in American newspapers due to possible reprisal from american defense, please break the news in some large British or European newspaper. I am sure there is enough evidence about existence of mind control and related inhuman and cruel practices that any serious journalist can firmly stand behind a well written article that completely exposes mind control and associated cruelties and retarding of thousands of innocent persons by corrupt people in American defense. As I earlier suggested Patriotism does not consist in saving any false pride or hiding evil practices of our countrymen, it consists in rationally doing what is broader good for people in our societies by trying to put an end to evil practices.

I will also request readers, for good of our societies, to let at least one journalist you know in the following American, European and Pakistani newspapers about true nature of mind control and about my thread.

NewYork Times.
Washington Post.
CNN.
Fox News.
USA Today.
Los Angeles Times.
Chicago Tribune.
Boston Globe.
Minnesota Star Tribune
NewYork Post
Newsday
San Fransisco Chronicle
The Dallas Morning News
Tampa Bay Times
Hartford Courant
Publick Occurrences Both Forreign and Domestick
The Michigan Daily
The Denver Post
Pittsburgh Post-Gazette
The Seattle Times
Las Vegas Sun
Bourbon County Citizen
The Register Star
NewYork Metro
The Queens Ledger
Philadelphia Tribune
Baltimore Afro-American
Portland Tribune
Indianapolis Recorder
The Patriot-News
El Diario NY
The Modesto Bee
Fairfax County Times
The Epoch Times
Barron's
Queens Chronicle
Phoenix New Times
El Informador
The Philadelphia Inquirer
NewYork Daily News
The Star-Ledger
St. Louis Post-Dispatch
The Mercury News
The Arizona Republic
The Plain Dealer
Detroit Free Press
The Press-Enterprise
The Tennessean
Sun Sentinel

The Copenhagen Post
The Slovak Spectator
The Budapest Times
Tirana Times
Menorca Sun
The Vienna Review
Le News
The Irish Times
Die Welt
Het Nieuwsblad
Le Soir
Albania Today
Algemeen Dagblad
De Telegraaf
de Volkskrant
The European
Het Laatste Nieuws
The Daily Telegraph
Daily Mirror
The Sun
Pravda
The Moscow Times
The Guardian
Le Figaro
Le Monde
ABC
la Repubblica
Il Messaggero
El Periódico de Catalunya
El Mundo
Libération
El País
Kronen Zeitung
La Stampa
Frankfurter Allgemeine Zeitung
Bild
20 minutes
Sport
Aftenposten
Marca
Gazeta Wyborcza
Izvestia
Die Presse
Süddeutsche Zeitung
La Tribune
Novaya Gazeta
Neue Zürcher Zeitung
Dagens Nyheter
Expressen
La Gazzetta dello Sport
Mladá fronta Dnes

Urdu News
Daily Pakistan
Dawn
Daily Express
Daily Jang
The Express Tribune
The News International
The Nation
Pakistan Today
Daily Dunya
Daily Ummat
Nawaiwaqt
The Regional Times of Sindh
Pakistan Observer
Khabrain
Daily Ausaf
Daily Jasarat
Daily Nai Baat
Daily Basharat
Daily Times
The Frontier Post
Daily Sarhad
Daily Imroze
The Pakistan Times
The Financial Daily
Daily Aaj
Daily Mashriq
Millat
Business Recorder
Daily Inqilab
Kawish
Daily Lokaai
Talár
Maidan Daily
Daily Asas
Christian Voice
The Friday Times
Awam
Daily Awami Awaz
The Star
Al Akhbar
Khyber Mail
Nawai Watan
Sajjan
Daily Hilal Pakistan
Daily Sindhu
Daily Koshish
Daily Mehran
Bhulekha
Civil and Military Gazette

Amin · 40 minutes ago

Friends, this is an old post I had written more than two years ago and I am copying it due to its relevance to giving pain in my heart and other related factors.

"When I started telling people about mind control nobody believed me and only later most of them realized the truth about what I was trying to tell them. I want to tell friends that I would have been killed/shot in NewYork 25 years ago by some planted mugger (with complete planning) or would have been given poison when my talent was discovered and they failed to retard me since it was impossible for many people to stand or accept a Muslim with special neurotransmitters. I remained alive only because a large number of people in the university had already come to know about it and if I were killed, most of them would have known who was behind it. I bitterly wanted to remain in United States and work there but once I realized that most people in government and defense wanted to retard me with an iron fist, I knew I had no future there and I decided to return to Pakistan and, today, I know I made the right decision.

Later I remained alive since my family was assured when they cooperated that I would not be killed or made "special" and again a large number of people in my broader family, army and several influential people in Pakistan got to know about it. But the way people behind my persecution are becoming impatient again, I am seriously afraid they might do what they never did 25 years ago in New York. I am very afraid they might give money to some mugger to shoot me or create some other way to end my life with full cooperation and planning of Pakistan army. I want to tell friends very clearly that if I had premature seemingly accidental death, it would never be accidental, it would be a completely pre-planned murder.

I have copied above post since I really felt during past ten days that mind control agents are becoming extremely impatient and have started to use tactics that they rarely use otherwise. For example, for several days in a row, I would feel pain in my heart and I was sure that mind control agents were doing it. They can in fact do it (have the technology for it) and have tried such things on me on some rare occasions in the past. In fact, I wrote it in one or two posts several years ago when they were trying to give pain in my heart.

So they might as well be saying in a few weeks that Amin was very obese and have died due to heart attack.
Though anybody can have a heart attack anytime, I would be fifty this September."

I was 52 last September of 2025.

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