function [IR IRu IRd] = var_ir(Ainv,B,Bs,Vs,shock,h,ci,tr)
% Compute impulse responses and bootstrap confidence bands
%   Input:  Ainv  = inverse of structural restriction matrix
%           B     = OLS estimates
%           Bs    = bs OLS estimates (enter scalar for no bs option)
%           Vs    = bs covariance estimates
%           shock = 1-by-k vector of structural shocks
%           h     = horizon for impulse response (default = 20)
%           ci    = confidence level (default = 0.95)
%           tr    = type of trend: const,linear,quadratic = (0,1,2)
%   Output: IR    = impulse response point estimate
%           IRu   = upper confidence band (0 if no bs option)
%           IRd   = lower confidence band (0 if no bs option)
%   Author: Kyle Jurado (kej2108@columbia.edu)

%% Options
if nargin == 5; h  = 20; ci = 0.95; tr = 0; end
if nargin == 6; ci = 0.95; tr = 0; end
if nargin == 7; tr = 0; end

%% Impulse responses
k  = size(B,2);                                     %number of variables
p  = (length(B)-(tr+1))/k;                          %number of lags
I  = [eye(k*(p-1)),zeros(k*(p-1),k)];               %block diagonal
C  = [B(tr+2:end,:)';I];                            %companion matrix
E  = [shock(:);zeros((p-1)*k,1)];                   %state space shocks
D  = blkdiag(Ainv,zeros((p-1)*k,(p-1)*k));          %structural matrix
IR = zeros(p*k,h);                                  %initialize irf
IR(:,1) = D*E;                                      %first response
for t = 2:h
    IR(:,t) = C*IR(:,t-1);                          %remaining responses
end
IR = IR(1:k,:)';                                    %remove st space part

%% Confidence bands (optional)
if (1-length(Bs))==0; IRu = 0; IRd = 0; end;
if (1-length(Bs))~=0;
S   = length(Bs);                                   %number of bs samples
IRb = zeros(p*k,h,S);                               %initialize array
for s = 1:S
    Cb    = [Bs{s}(tr+2:end,:)';I];                     %companion matrix
    Pb    = chol(Vs{s},'lower');
    Omb   = diag(diag(Pb));
    Ainvb = tril(Pb/Omb);    
    Db    = blkdiag(Ainvb,zeros((p-1)*k,(p-1)*k));  %structural matrix
    IRb(:,1,s) = Db*E;                              %first response
    for t = 2:h
        IRb(:,t,s) = Cb*IRb(:,t-1,s);               %remaining responses
    end
end
IRs = sort(IRb,3);                                  %sort impulse responses
IRu = IRs(:,:,ceil(S*(1+ci)/2));                    %upper band
IRd = IRs(:,:,ceil(S*(1-ci)/2));                    %lower ('down') band
IRu = IRu(1:k,:)'; IRd = IRd(1:k,:)';               %remove st space part
end;

end