function [X1,X2,epsd,epsz,se1_1,se1_2,se2_1,se2_2]=bquahboot(x,p,nt,nrep,frac)
%BQUAHBOOT   Compute impulse responses and confidence intervals for a
%            2 x 1 structural VMA with the identification scheme of 
%            Blanchard-Quah, that is
%
%                   X(t) = A(0)e(t)+A(1)e(t-1)+...            (1)
%
%            with the restrictions that E e e'= I, and sum_j a11(j)=0
%            where a11(j) is the upper left hand entry of A(j).
%            [X1,X2,epsd,epsz,se1_1,se1_2,se2_1,se2_2]=bquahboot(X,p,
%            nt,nrep,frac) requires a data matrix X with 2 columns, a 
%            choice for the number of lags p in the estimated VAR: 
%            (I-B1 L-...Bp L^p) X(t) = B0+ v(t), a choice nt 
%            determining the horizon for the impulse responses, 
%            a choice nrep for the number of replications in the 
%            bootstrapped standard errors, and a choice frac for 
%            the percentage confidence level (e.g., 95).  The responses 
%            to a ``demand shock'' are in X1 and the responses to a 
%            ``technology shock'' are in X2.  The empirical e(t)'s are 
%            in epsd and epsz.  All bootstrapped simulations are 
%            in the se matrices (with nrep columns each) if frac is
%            is set to 0.  The call bquahboot(X,p,nt) uses nrep=1000 
%            and frac=90 as defaults.  If standard errors are not 
%            required, use bquah.m.
%

%            Ellen McGrattan, 3-12-04

[X1,X2,epsd,epsz] = bquah(x,p,nt);

%
% Draw nrep times uniformly from T-p residuals
%
if nargin==3;
  nrep = 1000;
  frac = 95;
end;
[T,n]        = size(x);
int1         = [0:(1/(T-p)):((T-1-p)/(T-p))]';
int2         = [int1(2:T-p);1];
[B,se,resid] = autoreg(x,p);
xhat(1:p,:)  = x(1:p,:);
for i=1:nrep;
  for j=p+1:T;
    r        = rand;
    ind      = find(r>=int1 & r<int2);
    if isempty(ind); 
      disp('Warning: found no interval for ind')
      ind    = T-p; 
    end;
    u        = resid(ind,:);
    sum1     = 0;
    for k=1:p;
      sum1     = sum1+B(:,(k-1)*n+2:k*n+1)*xhat(j-k,:)';
    end;
    xhat(j,:)  = B(:,1)'+ sum1' + u;
  end;
  [x1,x2]    = bquah(xhat,p,nt);
  se1_1(:,i) = x1(:,1);
  se1_2(:,i) = x1(:,2);
  se2_1(:,i) = x2(:,1);
  se2_2(:,i) = x2(:,2);
end;
if frac > 0;
  cut        = (100-frac)/200;
  tem        = sort(se1_1');
  se1_1      = [tem(round(cut*nrep),:)',tem(round((1-cut)*nrep),:)'];  
  tem        = sort(se1_2');
  se1_2      = [tem(round(cut*nrep),:)',tem(round((1-cut)*nrep),:)'];  
  tem        = sort(se2_1');
  se2_1      = [tem(round(cut*nrep),:)',tem(round((1-cut)*nrep),:)'];  
  tem        = sort(se2_2');
  se2_2      = [tem(round(cut*nrep),:)',tem(round((1-cut)*nrep),:)'];  
end;

