% FAKE_GR   Determines if Gali-Rabanal can uncover the theoretical 
%           impulse response working with data in fakedat.dat.
%           Edit this code to change options.
%

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

%
% User Options:
%   nt     = number of quarters of response to compute
%
%   p      = number of lags in VAR
%
%   alpha in [0,1] with 0 implying hours in levels, 
%                       1 implying hours differenced
%                      <1 implying quasi-differenced
%
%
p     = 4;
alpha = 1;
filter= 2;

%
% Load data on hours and output
%
load fakedat.dat
prd  = fakedat(:,1);
hpc  = fakedat(:,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 
%
[X1,X2,epsd,epsz,A0,B0,B] = bquah(x,p,13);
for i=1:p;
  eval(['B',int2str(i),'=B(:,2*(',int2str(i),'-1)+1:2*',int2str(i),');']);
end;
v      = [zeros(p,2); [epsd,epsz]*A0'];
vz     = [zeros(p,2); [zeros(nobs-1-p,1),epsz]*A0'];
vd     = [zeros(p,2); [epsd,zeros(nobs-1-p,1)]*A0'];
xn     = x;
xz     = x;
xd     = x;
for i=p+1:nobs-1;
  sum1 = [0;0];
  sum2 = [0;0];
  sum3 = [0;0];
  for j=1:p;
    eval([' sum1 = sum1 + B',int2str(j),'*xn(i-j,:)'' ;'])
    eval([' sum2 = sum2 + B',int2str(j),'*xz(i-j,:)'' ;'])
    eval([' sum3 = sum3 + B',int2str(j),'*xd(i-j,:)'' ;'])
  end;
  xn(i,:) = ( B0+sum1+v(i,:)' )';
  xz(i,:) = ( B0+sum2+vz(i,:)' )';
  xd(i,:) = ( B0+sum3+vd(i,:)' )';
end;
for i=1:nobs-1;
  lh0(i,1)  = (alpha.^[i-1:-1:0])*x(1:i,2)*100;
  lh(i,1)   = (alpha.^[i-1:-1:0])*xn(1:i,2)*100;
  lhz(i,1)  = (alpha.^[i-1:-1:0])*xz(1:i,2)*100;
  lhd(i,1)  = (alpha.^[i-1:-1:0])*xd(1:i,2)*100;
end;
ly0    = cumsum(x(:,1))*100+lh0;
ly     = cumsum(xn(:,1))*100+lh;
lyz    = cumsum(xz(:,1))*100+lhz;
lyd    = cumsum(xd(:,1))*100+lhd;
if filter==1;
  vy0  = var(ly0-hptrend(ly0,1600));
  vy   = var(ly-hptrend(ly,1600));
  vyz  = var(lyz-hptrend(lyz,1600));
  vyd  = var(lyd-hptrend(lyd,1600));
  vh0  = var(lh0-hptrend(lh0,1600));
  vh   = var(lh-hptrend(lh,1600));
  vhz  = var(lhz-hptrend(lhz,1600));
  vhd  = var(lhd-hptrend(lhd,1600));
else;
  vy0  = var(baxter(ly0,6,32,8));
  vy   = var(baxter(ly,6,32,8));
  vyz  = var(baxter(lyz,6,32,8));
  vyd  = var(baxter(lyd,6,32,8));
  vh0  = var(baxter(lh0,6,32,8));
  vh   = var(baxter(lh,6,32,8));
  vhz  = var(baxter(lhz,6,32,8));
  vhd  = var(baxter(lhd,6,32,8));
end;


