%% Compute variance decomposition


H_GRID=[0 10 100]; % horizon of decomposition
OutPutTable=zeros(p+pp,pv1+1, NoGroups, B_REPS,length(H_GRID));
for b_idx=1:B_REPS
    for gg=1:NoGroups
        K_IV_b=[BootstrapData.A_IV(:,:,gg,b_idx)'; eye( (k-1)*p,p*k ) ];
        Q1=chol(BootstrapData.V_opt(1:p,1:p,gg,b_idx), 'lower')\BootstrapData.IV_IRF_opt(1:p,:,gg,b_idx);
        Q=[Q1 null(Q1')];
        A_periphery=BootstrapData.Ap_IV(:,:,gg,b_idx);
        B=[eye(p,p+pp); A_periphery( k+ (1:p) ,:)', eye(pp)]*[chol(BootstrapData.V_opt(:,:,gg,b_idx), 'lower')*Q, zeros(p,pp); zeros(pp,p) chol(BootstrapData.pV_opt(:,:,gg,b_idx), 'lower')];

        OldRli=@(l) eye(p+pp,(p+pp)*k) * K_IV_original{gg}^l * [ B; zeros((k-1)*(p+pp),p+pp)];
        h_min=0;
        oldS=zeros(p+pp,p+pp,pv1+2, length(H_GRID));
        S=zeros(p+pp,p+pp,p+pp, length(H_GRID));
        h_count=0;
        for h_count=1:length(H_GRID)
            HH=H_GRID(h_count);
            for h_idx=0:HH
%                 h_min = HH;
                CurrentRli=OldRli(h_idx);
                for pz_idx=1:pv1
                    oldS(:,:,pz_idx,h_count)=oldS(:,:,pz_idx,h_count)+(HH+1-max(h_idx,h_min))*CurrentRli(:,pz_idx)*CurrentRli(:,pz_idx)';
                end;
                oldS(:,:,pv1+1,h_count)=oldS(:,:,pv1+1,h_count)+(HH+1-max(h_idx,h_min))*CurrentRli(:,pv1+1:p)*CurrentRli(:,pv1+1:p)';
                oldS(:,:,pv1+2,h_count)=oldS(:,:,pv1+2,h_count)+(HH+1-max(h_idx,h_min))*CurrentRli(:,p+1:p+pp)*CurrentRli(:,p+1:p+pp)';

            end;
            for pz_idx=1:pv1+2
                OutPutTable(:,pz_idx,gg, b_idx, h_count)=100*diag(squeeze(oldS(:,:,pz_idx,h_count)))./diag(squeeze(sum(oldS(:,:,:,h_count),3)));
            end;


        end;
    end;
end;

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

for alpha=[10 32]
    for h_count=1:length(H_GRID)

        for gg=1:NoGroups
            HH=H_GRID(h_count);
            MeanOutPutTable=squeeze(mean(OutPutTable(:,:,gg,:,h_count),4));
            OutPutTable_p5=prctile(OutPutTable(:,:,gg,:,h_count),alpha/2,4);
            OutPutTable_p95=prctile(OutPutTable(:,:,gg,:,h_count),100-alpha/2,4);
            CellTable=cellstr([char(PlotLabels ) ...
                num2str(MeanOutPutTable(:,1), '& %1.1f') repmat(' ', p+pp,1) num2str(OutPutTable_p5(:,1), '& (%1.1f') num2str(OutPutTable_p95(:,1), ', %1.1f)') ... 
                num2str(MeanOutPutTable(:,2), '& %1.1f') repmat(' ', p+pp,1) num2str(OutPutTable_p5(:,2), '& (%1.1f') num2str(OutPutTable_p95(:,2), ', %1.1f)') ... 
                num2str(MeanOutPutTable(:,pv1+1), '& %1.1f') num2str(OutPutTable_p5(:,pv1+1), '& (%1.1f') num2str(OutPutTable_p95(:,pv1+1), ', %1.1f)') ... 
                num2str(MeanOutPutTable(:,pv1+2), '& %1.1f') num2str(OutPutTable_p5(:,pv1+2), '& (%1.1f') num2str(OutPutTable_p95(:,pv1+2), ', %1.1f)') ... 
                repmat(' \\', p+pp,1)]);


            my_fid=fopen(['../Tables/TheoryVariance_BStrap_' Suffix num2str(gg, '%1.0f') '_HH' num2str(HH) num2str(alpha,'_alpha%1.0f') '.tex'],'w');
            fprintf(my_fid, '%s\n',CellTable{:});
            fclose(my_fid);
            display(num2str(HH, 'Horizon %1.0f'));
            display(CellTable(:));
        end;
    end;
end;