% FIGUREIRF0_CHOL.M 

clear;

% Load posterior draws
load posterior0_chol; 

% Shells for impulse responses of admissible structural models
IRF11=zeros(size(z,1),H+1); IRF21=zeros(size(z,1),H+1); IRF31=zeros(size(z,1),H+1); 

for j=1:size(z,1)
        
    A      = (reshape(z(j,1:(1+n*p)*n),1+n*p,n))';   % RF level parameters
    B0inv  = reshape(z(j,(1+n*p)*n+n*n+1:(1+n*p)*n+2*n*n),n,n);    
    irf=irfvar([A(:,2:p*n+1); eye((p-1)*n,(p-1)*n) zeros((p-1)*n,n)],B0inv,p,H);	                    
    IRF11(j,:)=cumsum(irf(1,:)); IRF21(j,:)=irf(2,:); IRF31(j,:)=irf(3,:); 

end

% Stack posterior draws for impulse responses
IRF=[IRF11 IRF21 IRF31]; 

clear A B0inv irf;

% Bayes estimate of vector of inpulse responses for joint inference and 68% joint credible set
[irfbayes,jcs,index]=jointL1(IRF); clear IRF;    

IRF11bayes=irfbayes(1,0*(H+1)+1:1*(H+1));
IRF21bayes=irfbayes(1,1*(H+1)+1:2*(H+1));
IRF31bayes=irfbayes(1,2*(H+1)+1:3*(H+1));

IRF11jcs=jcs(:,0*(H+1)+1:1*(H+1));
IRF21jcs=jcs(:,1*(H+1)+1:2*(H+1));
IRF31jcs=jcs(:,2*(H+1)+1:3*(H+1));

z=z(index,:);
save models_posterior0_chol index M z p n H t Ydep X;
clear M rep z p n t Ydep X index j jcs irfbayes; 

IRF11bayes=IRF11bayes*100; IRF11jcs=IRF11jcs*100;

subplot(3,3,1); plot(0:1:H,IRF11jcs,'r-',0:1:H,IRF11bayes,'k-',0:1:H,zeros(1,H+1),'k-','linewidth',3); axis([0 H  -9 9])
ylabel('Nominal gasoline price','fontsize',15)
title('Nominal gas price shock','fontsize',15)
subplot(3,3,2); plot(0:1:H,IRF21jcs,'r-',0:1:H,IRF21bayes,'k-',0:1:H,zeros(1,H+1),'k-','linewidth',3); axis([0 H -0.1 0.25])
ylabel('Headline inflation','fontsize',15)
title('Nominal gas price shock','fontsize',15)
subplot(3,3,3); plot(0:1:H,IRF31jcs,'r-',0:1:H,IRF31bayes,'k-',0:1:H,zeros(1,H+1),'k-','linewidth',3); axis([0 H -0.1 0.5])
title('Nominal gas price shock','fontsize',15)
ylabel('Inflation expectation','fontsize',15)


% Impact effect on inflation expectations
10*IRF31jcs(1,1)./(IRF11jcs(1,1))
min(10*IRF31jcs(:,1)./(IRF11jcs(:,1)))
max(10*IRF31jcs(:,1)./(IRF11jcs(:,1)))


% Impact effect inflation
10*IRF21jcs(1,1)./(IRF11jcs(1,1))
min(10*IRF21jcs(:,1)./(IRF11jcs(:,1)))
max(10*IRF21jcs(:,1)./(IRF11jcs(:,1)))

