% IRChk     Compute and plot empirical impulse responses for a SVAR
%           with productivity and per-capita hours using one of the
%           following datasets:
%
%             (1) US GDP and total hours (Prescott-Ueberfeldt)
%             (2) Francis-Ramey's series with adjusted hours per capita
%             (3) Christiano-Eichenbaum-Vigfusson's business sector
%             (4) Gali-Rabanal's series for nonfarm business sector
%
%           Also allow for alternative differencing of hours by
%           changing alpha in the code.
%

%           Ellen McGrattan, 3-13-04
%           Revised, ERM, 9-2-04

randn('seed',1004)
rand('seed',912)
%
% User Options:
%   data   = 1   if aggregate output and hours used (default)
%          = 2   if Francis-Ramey's series used
%          = 3   if Christiano-Eichenbaum-Vigfusson's series used
%          = 4   if Gali-Rabanal's series 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  = input(' Which dataset? ("help ir" shows 4 options) ');
if isempty(data);
  data  = 1;
end;
nt    = 13;
p     = 4;
nrep  = 1000;
alpha = 0;

%
% Load data on hours and output
%
if data==1;
  load ../data/BEA/NIPA/nipa116.csv
  load ../data/prescott/prescott.csv
  rGDP = nipa116(3,49:228)';
  Pop  = prescott(:,4);
  Hrs  = prescott(:,8);
  prd  = rGDP./Hrs*10^9;
  hpc  = Hrs./Pop;
  t    = ptime(1959.1,length(hpc),4);
end;
if data==2;
  load ../data/literature/fr/ramey.dat
  prd  = ramey(:,4);
  hpc  = ramey(:,3)./ramey(:,2);
  t    = ptime(1948.1,length(hpc),4);
end;
if data==3;
  load ../data/literature/cev/cev.dat
  prd  = cev(:,2);
  hpc  = cev(:,3)./cev(:,4);
  t    = ptime(1948.1,length(hpc),4);
end;
if data==4;
  load ../data/literature/gr/gali.csv
  prd  = gali(:,2)./gali(:,3); 
  hpc  = gali(:,3)./gali(:,4);
  t    = ptime(1948.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] = bquahbootchk(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')
