%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute counterfactual histories for MSAs of choice
%
% Also output shock history for use in Stata
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


TargetMSA = { ...
    37980; ... Philly
    };

close all;

%%
% Compute / back out coefficient matrix
% later, could do for all coefficients

OLD_A_IV=A_IV;
OLD_Ap_IV=Ap_IV;

Mean_A_IV = cell(NoGroups,1);
Mean_Ap_IV = cell(NoGroups,pp);
% Mean_RHOs = ;

OLD_DATA=DATA;
B_REPS_current=B_REPS;

for b_idx=1:B_REPS
    for gg=1:NoGroups
    %     Mean_K{gg} = mean(squeeze(BootstrapData.K_IV_combined(:,:,gg,:)),3);
        Mean_A_IV{gg} = squeeze(BootstrapData.A_IV(:,:,gg,b_idx));
    %     Mean_Ap_IV{gg} = mean(squeeze(BootstrapData.Ap_IV(:,:,gg,:)),3);
        for p_idx = 1:Params.pp
            Mean_Ap_IV{gg,p_idx } = squeeze(BootstrapData.Ap_IV(:,p_idx ,gg,b_idx));
        end

        Mean_beta{gg} = squeeze(IRF_pooled_b{1,gg}(:,:,1, b_idx));
        Mean_V{gg} = squeeze(BootstrapData.V_opt(:,:,gg,b_idx));

        chol_Mean_V{gg} = chol(Mean_V{gg}(1:p,1:p), 'lower');
    %     chol_Mean_V{gg}*chol_Mean_V{gg}'-Mean_V{gg};
        Q_IV = inv(chol_Mean_V{gg})*Mean_beta{gg}(1:p,:);
        Q_orth = null(Q_IV');
        Q = [Q_IV, Q_orth];
        Q_IV'*Q_orth;
        beta{gg} = chol_Mean_V{gg}*Q;
    %     Q' * inv(chol_Mean_V{gg}) * beta;
        if b_idx==1
            Q_IV*Q_IV'
            Q'*Q
            beta{gg} * beta{gg}' - Mean_V{gg}
        end

    end
        Mean_RHOs = reshape( BootstrapData.RHOS_opt(:,:,b_idx), (p+pp+pz),NoGroups);
        % (4) Compute residuals 

        CF_DATA=OLD_DATA;

        % (4a) Residuals net of trend or year fixed effects
        Res_Lvl=nan(size(Y));
        pRes_Lvl=nan(size(pY));
        if Params.DO_YEAR_FE
            for gg=1:NoGroups
                Res_Lvl(:,:,1==IndexM(:,gg))=reshape( ...
                    reshape(Y(:,:,1==IndexM(:,gg)),Params.p,(Params.T-Params.MissingObs-Params.k)*sum(IndexM(:,gg)) ) ...
                            -Mean_A_IV{gg}'*reshape(X(1:Params.p*Params.k,:,1==IndexM(:,gg)), ...
                    Params.p*Params.k,(Params.T-Params.MissingObs-Params.k)* sum(IndexM(:,gg))), [Params.p, Params.T-Params.MissingObs-Params.k, sum(IndexM(:,gg)) ]);
                for p_idx = 1:Params.pp
                    pRes_Lvl(p_idx,:,1==IndexM(:,gg))=reshape( ...
                        reshape(pY(p_idx,:,1==IndexM(:,gg)),1,(Params.T-Params.MissingObs-Params.k)*sum(IndexM(:,gg)) ) ...
                                    -Mean_Ap_IV{gg,p_idx}'*reshape(pX(:,:,1==IndexM(:,gg),p_idx), Params.ppx,(Params.T-Params.MissingObs-Params.k)* sum(IndexM(:,gg))), ...
                                [1, Params.T-Params.MissingObs-Params.k, sum(IndexM(:,gg)) ]);
                end
            end
        %     EXO_YEAR=ones(Params.N,1);
        %     RES_MAKER=eye(Params.N)-EXO_YEAR*( (EXO_YEAR'*EXO_YEAR)\EXO_YEAR' );
        %     Normalize time 1 year fixed effect to zero?
            CF_DATA.Year_FE_res=zeros(Params.p, Params.T-Params.MissingObs-Params.k, NoGroups);
            CF_DATA.Year_FE_pres=zeros(Params.pp, Params.T-Params.MissingObs-Params.k, NoGroups);
            CF_DATA.Year_FE_Z=zeros(Params.pz, Params.T-Params.MissingObs-Params.k, NoGroups);
            for t_idx=1:Params.T-Params.k-Params.MissingObs
                for gg=1:NoGroups
                    CF_DATA.Year_FE_res(:,t_idx, gg) = mean(squeeze(Res_Lvl(:,t_idx,1==IndexM(:,gg))),2);
                    CF_DATA.Year_FE_pres(:,t_idx, gg) = mean(reshape(pRes_Lvl(:,t_idx,1==IndexM(:,gg)), [Params.pp, 1, sum(IndexM(:,gg))]),3);
                    Res_Lvl(:,t_idx,1==IndexM(:,gg))=Res_Lvl(:,t_idx,1==IndexM(:,gg)) - repmat( CF_DATA.Year_FE_res(:,t_idx,gg) , [1, 1, sum(IndexM(:,gg))]);
                    pRes_Lvl(:,t_idx,1==IndexM(:,gg))=pRes_Lvl(:,t_idx,1==IndexM(:,gg)) - repmat( CF_DATA.Year_FE_pres(:,t_idx,gg) , [1, 1, sum(IndexM(:,gg))]);
            %         CF_DATA.Year_FE_Z(:,t_idx, :) = mean(squeeze(Z(:,t_idx,:)),2);
                    CF_DATA.Year_FE_Z(:,t_idx, gg) = mean(reshape(Z(:,t_idx,1==IndexM(:,gg)), [Params.pz, sum(IndexM(:,gg))]),2);
                    Z(:,t_idx,1==IndexM(:,gg))=Z(:,t_idx,1==IndexM(:,gg)) - repmat( CF_DATA.Year_FE_Z(:,t_idx, gg) , [1, 1, sum(IndexM(:,gg))]);
                end
            end
        else
            for gg=1:NoGroups
                Res_Lvl(:,:,1==IndexM(:,gg))=reshape( reshape(Y(:,:,1==IndexM(:,gg)),Params.p,(Params.T-Params.MissingObs-Params.k)*sum(IndexM(:,gg)))- ...
                    [A_IV{gg}', dY_Mean{gg}]*reshape(X(1:Params.px,:,1==IndexM(:,gg)),Params.px,(Params.T-Params.MissingObs-Params.k)*sum(IndexM(:,gg))), [Params.p, Params.T-Params.MissingObs-Params.k, sum(IndexM(:,gg))]);
            end
        end

        % (4b) Residuals net of MSA fixed effects
        EXO=ones(Params.T-Params.k-Params.MissingObs,1);
        PROJECTION=EXO*( (EXO'*EXO)\EXO' );
        CF_DATA.MSA_FE_res=zeros(Params.p, Params.N);
        CF_DATA.MSA_FE_pres=zeros(Params.pp, Params.N);
        CF_DATA.MSA_FE_Z=zeros(Params.pz, Params.N);
        for n_idx=1:Params.N
            CF_DATA.MSA_FE_res(:,n_idx)=mean(squeeze(Res_Lvl(:,:,n_idx)),2);
            Res_Lvl(:,:,n_idx)=Res_Lvl(:,:,n_idx)-repmat(reshape(CF_DATA.MSA_FE_res(:,n_idx), [Params.p, 1,1]),[1,Params.T-Params.k-Params.MissingObs,1]);

            CF_DATA.MSA_FE_Z(:,n_idx)=mean(squeeze(Z(:,:,n_idx)),2);
            Z(:,:,n_idx)=Z(:,:,n_idx) - repmat(reshape(CF_DATA.MSA_FE_Z(:,n_idx), [Params.pz, 1,1]),[1,Params.T-Params.k-Params.MissingObs,1]);

            CF_DATA.MSA_FE_pres(:,n_idx)=mean(squeeze(pRes_Lvl(:,:,n_idx)),2);
            pRes_Lvl(:,:,n_idx)=pRes_Lvl(:,:,n_idx)-repmat(reshape(CF_DATA.MSA_FE_pres(:,n_idx), [Params.pp, 1,1]),[1,Params.T-Params.k-Params.MissingObs,1]);

        end

        [~,~, CF_E_opt, ~, CF_pE_opt]=my_spatial_LL_split_Periphery_jae(Res_Lvl, pRes_Lvl, Z, Params.T-Params.MissingObs-Params.k, Params.N, DistanceMatrix1, Mean_RHOs,[], IndexM);

    %% Counterfactual: without any start-up shock
    % 1) Re-produce bootstrap code: Given E_opt everywhere, make sure I obtain
    % the original time series
    % 2) Transform forecast errors into structural shocks everywhere
    % 3) 1st decomposition: Philly without identified shocks everywhere
    % 4) 2nd decomposition: Philly without identified shocks in Philly, but
    % with everywhere

%     OLD_DATA=DATA;
    DATA=CF_DATA;

    A_IV=Mean_A_IV ;
    Ap_IV=Mean_Ap_IV ;

%     B_REPS_current=1;

    GroupedRHO_opt =Mean_RHOs(1:(p+pz+pp), 1:NoGroups );

    if b_idx==1
        
        CF_BootstrapData=matfile(['CF_BootstrapD' Suffix '.mat'], 'writable', true);

        CF_BootstrapData.K_IV_combined=zeros([size(K_IV_combined{1}),NoGroups, B_REPS_current]);
        CF_BootstrapData.A_IV=zeros(p*k,p,NoGroups,B_REPS_current);
        CF_BootstrapData.RHOS_opt=zeros(p+pz+pp,NoGroups,B_REPS_current);
        CF_BootstrapData.RHOS_cons=zeros(1,NoGroups,B_REPS_current);
        CF_BootstrapData.V_opt=zeros(p,p,NoGroups,B_REPS_current);
        CF_BootstrapData.Ap_IV=zeros(k+(1+k)*p,pp,NoGroups,B_REPS_current);
        CF_BootstrapData.pV_opt=zeros(pp,pp,NoGroups,B_REPS_current);

        CF_BootstrapData.Shocks = zeros([pz, size(Y,2), size(Y,3), B_REPS_current]);
        
        CF_BootstrapData.Y_b_first = zeros([size(Y), B_REPS_current]);
        CF_BootstrapData.Y_b_second= zeros([size(Y), B_REPS_current]);
        CF_BootstrapData.pY_b_first= zeros([size(pY), B_REPS_current]);
        CF_BootstrapData.pY_b_second = zeros([size(pY), B_REPS_current]);
        
        
        CF_BootstrapData.Y_b_first_no_sp = zeros([size(Y), B_REPS_current]);
        CF_BootstrapData.Y_b_second_no_sp= zeros([size(Y), B_REPS_current]);
        CF_BootstrapData.pY_b_first_no_sp= zeros([size(pY), B_REPS_current]);
        CF_BootstrapData.pY_b_second_no_sp = zeros([size(pY), B_REPS_current]);

        CF_BootstrapData.IV_IRF_opt=zeros(p+pp,pv1,NoGroups,B_REPS_current);
        CF_BootstrapData.IV_IRF_iid=zeros(p+pp,pv1,NoGroups,B_REPS_current);

        CF_BootstrapData.IV_F_opt=zeros(pz+pv1,NoGroups,B_REPS_current);

        BS_BLOCK_LENGTH=1;
        BlocksIntoObs = ones(T-k-MissingObs,EffectiveObs);
        BlockToObsMap = @(block_indx) t;
        for bb=1:EffectiveObs 
            BlocksIntoObs(:,bb)= ceil( (1:T-k-MissingObs)' / BS_BLOCK_LENGTH )==bb;
        end
    end

    DO_TEST=1;

%     for b_idx=1:B_REPS_current
        SIGN_FLIP=ones( EffectiveObs,N );
        IDX=( 1:(T-k-MissingObs))' * ones(1,N);

        Y_b=NaN(p,(T-k-MissingObs),N);
        Z_b=NaN(pz,(T-k-MissingObs),N);
        X_b=NaN(px,(T-k-MissingObs),N);
        IV_Y_b=zeros( p_IV ,T-MissingObs-NoIVLags-1,N);    
    %         dX_b=NaN(px,(T-k-MissingObs)-1,N);

        pY_b=NaN(pp,(T-k-MissingObs),N);
        pX_b=NaN(ppx,(T-k-MissingObs),N,pp);
        IV_Yp_b=zeros( p_IVp ,T-MissingObs-NoIVLags-1,N);    
    %         dXp_b=NaN(ppx,(T-k-MissingObs)-1,N,pp);

        % Initialize: For the first k observations, the data is just data...
    %     Y_b(:, 1:k, :) = DATA.Y(:, 1:k, :);
        X_b(:, 1, :)  = DATA.X(:, 1, :);
    %         dX_b(:, 1:k, :) = DATA.dX(:, 1:k, :);
        pX_b(:, 1, :,:)  = DATA.pX(:, 1, :,:);
    %         dXp_b(:, 1:(1+k)*pp, :,:) = DATA.dXp(:, 1:(1+k)*pp, :,:);
    %     dY_b(:, 1:k, :) = DATA.dY(:, 1:k, :);
    %     IV_Y_b(:, 1:NoIVLags, :) = DATA.IV_Y(:, 1:NoIVLags , :);    

    %         ErrorTerm=zeros(p+pz,N);    
        ErrorTerm=zeros(p+pz,N);
        pErrorTerm=zeros(pp,N);
    %         ErrorTerm=zeros(p+pz,N);    
    
        DO_SPATIAL=1;

        if b_idx<=2
            my_cf_construction;

            sum(abs(reshape(Y-Y_b,1,[])))
            sum(abs(reshape(pY-pY_b,1,[])))
        end
        
        % without 2nd shock;
        shock_E_opt=CF_E_opt*0;
        old_CF_E_opt=CF_E_opt;
        old_CF_pE_opt=CF_E_opt;
        for gg=1:NoGroups
            for t_idx = 1 : T-MissingObs-k
                shock_E_opt(1:p,t_idx ,IndexM==gg)=beta{gg} \ squeeze(old_CF_E_opt(1:p,t_idx ,IndexM==gg));
                CF_BootstrapData.Shocks(1:pz, t_idx, find(IndexM==gg), b_idx) = shock_E_opt(1:pz,t_idx ,IndexM==gg);
%                 CF_BootstrapData.Shocks(1, t_idx, find(IndexM==gg), b_idx) = shock_E_opt(1,t_idx ,IndexM==gg);
                shock_E_opt(pz,t_idx ,IndexM==gg)=shock_E_opt(pz,t_idx ,IndexM==gg)*0;
                CF_E_opt(1:p,t_idx ,IndexM==gg)=beta{gg} * squeeze(shock_E_opt(1:p,t_idx ,IndexM==gg));
            end
        end

        my_cf_construction;

        Demeaned_Y_b= Y_b - repmat(mean(Y_b,2) ,[1,T-MissingObs-k,1]); % year FE
        Demeaned_Y_b= Demeaned_Y_b - repmat(mean(Demeaned_Y_b,3) ,[1,1,N]);
        
        Demeaned_pY_b= pY_b - repmat(mean(pY_b,2) ,[1,T-MissingObs-k,1]); % year FE
        Demeaned_pY_b= Demeaned_pY_b - repmat(mean(Demeaned_pY_b,3) ,[1,1,N]);
        
        CF_BootstrapData.Y_b_second(:,:,:,b_idx) = Demeaned_Y_b;
        CF_BootstrapData.pY_b_second(:,:,:,b_idx)= Demeaned_pY_b;


        CF_E_opt=old_CF_E_opt;
        CF_E_opt=old_CF_pE_opt;

        % without first shock;
        shock_E_opt=CF_E_opt*0;
        old_CF_E_opt=CF_E_opt;
        old_CF_pE_opt=CF_E_opt;
        for gg=1:NoGroups
            for t_idx = 1 : T-MissingObs-k
                shock_E_opt(1:p,t_idx ,IndexM==gg)=beta{gg} \ squeeze(old_CF_E_opt(1:p,t_idx ,IndexM==gg));
%                 CF_BootstrapData.Shocks(pz, t_idx, find(IndexM==gg), b_idx) = shock_E_opt(pz,t_idx ,IndexM==gg);
                shock_E_opt(1,t_idx ,IndexM==gg)=shock_E_opt(1,t_idx ,IndexM==gg)*0;
                CF_E_opt(1:p,t_idx ,IndexM==gg)=beta{gg} * squeeze(shock_E_opt(1:p,t_idx ,IndexM==gg));
            end
        end
        my_cf_construction;

        Demeaned_Y_b= Y_b - repmat(mean(Y_b,2) ,[1,T-MissingObs-k,1]); % year FE
        Demeaned_Y_b= Demeaned_Y_b - repmat(mean(Demeaned_Y_b,3) ,[1,1,N]);
        
        Demeaned_pY_b= pY_b - repmat(mean(pY_b,2) ,[1,T-MissingObs-k,1]); % year FE
        Demeaned_pY_b= Demeaned_pY_b - repmat(mean(Demeaned_pY_b,3) ,[1,1,N]);
        
        CF_BootstrapData.Y_b_first(:,:,:,b_idx) = Demeaned_Y_b;
        CF_BootstrapData.pY_b_first(:,:,:,b_idx)= Demeaned_pY_b;

%**************************************************************************
% Switch of spatial channel!
%**************************************************************************
        DO_SPATIAL=0;

        % without 2nd shock;
        shock_E_opt=CF_E_opt*0;
        old_CF_E_opt=CF_E_opt;
        old_CF_pE_opt=CF_E_opt;
        for gg=1:NoGroups
            for t_idx = 1 : T-MissingObs-k
                shock_E_opt(1:p,t_idx ,IndexM==gg)=beta{gg} \ squeeze(old_CF_E_opt(1:p,t_idx ,IndexM==gg));
                shock_E_opt(pz,t_idx ,IndexM==gg)=shock_E_opt(pz,t_idx ,IndexM==gg)*0;
                CF_E_opt(1:p,t_idx ,IndexM==gg)=beta{gg} * squeeze(shock_E_opt(1:p,t_idx ,IndexM==gg));
            end;
        end;        
        
        my_cf_construction;

        Demeaned_Y_b= Y_b - repmat(mean(Y_b,2) ,[1,T-MissingObs-k,1]); % year FE
        Demeaned_Y_b= Demeaned_Y_b - repmat(mean(Demeaned_Y_b,3) ,[1,1,N]);
        
        Demeaned_pY_b= pY_b - repmat(mean(pY_b,2) ,[1,T-MissingObs-k,1]); % year FE
        Demeaned_pY_b= Demeaned_pY_b - repmat(mean(Demeaned_pY_b,3) ,[1,1,N]);
        
        CF_BootstrapData.Y_b_second_no_sp(:,:,:,b_idx) = Demeaned_Y_b;
        CF_BootstrapData.pY_b_second_no_sp(:,:,:,b_idx)= Demeaned_pY_b;


        CF_E_opt=old_CF_E_opt;
        CF_E_opt=old_CF_pE_opt;

        % without first shock;
        shock_E_opt=CF_E_opt*0;
        old_CF_E_opt=CF_E_opt;
        old_CF_pE_opt=CF_E_opt;
        for gg=1:NoGroups
            for t_idx = 1 : T-MissingObs-k
                shock_E_opt(1:p,t_idx ,IndexM==gg)=beta{gg} \ squeeze(old_CF_E_opt(1:p,t_idx ,IndexM==gg));
                shock_E_opt(1,t_idx ,IndexM==gg)=shock_E_opt(1,t_idx ,IndexM==gg)*0;
                CF_E_opt(1:p,t_idx ,IndexM==gg)=beta{gg} * squeeze(shock_E_opt(1:p,t_idx ,IndexM==gg));
            end;
        end;
        my_cf_construction;

        Demeaned_Y_b= Y_b - repmat(mean(Y_b,2) ,[1,T-MissingObs-k,1]); % year FE
        Demeaned_Y_b= Demeaned_Y_b - repmat(mean(Demeaned_Y_b,3) ,[1,1,N]);
        
        Demeaned_pY_b= pY_b - repmat(mean(pY_b,2) ,[1,T-MissingObs-k,1]); % year FE
        Demeaned_pY_b= Demeaned_pY_b - repmat(mean(Demeaned_pY_b,3) ,[1,1,N]);
        
        CF_BootstrapData.Y_b_first_no_sp(:,:,:,b_idx) = Demeaned_Y_b;
        CF_BootstrapData.pY_b_first_no_sp(:,:,:,b_idx)= Demeaned_pY_b;
        
%     end;
    display(b_idx);
    
    if NoGroups>1
        error('year FE currently not group specific');
    end;
    % Demeaned_Y_b_cf = Y_b - repmat(mean(Y_b,2) ,[1,T-MissingObs-k,N]); % year FE
    % Demeaned_Y_b_cf = Demeaned_Y_b_cf - repmat(mean(Demeaned_Y_b_cf,3) ,[1,1,N]);
    % 
    % Demeaned_Y_b_actual = Y - repmat(mean(Y,2) ,[1,T-MissingObs-k,N]); % year FE
    % Demeaned_Y_b_actual= Demeaned_Y_b_actual- repmat(mean(Demeaned_Y_b_actual,3) ,[1,1,N]);
end;
%%
Demeaned_Y= Y - repmat(mean(Y,2) ,[1,T-MissingObs-k,1]); % year FE
Demeaned_Y= Demeaned_Y - repmat(mean(Demeaned_Y,3) ,[1,1,N]);

Demeaned_pY= pY - repmat(mean(pY,2) ,[1,T-MissingObs-k,1]); % year FE
Demeaned_pY= Demeaned_pY - repmat(mean(Demeaned_pY,3) ,[1,1,N]);

PlotScale= [VarScale(1:p), PerScale(1:pp)];
PlotNames  = [VarNames(1:p); PerNames(:)];
PlotLabels = [VarLabels(1:p); PerLabels(:)];


for target_idx=1:length(TargetMSA)
    close all; 
    for p_idx=1:p+pp
        if p_idx<=p
            MSA_idx=find(vMSA==TargetMSA{target_idx});
            YData=prctile(squeeze(CF_BootstrapData.Y_b_first(p_idx,:,MSA_idx,:)), [16 50 84],2 );
            
            YData2=prctile(squeeze(CF_BootstrapData.Y_b_second(p_idx,:,MSA_idx,:)), [16 50 84],2 ); 
        else
            YData= prctile(squeeze(CF_BootstrapData.pY_b_first(p_idx-p,:,MSA_idx,:)), [16 50 84],2 );

            YData2= prctile(squeeze(CF_BootstrapData.pY_b_second(p_idx-p,:,MSA_idx,:)), [16 50 84],2 ); 

        end;
        
        if PlotScale(p_idx)==1
            DisplayUnit='p.p.';
        else
            DisplayUnit='%';
        end;
        
        hndl=my_make_IRF(vYear(MissingObs+k+1:T),YData,YData2, 'DoAreaPlot', true, 'YLabel',DisplayUnit,'XLabel','years');
        hold all
        if p_idx<=p
            hndl_data=plot( vYear(MissingObs+k+1:T), Demeaned_Y(p_idx,:,MSA_idx) ); 
        else
            hndl_data=plot( vYear(MissingObs+k+1:T), Demeaned_pY(p_idx-p,:,MSA_idx) ); 
        end;
        set(gcf, 'name', [PlotNames{p_idx} '::' num2str(TargetMSA{target_idx})]);
        set(hndl_data, 'Color', 'b', 'LineWidth', 2);
%         plot(xlim, [0 0], '-k');
        axis([vYear([MissingObs+k+1,T])'+[-1, 1] ylim]);
        ylabel(DisplayUnit);
        legend([hndl_data;hndl([1 3])],{'data';'no startup shock';'no lab. dmd shock'}, 'location', 'best')
        legend('boxoff');
        set(findall(gcf,'type','axes'),'fontsize',16);
        set(findall(gcf,'type','text'),'fontSize',16);   
    %     axis([min(vYear) EndYear ylim]);
        grid on;
        
        print('-depsc', ['../Graphs/Hist_' PlotNames{p_idx} Suffix num2str(TargetMSA{target_idx},'_%1.0f') '.eps']);
 
    end;
    
    % plot shock history -- works only with 2 shocks
    YData= prctile(squeeze(CF_BootstrapData.Shocks(1,:,MSA_idx,:)), [16 50 84],2 );
    YData2= prctile(squeeze(CF_BootstrapData.Shocks(2,:,MSA_idx,:)), [16 50 84],2 );

    DisplayUnit='st.devs.';

    hndl=my_make_IRF(vYear(MissingObs+k+1:T),YData,YData2, 'DoAreaPlot', true, 'YLabel',DisplayUnit,'XLabel','years');
    hold all
%         plot(xlim, [0 0], '-k');
    axis([vYear([MissingObs+k+1,T])'+[-1, 1] ylim]);
    ylabel(DisplayUnit);
    legend([hndl([2 3])],{'startup shock';'lab. dmd shock'}, 'location', 'best')
    legend('boxoff');
    set(findall(gcf,'type','axes'),'fontsize',16);
    set(findall(gcf,'type','text'),'fontSize',16);   
%     axis([min(vYear) EndYear ylim]);
    grid on;

    print('-depsc', ['../Graphs/Hist_shock_' Suffix num2str(TargetMSA{target_idx},'_%1.0f') '.eps']);
 
end;
%% same, but no spatial effects
Demeaned_Y= Y - repmat(mean(Y,2) ,[1,T-MissingObs-k,1]); % year FE
Demeaned_Y= Demeaned_Y - repmat(mean(Demeaned_Y,3) ,[1,1,N]);

Demeaned_pY= pY - repmat(mean(pY,2) ,[1,T-MissingObs-k,1]); % year FE
Demeaned_pY= Demeaned_pY - repmat(mean(Demeaned_pY,3) ,[1,1,N]);

PlotScale= [VarScale(1:p), PerScale(1:pp)];
PlotNames  = [VarNames(1:p); PerNames(:)];
PlotLabels = [VarLabels(1:p); PerLabels(:)];

for target_idx=1:length(TargetMSA)
    close all;
    for p_idx=1:p+pp
        if p_idx<=p
            MSA_idx=find(vMSA==TargetMSA{target_idx});
            YData=prctile(squeeze(CF_BootstrapData.Y_b_first_no_sp(p_idx,:,MSA_idx,:)), [16 50 84],2 );
            
            YData2=prctile(squeeze(CF_BootstrapData.Y_b_second_no_sp(p_idx,:,MSA_idx,:)), [16 50 84],2 ); 
        else
            YData= prctile(squeeze(CF_BootstrapData.pY_b_first_no_sp(p_idx-p,:,MSA_idx,:)), [16 50 84],2 );

            YData2= prctile(squeeze(CF_BootstrapData.pY_b_second_no_sp(p_idx-p,:,MSA_idx,:)), [16 50 84],2 ); 

        end;
        
        if PlotScale(p_idx)==1
            DisplayUnit='p.p.';
        else
            DisplayUnit='%';
        end;
        
        hndl=my_make_IRF(vYear(MissingObs+k+1:T),YData,YData2, 'DoAreaPlot', true, 'YLabel',DisplayUnit,'XLabel','years');
        hold all
        if p_idx<=p
            hndl_data=plot( vYear(MissingObs+k+1:T), Demeaned_Y(p_idx,:,MSA_idx) ); 
        else
            hndl_data=plot( vYear(MissingObs+k+1:T), Demeaned_pY(p_idx-p,:,MSA_idx) ); 
        end;

        set(hndl_data, 'Color', 'b', 'LineWidth', 2);
%         plot(xlim, [0 0], '-k');
        axis([vYear([MissingObs+k+1,T])'+[-1, 1] ylim]);
        ylabel(DisplayUnit);
        legend([hndl_data;hndl([1 3])],{'data';'no startup shock, R=0';'no lab. dmd shock, R=0'}, 'location', 'best')
        legend('boxoff');
        set(findall(gcf,'type','axes'),'fontsize',16);
        set(findall(gcf,'type','text'),'fontSize',16);   
    %     axis([min(vYear) EndYear ylim]);
        grid on;
        
        print('-depsc', ['../Graphs/Hist_no_sp_' PlotNames{p_idx} Suffix num2str(TargetMSA{target_idx},'_%1.0f') '.eps']);

    end
end
%%
DATA=OLD_DATA;
A_IV=OLD_A_IV;
Ap_IV=OLD_Ap_IV;

save(['tmp_results' Suffix '_cf'],'CF_BootstrapData', 'Demeaned_*')

MSA_Shocks=reshape(prctile(CF_BootstrapData.Shocks(1,:,:,:),50, 4), 28, 354)';
ToWrite=[num2cell(vMSA), cMSA, num2cell(MSA_Shocks)];
filename = ['../Data/ShockSeries' Suffix '_cf2.xlsx'];
xlswrite(filename,ToWrite)

MSA_Shocks=reshape(mean(CF_BootstrapData.Shocks(1,:,:,:),4), 28, 354)';
ToWrite=[num2cell(vMSA), cMSA, num2cell(MSA_Shocks)];
filename = ['../Data/MeanShockSeries' Suffix '_cf2.xlsx'];
xlswrite(filename,ToWrite)