function [IR,CIlb,CIub] = impres(mdata,alpha,p,nt,nrep);
% IMPRES    Computes impulse responses to technology and boot-strapped 
%           confidence bands given data on productivity and per-capita hours.
%           [IR,CIlb,CIub]=impres(mdata,alpha,p,nt,nrep) takes
%           productivity and hours per capita data in the Tx2 matrix 
%           mdata, along with options alpha, p, nt, and nrep and
%           applies the Blanchard-Quah SVAR procedure. The options are:
%             * alpha in [0,1] with 0 implying hours in levels, 
%                                   1 implying hours differenced
%                                  <1 implying quasi-differenced
%             * p      = number of lags in VAR
%             * nt     = number of quarters of response to compute
%             * nrep   = number of replications for the bootstrapped bands.
%           If nrep=1, the confidence bands are not computed.
% 

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

prd    = mdata(:,1);
hpc    = mdata(:,2);

%
% 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 
%
if nrep==1;
  [X1,X2,epsd,epsz] = bquah(x,p,nt);
  prdz   = cumsum(X2(:,1))*100;
  for i=1:nt;
    hrsz(i,1)  = (alpha.^[i-1:-1:0])*X2(1:i,2)*100;
  end;
  IR     = [prdz,hrsz];
  CIlb   = 0;
  CIub   = 0;
else
  [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')';
  IR     = [prdz,hrsz];
  lb     = round(.025*nrep);
  ub     = round(.975*nrep);
  if lb>0;
    CIub   = [sepz(:,ub),sehz(:,ub)];
    CIlb   = [sepz(:,lb),sehz(:,lb)];
  else
    CIub   = [sepz(:,nrep),sehz(:,nrep)];
    CIlb   = [sepz(:,1),sehz(:,1)];
  end;
end;  
