% GR_IR     Compute and plot empirical impulse responses using either
%             (1) Gali-Rabanal's series for nonfarm business, or
%             (2) US GDP and total hours 
%
%           Also allow for alternative differencing of hours by
%           changing alpha in the code.
%

%           Ellen McGrattan, 3-13-04
%           Revised, ERM, 6-16-04

%
% User Options:
%   data   = 1   if Gali-Rabanal's series used
%          = 2   if aggregate output and hours used
%
%   nt     = number of quarters to plot 
%
%   p      = number of lags in VAR
%
%   nrep   = number of replications for boot-strapped confidence bands
%
%   alpha in [0,1] with 0 implying hours in levels, 
%                       1 implying hours differenced
%                      <1 implying quasi-differenced
%
%
data  = 2;
nt    = 13;
p     = 4;
nrep  = 1000;
alpha = 1;

%
% Load data on hours and output
%
if data==1;
  load gali.csv
  prd  = gali(:,2)./gali(:,3); 
  hpc  = gali(:,3)./gali(:,4);
  t    = ptime(1948.1,length(hpc),4);
else
  load nipa116.csv
  load prescott.csv
  rGDP = nipa116(2,:)';
  Pop  = prescott(:,4);
  Hrs  = prescott(:,8);
  prd  = rGDP./Hrs*10^9;
  hpc  = Hrs./Pop;
  t    = ptime(1959.1,length(hpc),4);
end;

%
% Construct the inputs to Blanchard-Quah
%
nobs   = length(hpc);
x      = log(prd);
n      = log(hpc);
x      = x-x(1);
n      = n-n(1);
dx     = x(2:nobs)-x(1:nobs-1);
dn     = n(2:nobs)-alpha*n(1:nobs-1); 
x      = [dx,dn];

%
% Call Blanchard-Quah routine which also generates boot-strapped errors
%
[X1,X2,epsd,epsz,se1_1,se1_2,se2_1,se2_2] = bquahboot(x,p,nt,nrep,0);
prdz   = cumsum(X2(:,1))*100;
sepz   = cumsum(se2_1)*100;
sepz   = sort(sepz')';
for i=1:nt;
 hrsz(i,1)  = (alpha.^[i-1:-1:0])*X2(1:i,2)*100;
 sehz(i,:)  = (alpha.^[i-1:-1:0])*se2_2(1:i,:)*100;
end;
sehz   = sort(sehz')';

%
% Plot results with 95% confidence bands
%
bnds   = [round(.025*nrep),round(.975*nrep)];
figure(1)
plot(0:nt-1,prdz,0:nt-1,sepz(:,bnds),'k--')
title('Response of Productivity to Technology Shock')

figure(2)
plot(0:nt-1,hrsz,0:nt-1,sehz(:,bnds),'k--')
title('Response of Hours to Technology Shock')
