% FIGUREF1.M 

clear;

% Load posterior draws
load posterior_arrw0; z=za; clear za;

% Shells for impulse responses of admissible structural models
IRF11=zeros(size(z,1),H+1); IRF12=zeros(size(z,1),H+1); IRF13=zeros(size(z,1),H+1); 
IRF21=zeros(size(z,1),H+1); IRF22=zeros(size(z,1),H+1); IRF23=zeros(size(z,1),H+1); 
IRF31=zeros(size(z,1),H+1); IRF32=zeros(size(z,1),H+1); IRF33=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,:)=irf(1,:);   IRF21(j,:)=irf(2,:);   IRF31(j,:)=irf(3,:); 
    IRF12(j,:)=irf(4,:);   IRF22(j,:)=irf(5,:);   IRF32(j,:)=irf(6,:);  
    IRF13(j,:)=irf(7,:);   IRF23(j,:)=irf(8,:);   IRF33(j,:)=irf(9,:);

end

% Stack posterior draws for impulse responses
IRF=[IRF11 IRF12 IRF13 IRF21 IRF22 IRF23 IRF31 IRF32 IRF33]; 

clear A B0inv irf IRF11 IRF12 IRF13 IRF21 IRF22 IRF23 IRF31 IRF32 IRF33;

% 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));
IRF12bayes=irfbayes(1,1*(H+1)+1:2*(H+1));
IRF13bayes=irfbayes(1,2*(H+1)+1:3*(H+1));
IRF21bayes=irfbayes(1,3*(H+1)+1:4*(H+1));
IRF22bayes=irfbayes(1,4*(H+1)+1:5*(H+1));
IRF23bayes=irfbayes(1,5*(H+1)+1:6*(H+1));
IRF31bayes=irfbayes(1,6*(H+1)+1:7*(H+1));
IRF32bayes=irfbayes(1,7*(H+1)+1:8*(H+1));
IRF33bayes=irfbayes(1,8*(H+1)+1:9*(H+1));

IRF11jcs=jcs(:,0*(H+1)+1:1*(H+1));
IRF12jcs=jcs(:,1*(H+1)+1:2*(H+1));
IRF13jcs=jcs(:,2*(H+1)+1:3*(H+1));
IRF21jcs=jcs(:,3*(H+1)+1:4*(H+1));
IRF22jcs=jcs(:,4*(H+1)+1:5*(H+1));
IRF23jcs=jcs(:,5*(H+1)+1:6*(H+1));
IRF31jcs=jcs(:,6*(H+1)+1:7*(H+1));
IRF32jcs=jcs(:,7*(H+1)+1:8*(H+1));
IRF33jcs=jcs(:,8*(H+1)+1:9*(H+1));

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

% Rescale responses of first variable to facilitate interpretation of IRF plot
IRF11jcs=IRF11jcs*100; IRF11bayes=IRF11bayes*100;
IRF12jcs=IRF12jcs*100; IRF12bayes=IRF12bayes*100;
IRF13jcs=IRF13jcs*100; IRF13bayes=IRF13bayes*100;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot structural impulse responses 
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('Real price of gasoline','fontsize',15)
title('Expectations 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('Expectations 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])
ylabel('Inflation expectation','fontsize',15)
title('Expectations shock','fontsize',15)

subplot(3,3,4); plot(0:1:H,IRF12jcs,'r-',0:1:H,IRF12bayes,'k-',0:1:H,zeros(1,H+1),'k-','linewidth',3); axis([0 H  -9 9])
ylabel('Real price of gasoline','fontsize',15)
title('Nominal gas price shock','fontsize',15)
subplot(3,3,5); plot(0:1:H,IRF22jcs,'r-',0:1:H,IRF22bayes,'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,6); plot(0:1:H,IRF32jcs,'r-',0:1:H,IRF32bayes,'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)

subplot(3,3,7); plot(0:1:H,IRF13jcs,'r-',0:1:H,IRF13bayes,'k-',0:1:H,zeros(1,H+1),'k-','linewidth',3); axis([0 H  -9 9])
ylabel('Real price of gasoline','fontsize',15)
title('Core CPI shock','fontsize',15)
subplot(3,3,8); plot(0:1:H,IRF23jcs,'r-',0:1:H,IRF23bayes,'k-',0:1:H,zeros(1,H+1),'k-','linewidth',3); axis([0 H -0.1 0.25])
ylabel('Headline inflation','fontsize',15)
title('Core CPI shock','fontsize',15)
subplot(3,3,9); plot(0:1:H,IRF33jcs,'r-',0:1:H,IRF33bayes,'k-',0:1:H,zeros(1,H+1),'k-','linewidth',3); axis([0 H -0.1 0.5])
ylabel('Inflation expectation','fontsize',15)
title('Core CPI shock','fontsize',15)

% Impact effect inflation expectations
10*IRF32jcs(1,1)./(IRF12jcs(1,1)+IRF22jcs(1,1))
min(10*IRF32jcs(:,1)./(IRF12jcs(:,1)+IRF22jcs(:,1)))
max(10*IRF32jcs(:,1)./(IRF12jcs(:,1)+IRF22jcs(:,1)))

% Impact effect inflation
10*IRF22jcs(1,1)./(IRF12jcs(1,1)+IRF22jcs(1,1))
min(10*IRF22jcs(:,1)./(IRF12jcs(:,1)+IRF22jcs(:,1)))
max(10*IRF22jcs(:,1)./(IRF12jcs(:,1)+IRF22jcs(:,1)))

save posterior_output IRF11jcs IRF12jcs IRF13jcs IRF21jcs IRF22jcs IRF23jcs IRF31jcs IRF32jcs IRF33jcs;
