%% Check moments in process
clear

U=csvread('../Processes/eta_space.csv',1,0);
U2=csvread('../Processes/eta_tranmatrices.csv',1,0);
U3=csvread('../Processes/eps_space.csv',1,0);
dy=18;TL=36;dye=8;

U4=csvread('../Processes/eta_space_36.csv',1,0);
U5=csvread('../Processes/eta_tranmatrices_36.csv',1,0);
U6=csvread('../Processes/eps_space_16.csv',1,0);
dy2=36;dye2=16;

invher=[0.02*ones(5,1); 0.1*ones(8,1); 0.02*ones(5,1)];
invher2=[0.0005;0.0005;0.002;0.002;0.005;0.01;0.01;0.01;0.01;...
        0.05*ones(18,1);0.01;0.01;0.01;0.01;0.005;0.002;0.002;0.0005;0.0005];


ylifelog=reshape(U(:,3),dy,TL);
Qylife=permute(reshape(U2(:,4),dy,dy,TL-1),[2 1 3]);
    
ye=exp(reshape(U3(:,3),dye,TL));
yedtb=repmat([0.025;0.025;0.05;0.4;0.4;0.05;0.025;0.025],1,TL);

ylifelog2=reshape(U4(:,3),dy2,TL);
Qylife2=permute(reshape(U5(:,4),dy2,dy2,TL-1),[2 1 3]);

ye2=exp(reshape(U6(:,3),dye2,TL));
yedtb2=repmat([0.001;0.004;0.005;0.02;0.02;0.05;0.2;0.2;0.2;0.2;0.05;0.02;0.02;...
        0.005;0.004;0.001],1,TL);
    
% My interest is in ylifelog and Qylife

sdlife=zeros(dy,TL-1);
sklife=zeros(dy,TL-1);
ktlife=zeros(dy,TL-1);
cslife=zeros(dy,TL-1);
kelife=cslife;

for t=1:TL-1
    for i=1:dy
    yhoy=kron(ones(dy,1),ylifelog(i,t));
    yman=ylifelog(:,t+1);
    ng=yman-yhoy;
    probs=reshape(Qylife(i,:,t),dy,1);
    
    [ng,I]=sort(ng);
    probs=probs(I);
    cumprobs=zeros(length(probs));
    cumprobs(1)=probs(1)/2;
    for j=2:length(probs)
        cumprobs(j)=cumprobs(j-1)+probs(j-1)/2+probs(j)/2;
    end
    
    J=(probs<eps);
    cumprobs=cumprobs(~J);
    ng=ng(~J);
    probs=probs(~J);
    
    
    cspoints=[0.025 0.25 0.75 0.975 0.1 0.5 0.9];
    if length(ng)>1
    cspieces=interp1(cumprobs,ng,cspoints);
    cspieces(isnan(cspieces))=interp1(cumprobs,ng,cspoints(isnan(cspieces)),...
       'nearest','extrap');
    end
    
    kelife(i,t)=((cspieces(7)-cspieces(6))-(cspieces(6)-cspieces(5)))./(cspieces(7)-cspieces(5));
    cslife(i,t)=(cspieces(4)-cspieces(1))/(cspieces(3)-cspieces(2));
    
    condmean=sum(ng.*probs);
    sdlife(i,t)=sqrt(var(ng,probs));
    
    sklife(i,t)=sum((((ng-condmean)/sdlife(i,t)).^3).*probs);
    ktlife(i,t)=sum((((ng-condmean)/sdlife(i,t)).^4).*probs);
    
    
    end
end

csgraph=zeros(dy,3);
kegraph=zeros(dy,3);
sdgraph=csgraph;
skgraph=csgraph;
ktgraph=csgraph;
for i=1:3
   csgraph(:,i)=mean(cslife(:,(i-1)*10+1:10*i),2);
   kegraph(:,i)=mean(kelife(:,(i-1)*10+1:10*i),2);
   sdgraph(:,i)=mean(sdlife(:,(i-1)*10+1:10*i),2);
   skgraph(:,i)=mean(sklife(:,(i-1)*10+1:10*i),2);
   ktgraph(:,i)=mean(ktlife(:,(i-1)*10+1:10*i),2);
end
csvwrite('Outputs/mm_sd.csv',sdgraph)
csvwrite('Outputs/mm_sk.csv',skgraph)
csvwrite('Outputs/mm_kt.csv',ktgraph)
csvwrite('Outputs/mm_ke.csv',kegraph)
csvwrite('Outputs/mm_cs.csv',csgraph)

% And unconditional moments of income

sdunc=zeros(TL,1);
ktunc=zeros(TL,1);
skunc=zeros(TL,1);

for i=1:TL
   condmean=sum(ylifelog(:,i).*invher);
   sdunc(i)=sqrt(var(ylifelog(:,i),invher));
   skunc(i)=sum((((ylifelog(:,i)-condmean)/sdunc(i)).^3).*invher);
   ktunc(i)=sum((((ylifelog(:,i)-condmean)/sdunc(i)).^4).*invher);
end


csvwrite('Outputs/mm_sdunc.csv',sdunc)
csvwrite('Outputs/mm_skunc.csv',skunc)
csvwrite('Outputs/mm_ktunc.csv',ktunc)

% And also for the transitory component...

sdtrans=zeros(TL,1);
sktrans=zeros(TL,1);
kttrans=zeros(TL,1);

for i=1:TL
    cmean=sum(log(ye(:,i)).*yedtb(:,i));
    sdtrans(i)=sqrt(var(log(ye(:,i)),yedtb(:,i)));
    sktrans(i)=sum((((log(ye(:,i))-cmean)/sdtrans(i)).^3).*yedtb(:,i));
    kttrans(i)=sum((((log(ye(:,i))-cmean)/sdtrans(i)).^4).*yedtb(:,i));
end

csvwrite('Outputs/mm_sdtrans.csv',sdtrans);
csvwrite('Outputs/mm_sktrans.csv',sktrans);
csvwrite('Outputs/mm_kttrans.csv',kttrans);

% And unconditional moments of income (extra gridpoints)

sdunc2=zeros(TL,1);
ktunc2=zeros(TL,1);
skunc2=zeros(TL,1);

for i=1:TL
   condmean2=sum(ylifelog2(:,i).*invher2);
   sdunc2(i)=sqrt(var(ylifelog2(:,i),invher2));
   skunc2(i)=sum((((ylifelog2(:,i)-condmean2)/sdunc2(i)).^3).*invher2);
   ktunc2(i)=sum((((ylifelog2(:,i)-condmean2)/sdunc2(i)).^4).*invher2);
end


csvwrite('Outputs/mm_sdunc2.csv',sdunc2)
csvwrite('Outputs/mm_skunc2.csv',skunc2)
csvwrite('Outputs/mm_ktunc2.csv',ktunc2)

% And also for the transitory component... (extra gridpoints)

sdtrans2=zeros(TL,1);
sktrans2=zeros(TL,1);
kttrans2=zeros(TL,1);

for i=1:TL
    cmean2=sum(log(ye2(:,i)).*yedtb2(:,i));
    sdtrans2(i)=sqrt(var(log(ye2(:,i)),yedtb2(:,i)));
    sktrans2(i)=sum((((log(ye2(:,i))-cmean2)/sdtrans2(i)).^3).*yedtb2(:,i));
    kttrans2(i)=sum((((log(ye2(:,i))-cmean2)/sdtrans2(i)).^4).*yedtb2(:,i));
end

csvwrite('Outputs/mm_sdtrans2.csv',sdtrans2);
csvwrite('Outputs/mm_sktrans2.csv',sktrans2);
csvwrite('Outputs/mm_kttrans2.csv',kttrans2);
