   % This code computes additional statistics in the
   % Code_earnings_nonparametric_bootstrap_parallel.m file provided in
   % Arellano, Blundell and Bonhomme (2017)

    % 1. Autocorrelation of persistent component
    % Figure 3c
    
    Mtt1=Mateta_trueB(:,1);
    Mtt2=Mateta_trueB(:,2);
    Mtt3=Mateta_trueB(:,3);
    rhoeta_boot=zeros(36,1);
    
    for i=25:60
        Perseta_temp=Mtt1(AGEB(:,1)==i);
        Perseta_temp2=Mtt2(AGEB(:,1)==i);
        Perseta_temp3=Mtt2(AGEB(:,2)==i);
        Perseta_temp4=Mtt3(AGEB(:,2)==i);
        Perseta_temp5=[Perseta_temp; Perseta_temp3];
        Perseta_temp6=[Perseta_temp2; Perseta_temp4];
        temp=(cov(Perseta_temp5,Perseta_temp6)./...
            var(Perseta_temp5));
        rhoeta_boot(i-24)=temp(1,2);
    end
    
    Resboot_pers_eta(:,jboot)=rhoeta_boot;
    
    % 1. Stdv of shocks to persistent component
    % Fig 3a
    % Also skewness and kurtosis, 4c and 4d
    
    Mateta_long2=reshape(Mateta_trueB(:,2:3),NB*(T-1),1);
    Mateta_long1=reshape(Mateta_trueB(:,1:2),NB*(T-1),1);
    Ageeta_trueB=reshape(AGEB(:,1:2),NB*(T-1),1);
    Shocketa_trueB=Mateta_long2-...
        rhoeta_boot(max(min(Ageeta_trueB,60),25)-24).*Mateta_long1; 
    Mateta_shock_trueB=reshape(Mateta_trueB(:,1:2),NB*(T-1),1);
    
    sd_eta_boot=zeros(36,1);
    skew_eta_boot=sd_eta_boot;
    kurt_eta_boot=sd_eta_boot;
    
    for i=25:60
        Shocketa_temp=Shocketa_trueB(Ageeta_trueB==i);
        sd_eta_boot(i-24)=std(Shocketa_temp);
        skew_eta_boot(i-24)=mean(((Shocketa_temp-mean(Shocketa_temp))./ ...
            sd_eta_boot(i-24)).^3);
        kurt_eta_boot(i-24)=mean(((Shocketa_temp-mean(Shocketa_temp))./ ...
            sd_eta_boot(i-24)).^4);
    end
    
    Resboot_sdshock_eta(:,jboot)=sd_eta_boot;
    Resboot_skew2_eta(:,jboot)=skew_eta_boot;
    Resboot_kurt_eta(:,jboot)=kurt_eta_boot;
    
    
    % 2. Std of transitory component (Fig 3b)
    % Also skewness (4a) and kurtosis (4b)
    
    
    Mateps_long=reshape(Mateps_trueB(:,1:3),NB*(T),1);
    sd_eps_boot=zeros(36,1);
    skew_eps_boot=sd_eps_boot;
    kurt_eps_boot=sd_eps_boot;
    
    for i=25:60
        Shockeps_temp=Mateps_long(reshape(AGEB,NB*T,1)==i);
        sd_eps_boot(i-24)=std(Shockeps_temp);
        skew_eps_boot(i-24)=mean(((Shockeps_temp-mean(Shockeps_temp))./ ...
            sd_eps_boot(i-24)).^3);
        kurt_eps_boot(i-24)=mean(((Shockeps_temp-mean(Shockeps_temp))./ ...
            sd_eps_boot(i-24)).^4);
    end
    
    Resboot_sdshock_eps(:,jboot)=sd_eps_boot;
    Resboot_skew_eps(:,jboot)=skew_eps_boot;
    Resboot_kurt_eps(:,jboot)=kurt_eps_boot;
    

    
    sdpc_boot=zeros(100,1);
    
    % 4. SD by quantile of previous persistent component
    % Fig 5a
    
    xpre=-10000;
    for i=1:100
        x=quantile(Mateta_shock_trueB,i/100);
        Shocketa_temp=Shocketa_trueB(Mateta_shock_trueB<x & ...
            Mateta_shock_trueB>xpre);
        sdpc_boot(i)=std(Shocketa_temp);
        xpre=x;
    end
    
    Resboot_sdpc_eta(:,jboot)=sdpc_boot;
    