function [] = CalculateDensitiesAndCorrsTechStocks02()


%Below, read AAPL, AMZN, GOOGL,META,MSFT,NVDA, NASDAQ100 daily open and
%close data files and extract daily open from the files. All files are
%exactly the same length and therefore corresponding array entries are 
%from the same dates. 
T = readtable('C:\Users\Hp\Documents\HistoricalPricesAAPL2022.csv');
[NDays,NVar]=size(T);
ApplPrices(1:NDays)=table2array(T(1:NDays,2));
T = readtable('C:\Users\Hp\Documents\HistoricalPricesAMZN2022.csv');
AmznPrices(1:NDays)=table2array(T(1:NDays,2));
T = readtable('C:\Users\Hp\Documents\HistoricalPricesGOOGL2022.csv');
GooglPrices(1:NDays)=table2array(T(1:NDays,2));
T = readtable('C:\Users\Hp\Documents\HistoricalPricesMETA2022.csv');
MetaPrices(1:NDays)=table2array(T(1:NDays,2));
T = readtable('C:\Users\Hp\Documents\HistoricalPricesMSFT2022.csv');
MsftPrices(1:NDays)=table2array(T(1:NDays,2));
T = readtable('C:\Users\Hp\Documents\HistoricalPricesNVDA2022.csv');
NvdaPrices(1:NDays)=table2array(T(1:NDays,2));
T = readtable('C:\Users\Hp\Documents\HistoricalPricesNASDAQ1002022.csv');
NasdaqPrices(1:NDays)=table2array(T(1:NDays,2));


%Below calculate daily log-returns of each stock and nasdaq 100.
DApplPrices(1:NDays-1)=log(ApplPrices(2:NDays))-log(ApplPrices(1:NDays-1));
DAmznPrices(1:NDays-1)=log(AmznPrices(2:NDays))-log(AmznPrices(1:NDays-1));
DGooglPrices(1:NDays-1)=log(GooglPrices(2:NDays))-log(GooglPrices(1:NDays-1));
DMetaPrices(1:NDays-1)=log(MetaPrices(2:NDays))-log(MetaPrices(1:NDays-1));
DMsftPrices(1:NDays-1)=log(MsftPrices(2:NDays))-log(MsftPrices(1:NDays-1));
DNvdaPrices(1:NDays-1)=log(NvdaPrices(2:NDays))-log(NvdaPrices(1:NDays-1));
DNasdaqPrices(1:NDays-1)=log(NasdaqPrices(2:NDays))-log(NasdaqPrices(1:NDays-1));

%Below calculate first four raw moments of daily log returns. Since I was 
%dealing with daily log-returns, I thought first four moments should be 
%enough and calculating higher moments for density and correlation 
%calculations for this data would be an over-kill.  
NDays0=NDays-1;
rMuAppl=CalculateFourRawMomentsFromData(DApplPrices,NDays0);
rMuAmzn=CalculateFourRawMomentsFromData(DAmznPrices,NDays0);
rMuGoogl=CalculateFourRawMomentsFromData(DGooglPrices,NDays0);
rMuMeta=CalculateFourRawMomentsFromData(DMetaPrices,NDays0);
rMuMsft=CalculateFourRawMomentsFromData(DMsftPrices,NDays0);
rMuNvda=CalculateFourRawMomentsFromData(DNvdaPrices,NDays0);
rMuNasdaq=CalculateFourRawMomentsFromData(DNasdaqPrices,NDays0);

clf;
%Below construct log-return densities from first four moments
[c10,c1] = CalculateZSeriesDensityFromRawMomentsM4(rMuAppl)
str=input('Look at density of Apple daily log returns')

[c20,c2] = CalculateZSeriesDensityFromRawMomentsM4(rMuAmzn)
str=input('Look at density of Amazon daily log returns')

[c30,c3] = CalculateZSeriesDensityFromRawMomentsM4(rMuGoogl)
str=input('Look at density of Google daily log returns')

[c40,c4] = CalculateZSeriesDensityFromRawMomentsM4(rMuMeta)
str=input('Look at density of Meta daily log returns')

[c50,c5] = CalculateZSeriesDensityFromRawMomentsM4(rMuMsft)
str=input('Look at density of Microsoft daily log returns')
[c60,c6] = CalculateZSeriesDensityFromRawMomentsM4(rMuNvda)
str=input('Look at density of Nvidia daily log returns')


[c70,c7] = CalculateZSeriesDensityFromRawMomentsM4(rMuNasdaq)
str=input('Look at density of Nasdaq100 daily log returns')

%Below calculate Z value for each realization of log-returns. This function
%takes an array of log-returns and returns an array of corresponding Z
%vlues in output array.

[ZAppl] = CalculateZgivenXAndZSeriesC3(DApplPrices,c10,c1);

DApplPrices-(c10+c1(1).*ZAppl+c1(2).*ZAppl.^2+c1(3).*ZAppl.^3)
%If the difference above does not go to zero, their will be inaccuracies in
%correlation calculations and following regressions.
str=input('Look at APPL difference between log-returns and their values constructed from Z-series');


[ZAmzn] = CalculateZgivenXAndZSeriesC3(DAmznPrices,c20,c2);

DAmznPrices-(c20+c2(1).*ZAmzn+c2(2).*ZAmzn.^2+c2(3).*ZAmzn.^3)
%If the difference above does not go to zero, their will be inaccuracies in
%correlation calculations and following regressions.
str=input('Look at AMZN difference between log-returns and their values constructed from Z-series');

[ZGoogl] = CalculateZgivenXAndZSeriesC3(DGooglPrices,c30,c3);

DGooglPrices-(c30+c3(1).*ZGoogl+c3(2).*ZGoogl.^2+c3(3).*ZGoogl.^3)
%If the difference above does not go to zero, their will be inaccuracies in
%correlation calculations and following regressions.
str=input('Look at Googl difference between log-returns and their values constructed from Z-series');


[ZMeta] = CalculateZgivenXAndZSeriesC3(DMetaPrices,c40,c4);
DMetaPrices-(c40+c4(1).*ZMeta+c4(2).*ZMeta.^2+c4(3).*ZMeta.^3)
%If the difference above does not go to zero, their will be inaccuracies in
%correlation calculations and following regressions.
str=input('Look at META difference between log-returns and their values constructed from Z-series');





[ZMsft] = CalculateZgivenXAndZSeriesC3(DMsftPrices,c50,c5);
DMsftPrices-(c50+c5(1).*ZMsft+c5(2).*ZMsft.^2+c5(3).*ZMsft.^3)
%If the difference above does not go to zero, their will be inaccuracies in
%correlation calculations and following regressions.
str=input('Look at MSFT difference between log-returns and their values constructed from Z-series');


[ZNvda] = CalculateZgivenXAndZSeriesC3(DNvdaPrices,c60,c6);
DNvdaPrices-(c60+c6(1).*ZNvda+c6(2).*ZNvda.^2+c6(3).*ZNvda.^3)
%If the difference above does not go to zero, their will be inaccuracies in
%correlation calculations and following regressions.
str=input('Look at NVDA difference between log-returns and their values constructed from Z-series');



[ZNasdaq] = CalculateZgivenXAndZSeriesC3(DNasdaqPrices,c70,c7);
DNasdaqPrices-(c70+c7(1).*ZNasdaq+c7(2).*ZNasdaq.^2+c7(3).*ZNasdaq.^3)
%If the difference above does not go to zero, their will be inaccuracies in
%correlation calculations and following regressions.
str=input('Look at NASDAQ100 difference between log-returns and their values constructed from Z-series');


%Below calculate hermite coefficients from Z-series coefficients
[ch10,ch1]=ConvertZCoeffsToHCoeffs(c10,c1,3);
[ch20,ch2]=ConvertZCoeffsToHCoeffs(c20,c2,3);
[ch30,ch3]=ConvertZCoeffsToHCoeffs(c30,c3,3);
[ch40,ch4]=ConvertZCoeffsToHCoeffs(c40,c4,3);
[ch50,ch5]=ConvertZCoeffsToHCoeffs(c50,c5,3);
[ch60,ch6]=ConvertZCoeffsToHCoeffs(c60,c6,3);
[ch70,ch7]=ConvertZCoeffsToHCoeffs(c70,c7,3);

%In what follows, I will calculate correlations of first three hermites
%between each of the six tech stocks and Nasdaq100 index. And finally I
%will draw scatter-plots of joint observations of daily log-returns of each of the
%tech stocks against Nasdaq100 index returns. And then construct the
%conditional regression based graph of each of the tech stock returns VS
%Nasdaq100 returns.

%Below calculate correlations of first three hermites between APPL
%log-returns and Nasdaq100 log rturns.
rho17_1=sum(ZAppl.*ZNasdaq)/sqrt(sum(ZAppl.^2).*sum(ZNasdaq.^2))
rho17_2=sum((ZAppl.^2-1).*(ZNasdaq.^2-1))/sqrt(sum((ZAppl.^2-1).^2).*sum(((ZNasdaq.^2)-1).^2))
rho17_3=sum((ZAppl.^3-3*ZAppl).*(ZNasdaq.^3-3*ZNasdaq))/sqrt(sum((ZAppl.^3-3.*ZAppl).^2).*sum((ZNasdaq.^3-3*ZNasdaq).^2))

%Below calculate correlations of first three hermites between AMZN
%log-returns and Nasdaq100 log rturns.
rho27_1=sum(ZAmzn.*ZNasdaq)/sqrt(sum(ZAmzn.^2).*sum(ZNasdaq.^2))
rho27_2=sum((ZAmzn.^2-1).*(ZNasdaq.^2-1))/sqrt(sum((ZAmzn.^2-1).^2).*sum(((ZNasdaq.^2)-1).^2))
rho27_3=sum((ZAmzn.^3-3*ZAmzn).*(ZNasdaq.^3-3*ZNasdaq))/sqrt(sum((ZAmzn.^3-3.*ZAmzn).^2).*sum((ZNasdaq.^3-3*ZNasdaq).^2))

%Below calculate correlations of first three hermites between Googl
%log-returns and Nasdaq100 log rturns.
rho37_1=sum(ZGoogl.*ZNasdaq)/sqrt(sum(ZGoogl.^2).*sum(ZNasdaq.^2))
rho37_2=sum((ZGoogl.^2-1).*(ZNasdaq.^2-1))/sqrt(sum((ZGoogl.^2-1).^2).*sum(((ZNasdaq.^2)-1).^2))
rho37_3=sum((ZGoogl.^3-3*ZGoogl).*(ZNasdaq.^3-3*ZNasdaq))/sqrt(sum((ZGoogl.^3-3.*ZGoogl).^2).*sum((ZNasdaq.^3-3*ZNasdaq).^2))

%Below calculate correlations of first three hermites between META
%log-returns and Nasdaq100 log rturns.

rho47_1=sum(ZMeta.*ZNasdaq)/sqrt(sum(ZMeta.^2).*sum(ZNasdaq.^2))
rho47_2=sum((ZMeta.^2-1).*(ZNasdaq.^2-1))/sqrt(sum((ZMeta.^2-1).^2).*sum(((ZNasdaq.^2)-1).^2))
rho47_3=sum((ZMeta.^3-3*ZMeta).*(ZNasdaq.^3-3*ZNasdaq))/sqrt(sum((ZMeta.^3-3.*ZMeta).^2).*sum((ZNasdaq.^3-3*ZNasdaq).^2))



%Below calculate correlations of first three hermites between MSFT
%log-returns and Nasdaq100 log rturns.

rho57_1=sum(ZMsft.*ZNasdaq)/sqrt(sum(ZMsft.^2).*sum(ZNasdaq.^2))
rho57_2=sum((ZMsft.^2-1).*(ZNasdaq.^2-1))/sqrt(sum((ZMsft.^2-1).^2).*sum(((ZNasdaq.^2)-1).^2))
rho57_3=sum((ZMsft.^3-3*ZMsft).*(ZNasdaq.^3-3*ZNasdaq))/sqrt(sum((ZMsft.^3-3.*ZMsft).^2).*sum((ZNasdaq.^3-3*ZNasdaq).^2))

%Below calculate correlations of first three hermites between NVDA
%log-returns and Nasdaq100 log rturns.ays0-1).^3-3*ZNasdaq(1:NDays0-1)))/sqrt(sum((ZMsft(2:NDays0).^3-3.*ZMsft(2:NDays0)).^2).*sum((ZNasdaq(1:NDays0-1).^3-3*ZNasdaq(1:NDays0-1)).^2))

rho67_1=sum(ZNvda.*ZNasdaq)/sqrt(sum(ZNvda.^2).*sum(ZNasdaq.^2))
rho67_2=sum((ZNvda.^2-1).*(ZNasdaq.^2-1))/sqrt(sum((ZNvda.^2-1).^2).*sum(((ZNasdaq.^2)-1).^2))
rho67_3=sum((ZNvda.^3-3*ZNvda).*(ZNasdaq.^3-3*ZNasdaq))/sqrt(sum((ZNvda.^3-3.*ZNvda).^2).*sum((ZNasdaq.^3-3*ZNasdaq).^2))


%Below do calculations to do regressions with daily log-returns of each
%tech stock as dependent variables and Nasdaq100 as independent variable.

dNn=.1/2;   % Normal density subdivisions width. would change with number of subdivisions
Nn=45*4;  % No of normal density subdivisions
NnMid=((1+Nn)/2)*dNn;
ZNasdaq0(1:Nn)=(((1:Nn)*dNn)-NnMid);

%Below reconstruct the Z-series representation of returns of each of the
%tech stock as a function of ZNasdaq, i.e the Z that underlies Nasdaq log
%returns. This is regression of each tech stock log-returns on log-returns of Nasdaq 100. 
XAppl(1:Nn)=(ch10)+rho17_1.*(ch1(1))*ZNasdaq0(1:Nn)+rho17_2.*ch1(2).*(ZNasdaq0(1:Nn).^2-1)+rho17_3.*ch1(3).*(ZNasdaq0(1:Nn).^3-3*ZNasdaq0(1:Nn));
XAmzn(1:Nn)=(ch20)+rho27_1.*(ch2(1))*ZNasdaq0(1:Nn)+rho27_2.*ch2(2).*(ZNasdaq0(1:Nn).^2-1)+rho27_3.*ch2(3).*(ZNasdaq0(1:Nn).^3-3*ZNasdaq0(1:Nn));
XGoogl(1:Nn)=(ch30)+rho37_1.*(ch3(1))*ZNasdaq0(1:Nn)+rho37_2.*ch3(2).*(ZNasdaq0(1:Nn).^2-1)+rho37_3.*ch3(3).*(ZNasdaq0(1:Nn).^3-3*ZNasdaq0(1:Nn));
XMeta(1:Nn)=(ch40)+rho47_1.*(ch4(1))*ZNasdaq0(1:Nn)+rho47_2.*ch4(2).*(ZNasdaq0(1:Nn).^2-1)+rho47_3.*ch4(3).*(ZNasdaq0(1:Nn).^3-3*ZNasdaq0(1:Nn));
XMsft(1:Nn)=(ch50)+rho57_1.*(ch5(1))*ZNasdaq0(1:Nn)+rho57_2.*ch5(2).*(ZNasdaq0(1:Nn).^2-1)+rho57_3.*ch5(3).*(ZNasdaq0(1:Nn).^3-3*ZNasdaq0(1:Nn));
XNvda(1:Nn)=(ch60)+rho67_1.*(ch6(1))*ZNasdaq0(1:Nn)+rho67_2.*ch6(2).*(ZNasdaq0(1:Nn).^2-1)+rho67_3.*ch6(3).*(ZNasdaq0(1:Nn).^3-3*ZNasdaq0(1:Nn));
XNasdaq(1:Nn)=ch70+ch7(1)*ZNasdaq0(1:Nn)+ch7(2).*(ZNasdaq0(1:Nn).^2-1)+ch7(3).*(ZNasdaq0(1:Nn).^3-3*ZNasdaq0(1:Nn))



clf
hold on

scatter(DNasdaqPrices(1:NDays0),DApplPrices(1:NDays0))
plot(XNasdaq(1:Nn),XAppl(1:Nn),'g')
title('Regression of daily log returns of Apple on Daily log-returns of NASDAQ 100'); 
str=input('Look at scatterplot and regression between Apple daily log returns and Nasdaq100 daily log-returns');

clf
hold on

scatter(DNasdaqPrices(1:NDays0),DAmznPrices(1:NDays0))
plot(XNasdaq(1:Nn),XAmzn(1:Nn),'g')
title('Regression of daily log returns of Amazon on Daily log-returns of NASDAQ 100');
str=input('Look at scatterplot and regression between Amazon daily log returns and Nasdaq100 daily log-returns');

clf
hold on

scatter(DNasdaqPrices(1:NDays0),DGooglPrices(1:NDays0))
plot(XNasdaq(1:Nn),XGoogl(1:Nn),'g')
title('Regression of daily log returns of Google on Daily log-returns of NASDAQ 100');
str=input('Look at scatterplot and regression between Google daily log returns and Nasdaq100 daily log-returns');

clf
hold on

scatter(DNasdaqPrices(1:NDays0),DMetaPrices(1:NDays0))
plot(XNasdaq(1:Nn),XMeta(1:Nn),'g')
title('Regression of daily log returns of Meta on Daily log-returns of NASDAQ 100');

str=input('Look at scatterplot and regression between Meta daily log returns and Nasdaq100 daily log-returns');

clf
hold on

scatter(DNasdaqPrices(1:NDays0),DMsftPrices(1:NDays0))
plot(XNasdaq(1:Nn),XMsft(1:Nn),'g')
title('Regression of daily log returns of Microsoft on Daily log-returns of NASDAQ 100');
str=input('Look at scatterplot and regression between Microsoft daily log returns and Nasdaq100 daily log-returns');

clf
hold on

scatter(DNasdaqPrices(1:NDays0),DNvdaPrices(1:NDays0))
plot(XNasdaq(1:Nn),XNvda(1:Nn),'g')
title('Regression of daily log returns of Nvidia on Daily log-returns of NASDAQ 100');
str=input('Look at scatterplot and regression between Nvidia daily log returns and Nasdaq100 daily log-returns');

clf
end

function [rMu] = CalculateFourRawMomentsFromData(DPrices,NDays0)

rMu(1)=sum(DPrices)/NDays0;
rMu(2)=sum(DPrices.^2)/NDays0;
rMu(3)=sum(DPrices.^3)/NDays0;
rMu(4)=sum(DPrices.^4)/NDays0;

end

function [c0,c] = CalculateZSeriesDensityFromRawMomentsM4(rMu)

    mOrder=4;
    [Mu1,cMu] = ConvertRawMomentsToCentralMoments(rMu,mOrder);

    sMu(1)=0;
    sMu(2)=1;
    sMu(3)=cMu(3)/cMu(2).^1.5;
    sMu(4)=cMu(4)/cMu(2).^2.0;
    %sMu(5)=cMu(5)/cMu(2).^2.5;
    %sMu(6)=cMu(6)/cMu(2).^3.0;
    %sMu(7)=cMu(7)/cMu(2).^3.5;
    %sMu(8)=cMu(8)/cMu(2).^4.0;
    
    
    
    iter=20;
    bGuess(1:3)=0;
    bGuess(1)=1;
    [c0,c] = PreSmoothingGuessAdvancedFromGuessBestNewC3(sMu,bGuess,iter);
    

    
    SeriesOrder=3;
    NMoments=4;    
    
    [F,dF] = CalculateMomentsAndDerivatives_0(sMu,c0,c,SeriesOrder,SeriesOrder,NMoments);

da(1,1)=c0;
da(2:SeriesOrder+1,1)=c(1:SeriesOrder);

[Moments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);

%Replace with your own more intelligent objective function if you like.
ObjBest=100000*(abs(sMu(1)-Moments(1)))+abs(sMu(2)-Moments(2))+abs((sMu(3)-Moments(3))^(1/1.5))+abs((sMu(4)-Moments(4))^(1/2.0));

b0Best=c0;
bBest(1:SeriesOrder)=c(1:SeriesOrder);

nn=0;
while((nn<20)&&((abs(F(1,1))>.000000000001) || (abs(F(2,1))>.000000000001) || (abs(F(3,1))>.000000000001) || (abs(F(4,1))>.00000000001)  ))

    nn=nn+1;
    %Below Newton matrix equation to improve the Z-series coefficients guess at previous step.
    da=da-dF\F;

    %b0=median
    b0=da(1,1);
    b(1:SeriesOrder)=da(2:SeriesOrder+1,1);

    %[F,dF] = CalculateCumulantsAndDerivativesFromMoments_0(C,b0,b,SeriesOrder,SeriesOrder,NoOfCumulants);
    [F,dF] = CalculateMomentsAndDerivatives_0(sMu,b0,b,SeriesOrder,SeriesOrder,NMoments);
    [IsValidFlag] =1;% CheckIsValidDensity(b0,b);
    [Moments] = CalculateMomentsOfZSeries(b0,b,SeriesOrder,NMoments);

    ObjNew=100000*(abs(sMu(1)-Moments(1)))+abs(sMu(2)-Moments(2))+abs((sMu(3)-Moments(3))^(1/1.5))+abs((sMu(4)-Moments(4))^(1/2.0));
  
    if((ObjBest>ObjNew) &&( IsValidFlag))
      
       ObjBest=ObjNew;
       b0Best=b0;
       bBest(1:SeriesOrder)=b(1:SeriesOrder);
    end

    da(1,1)=b0;
    da(2:SeriesOrder+1,1)=b(1:SeriesOrder);

end
c0=b0Best;%Best;
c(1:SeriesOrder)=bBest;%Best(1:SeriesOrder);
    
    
    
    c0=c0*sqrt(cMu(2));
    c=c*sqrt(cMu(2));
    
    c0=c0+Mu1;

[FittedMoments] = CalculateMomentsOfZSeries(c0,c,SeriesOrder,NMoments);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Please Uncomment this if you want to see a comparison of input moments and
%fitted moments
rMu
FittedMoments
str=input('Look at comparison of fitted moments with input moments');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%COMMENT lines below if you do not want to see the density of new Z-series random variable with
%fitted moments.
%First calculate normal random variable on a grid below
dNn=.1/2;   % Normal density subdivisions width. would change with number of subdivisions
Nn=45*4;  % No of normal density subdivisions
NnMid=((1+Nn)/2)*dNn;
Z(1:Nn)=(((1:Nn)*dNn)-NnMid);
%Z
%str=input('Look at Z');
%Now calculate Y as a function of normal random variable.
Y(1:Nn)=c0;
for nn=1:3
    Y(1:Nn)=Y(1:Nn)+c(nn)*Z(1:Nn).^nn;
end

%Now take change of densities derivative of Y with respect to normal
DfY(1:Nn)=0;
for nn=2:Nn-1
    DfY(nn) = (Y(nn + 1) - Y(nn - 1))/(Z(nn + 1) - Z(nn - 1));
    %Change of variable derivative for densities
end
DfY(Nn)=DfY(Nn-1);
DfY(1)=DfY(2);

%Now calculate the density of Y from density of normal random variable
%using change of probability derivative.
pY(1:Nn)=0;
for nn = 1:Nn    
    pY(nn) = (normpdf(Z(nn),0, 1))/abs(DfY(nn));
end

plot(Y(1:Nn),pY(1:Nn),'b');
title('Density of the Z-Series Variable with Fitted Moments');

%str=input('Look at the density of Z-series random variable')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

end

function [mu1,cMu] = ConvertRawMomentsToCentralMoments(Mu,mOrder)

%cu
%mu
% %mOrder
%str=input('Look at numbers');
mu1=Mu(1);
for mm=2:mOrder
    cMu(mm)=0;
    for jj=0:mm
        if(jj==0)
            cMu(mm)=cMu(mm)+ (-1)^(mm-jj).*factorial(mm)/factorial(jj)/factorial(mm-jj)*1*mu1.^(mm-jj);
        else
            cMu(mm)=cMu(mm)+ (-1)^(mm-jj).*factorial(mm)/factorial(jj)/factorial(mm-jj)*Mu(jj)*mu1.^(mm-jj);
        end
    end
end

function [c0,c] = PreSmoothingGuessAdvancedFromGuessBestNewC3(cmu,cin,iter)


%Mul=1.0;
SeriesOrder=3;
NMoments=4;

 EZ(1)=0;
 EZ(2)=1;
 for nn=3:NMoments*SeriesOrder+SeriesOrder+2
     if rem(nn,2)==1
         EZ(nn)=0;
     else
         EZ(nn)=EZ(nn-2)*(nn-1);
         EZ(nn);
     end
 end


c0=-(cin(2));  %This is for zero mean condition.
c=cin;
MaxIter=10;

for nn=1:iter     %increase the iterations over coefficients if neeeded


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

[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);

tolerance=.0000001*cmu(3);
mm=0;
while ((abs(M3-cmu(3)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(2)=c(2)-(M3-cmu(3))/dc2;

[M3,dc2] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,3,EZ);

end

c0=-(c(2));

[SecondMoment] = CalculateSecondMomentC3(c0,c);

c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);

%[c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,3,2,EZ);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);

tolerance=.0000001*cmu(4);
mm=0;
while ((abs(M4-cmu(4)) >tolerance ) && (mm<MaxIter))
mm=mm+1;
c(3)=c(3)-(M4-cmu(4))/dc3;

[M4,dc3] = CalculateParticularMomentAndDerivativeOfItsCoeff(c0,c,SeriesOrder,4,EZ);

end

[SecondMoment] = CalculateSecondMomentC3(c0,c);

c=c/sqrt(SecondMoment);
c0=c0/sqrt(SecondMoment);


%[c0,c] = IterateSmoothingGuessOverPreviousMoments(cmu,c,4,2,EZ);


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


end


end

function [EnMoment,dEnMomentdCoeff] = CalculateParticularMomentAndDerivativeOfItsCoeff(a0,a,SeriesOrder,nMoment,EZ)


% if(SeriesOrder>=8)
% a0=-(a(2)+3*a(4)+15*a(6)+105*a(8));% ---1
% end
% if(SeriesOrder<8)
% a0=-(a(2)+3*a(4)+15*a(6));% ---1
% end
% 
% EZ(1)=0;
% EZ(2)=1;
% for nn=3:NMoments*SeriesOrder+NZterms+2
%     if rem(nn,2)==1
%         EZ(nn)=0;
%     else
%         EZ(nn)=EZ(nn-2)*(nn-1);
%         EZ(nn);
%     end
% end
% EZ;

%EXZ(1,1)=1;
%for pp1=1:NZterms
%        EXZ(1,pp1+1)=EZ(pp1);
%end


a(SeriesOrder+1:nMoment*SeriesOrder+1+SeriesOrder)=0;
b0=a0;
b=a;

for mm=2:nMoment
    
    [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*mm+SeriesOrder+1);
     b(SeriesOrder*mm+1:nMoment*SeriesOrder+SeriesOrder+1)=0;
     if(mm==nMoment-1)
         dEnMomentdCoeff=b0.*EZ(nMoment-1);
         for pp2=1:SeriesOrder*nMoment-1
            
            dEnMomentdCoeff=dEnMomentdCoeff+b(pp2).*EZ(pp2+nMoment-1);
        end
        dEnMomentdCoeff=dEnMomentdCoeff*nMoment; 
     end
     if(mm==nMoment)
        EnMoment=b0;
        for pp2=1:SeriesOrder*nMoment
            EnMoment=EnMoment+b(pp2).*EZ(pp2);
    
        end
     end
end   

    

end

function [c0,c] =SeriesProduct(a0,a,b0,b,SeriesOrder)


c0=a0*b0;
c(1:SeriesOrder)=0.0;
for nn=1:SeriesOrder
    c(nn)=c(nn)+a0*b(nn);
    c(nn)=c(nn)+b0*a(nn);
    for kk=1:nn-1
        c(nn)=c(nn)+a(kk)*b(nn-kk);
    end
end

end

function [SecondMoment] = CalculateSecondMomentC3(a0,a)


SecondMoment=a0^2+a(1)^2 +2* a0* a(2) +3*(a(2).^2+2* a(1).* a(3))+ ...
    15*(a(3).^2) ;

end

function [F,dF] = CalculateMomentsAndDerivatives_0(Ms,a0,a,SeriesOrder,NZterms,NMoments)
%[F,dF] = CalculateCumulantsAndDerivativesFromMoments(C,a0,a,SeriesOrder,SeriesOrder,NoOfCumulants);

%a(1)=a1;
%a(2)=a2;
%a(3)=a3;
%a(4)=a4;
%a(5)=a5;
%a(6)=a6;
%a(7)=a7;



%aa0=a0;
if(NMoments>8)
a0=-(a(2)+3*a(4)+15*a(6)+105*a(8));% ---1
end
if(NMoments==8)
a0=-(a(2)+3*a(4)+15*a(6));% ---1
end
if(NMoments==6)
a0=-(a(2)+3*a(4));% ---1
end

EZ(1)=0;
EZ(2)=1;
for nn=3:NMoments*SeriesOrder+NZterms+2
    if rem(nn,2)==1
        EZ(nn)=0;
    else
        EZ(nn)=EZ(nn-2)*(nn-1);
        EZ(nn);
    end
end
EZ;

EXZ(1,1)=1;
for pp1=1:NZterms
        EXZ(1,pp1+1)=EZ(pp1);
end


a(SeriesOrder+1:NMoments*SeriesOrder+1)=0;
b0=a0;
b=a;

for mm=1:NMoments
    if(mm>1)
        [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*mm);
        b(SeriesOrder*mm+1:NMoments*SeriesOrder+1)=0;
    end
   % b0
   % b
%str=input('Look at numbers')    
    EXZ(mm+1,1)=b0;
    for pp2=1:SeriesOrder*mm
        
            EXZ(mm+1,1)=EXZ(mm+1,1)+b(pp2).*EZ(pp2);
    end
    for pp1=1:NZterms
        EXZ(mm+1,pp1+1)=b0.*EZ(pp1);
        for pp2=1:SeriesOrder*mm
        
            EXZ(mm+1,pp1+1)=EXZ(mm+1,pp1+1)+b(pp2).*EZ(pp2+pp1);
        end
    end
end

%u1=EXZ(2,1);
u2=EXZ(3,1);
u3=EXZ(4,1);
u4=EXZ(5,1);

u1=a0+a(2);
if(SeriesOrder>=4)
u1=a0+a(2)+3*a(4);
end
if(SeriesOrder>=6)
u1=u1+15*a(6);
end
if(SeriesOrder>=8)
u1=u1+105*a(8);
end


%k2=u2-u1^2;
%k3=u3-3*u2*u1+2*u1^3;
%k4=u4-4*u3*u1-3*u2^2+12*u2*u1^2-6*u1^4;


du1(1)=1;%----2
% du1(2)=0;%----2
% du1(3)=1;%----2
% du1(4)=0;%----2
% du1(5)=1;%----2
% du1(6)=0;%----2
% du1(7)=15;%----2
% du1(8)=0;%----2
% du1(9)=105;%----2
% du1(10)=0;%----2
 du2(1)=2*EXZ(2,1);
 du3(1)=3*EXZ(3,1);
 du4(1)=4*EXZ(4,1);


%du1(1)=1;%----2
%du2(1)=0;
%du3(1)=0;
%du4(1)=0;

for mm=2:SeriesOrder+1
du1(mm)=EXZ(1,mm);
du2(mm)=2*EXZ(2,mm);
du3(mm)=3*EXZ(3,mm);
du4(mm)=4*EXZ(4,mm);

end


if(NMoments>=5)
u5=EXZ(6,1);


%du5(1)=5*EXZ(5,1);

for mm=1:SeriesOrder+1
du5(mm)=5*EXZ(5,mm);
end

end

if(NMoments>=6)
u6=EXZ(7,1);
%du6(1)=0;

for mm=1:SeriesOrder+1
du6(mm)=6*EXZ(6,mm);
end

end


if(NMoments>=7)
u7=EXZ(8,1);
%str=input('Look at k7 and k71')

%du7(1)=0;

for mm=1:SeriesOrder+1
du7(mm)=7*EXZ(7,mm);
end

end
if(NMoments>=8)

u8=EXZ(9,1);
for mm=1:SeriesOrder+1
du8(mm)=8*EXZ(8,mm);
end

end

if(NMoments>=9)

u9=EXZ(10,1);
for mm=1:SeriesOrder+1
du9(mm)=9*EXZ(9,mm);
end

end
if(NMoments>=10)

u10=EXZ(11,1);
for mm=1:SeriesOrder+1
du10(mm)=10*EXZ(10,mm);
end

end




 F(1,1)=u1-Ms(1);
 F(2,1)=u2-Ms(2);
 F(3,1)=u3-Ms(3);
 F(4,1)=u4-Ms(4);
 if(NMoments>=5)
 F(5,1)=u5-Ms(5);
 end
 if(NMoments>=6)
 F(6,1)=u6-Ms(6);
 end
 if(NMoments>=7)
 F(7,1)=u7-Ms(7);
 end
 if(NMoments>=8)
 F(8,1)=u8-Ms(8);
 end
 if(NMoments>=9)
 F(9,1)=u9-Ms(9);
 end
 if(NMoments>=10)
 F(10,1)=u10-Ms(10);
 end
 
     for mm=1:SeriesOrder+1
        dF(1,mm)=du1(mm);
        dF(2,mm)=du2(mm);
        dF(3,mm)=du3(mm);
        dF(4,mm)=du4(mm);
        if(NMoments>=5)
        dF(5,mm)=du5(mm);
        end
        if(NMoments>=6)
        dF(6,mm)=du6(mm);
        end
        if(NMoments>=7)
        dF(7,mm)=du7(mm);
        end
        if(NMoments>=8)
        dF(8,mm)=du8(mm);
        end
        if(NMoments>=9)
        dF(9,mm)=du9(mm);
        end
        if(NMoments>=10)
        dF(10,mm)=du10(mm);
        end
     end

end

function [Moments] = CalculateMomentsOfZSeries(a0,a,SeriesOrder,NMoments)


%aa0=a0;
%a0=0;% ---1

EZ(1)=0;
EZ(2)=1;
%LogEZ(2)=0;
for nn=3:NMoments*SeriesOrder
    if rem(nn,2)==1
        EZ(nn)=0;
    else
        EZ(nn)=EZ(nn-2)*(nn-1);
        %LogEZ(nn)=log(EZ(nn-2))+log(nn-1);
        %EZ(nn);
    end
end
%EZ
%EZ(1:30)
%str=input('Look at numbers')
a(SeriesOrder+1:SeriesOrder*NMoments+1)=0;
b0=a0;
b=a;


for mm=1:NMoments
    if(mm>1)
        %[b0,b] =SeriesProductLogarithmic(a0,a,b0,b,SeriesOrder*mm);
        
        [b0,b] =SeriesProduct(a0,a,b0,b,SeriesOrder*mm);
        
        %b0-bb0
        
        %b-bb
        
        %str=input('Look at two products')
        
        
        b(SeriesOrder*mm+1:SeriesOrder*NMoments+1)=0;
    end
   % Logb=log(abs(b));
   % Signb=sign(b);
    EXZ(mm,1)=b0;
    for pp2=1:SeriesOrder*mm
        if rem(pp2,2)==0
            %EXZ(mm,1)=EXZ(mm,1)+exp(Logb(pp2)+LogEZ(pp2)).*Signb(pp2);
            EXZ(mm,1)=EXZ(mm,1)+b(pp2).*EZ(pp2);
 %           b(pp2).*EZ(pp2)
 %           exp(Logb(pp2)+LogEZ(pp2)).*Signb(pp2)
 %           mm
 %           pp2
            %str=input('Look at moment values')
        end
    end
end

    
for nn=1:NMoments
    Moments(nn)=EXZ(nn,1);
end



end

function [Z] = CalculateZgivenXAndZSeriesC3(X,c0,c)


Z=(X-c0)/c(1)-(c(2) * (X-c0).^2)/c(1).^3+(((2* c(2).^2)./c(1).^2-c(3)./c(1)).* (X-c0).^3)./c(1).^3;

for nn=1:30

    X0=c0+c(1)*Z+c(2)*Z.^2+c(3).*Z.^3;
    dX0dZ=c(1)+2.*c(2)*Z+3.*c(3).*Z.^2;
    Z=Z-(X0-X)./(dX0dZ);
end



end

function [aH0,aH] = ConvertZCoeffsToHCoeffs(a0,a,SeriesOrder)

if(SeriesOrder==3)
aH0=a0+a(2);
aH(1)=a(1)+3*a(3);
aH(2)=a(2);
aH(3)=a(3);
end


if(SeriesOrder==4)
aH0=a0+a(2)+3*a(4);
aH(1)=a(1)+3*a(3);
aH(2)=a(2)+6*a(4);
aH(3)=a(3);
aH(4)=a(4);
end


if(SeriesOrder==5)
aH0=a0+a(2)+3*a(4);
aH(1)=a(1)+3*a(3)+15*a(5);
aH(2)=a(2)+6*a(4);
aH(3)=a(3)+10*a(5);
aH(4)=a(4);
aH(5)=a(5);
end



if(SeriesOrder==6)
aH0=a0+a(2)+3*a(4)+15*a(6);
aH(1)=a(1)+3*a(3)+15*a(5);
aH(2)=a(2)+6*a(4)+45*a(6);
aH(3)=a(3)+10*a(5);
aH(4)=a(4)+15*a(6);
aH(5)=a(5);
aH(6)=a(6);
end


if(SeriesOrder==7)
aH0=a0+a(2)+3*a(4)+15*a(6);
aH(1)=a(1)+3*a(3)+15*a(5)+105*a(7);
aH(2)=a(2)+6*a(4)+45*a(6);
aH(3)=a(3)+10*a(5)+105*a(7);
aH(4)=a(4)+15*a(6);
aH(5)=a(5)+21*a(7);
aH(6)=a(6);
aH(7)=a(7);
end



end