% PREPHD_INFLEXP.M
clear;

% Load draws from 68% credible set
load models_posterior0;

% Shells for historical decompositions
num=size(z,1); 
yhat31=zeros(t-p,num); yhat32=zeros(t-p,num); yhat33=zeros(t-p,num); 

for j=1:num
j
    % Construct RF posterior parameters for jth admissible model   
    B      = (reshape(z(j,1:(1+n*p)*n),1+n*p,n))'; 
    B0inv  = reshape(z(j,(1+n*p)*n+n*n+1:(1+n*p)*n+2*n*n),n,n);

    % Compute structural multipliers for admissible model j
    IRF =irfvar([B(:,1+1:end);eye(n*(p-1)) zeros(n*(p-1),n)],B0inv,p,t-p-1);

    % Compute reduced-form residuals for modal model, given X and Ydep
    Uhat=Ydep-X*B';

    % Compute structural shocks Ehat from reduced form shocks Uhat
    Ehat=inv(B0inv)*Uhat';

    % Cross-multiply the weights for the effect of a given shock on the real
    % oil price (given by the relevant row of IRF) with the structural shock
    % in question
    for i=1:t-p
        yhat31(i,j)=dot(IRF(3,1:i),Ehat(1,i:-1:1));
        yhat32(i,j)=dot(IRF(6,1:i),Ehat(2,i:-1:1));
        yhat33(i,j)=dot(IRF(9,1:i),Ehat(3,i:-1:1));   

    end;
    
end; % and of j loop
clear IRF z X Ydep Ehat Uhat B0inv B n M rep;

% Time line for plots (adjust time axis later)
time=1981+7/12+p/12:1/12:2020+4/12;
yhat31=yhat31(1:end,:);
yhat32=yhat32(1:end,:);
yhat33=yhat33(1:end,:);

save hd_inflexp time p yhat31 yhat32 yhat33;
