% VARDEC    Compute variances of output and hours attributable to
%           the demand and technology shocks.
%
%

%
% User Options:
%   data   = 1   if Gali-Rabanal's series used
%          = 2   if aggregate output and hours used
%
%   p      = number of lags in VAR
%
%   alpha in [0,1] with 0 implying hours in levels, 
%                       1 implying hours differenced
%                      <1 implying quasi-differenced
%
%   filter = 1   if Hodrick-Prescott filtered
%          = 2   if band-pass filtered
%
%
data   = 2;
p      = 4;
alpha  = 1;
filter = 2;

%
% 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 to get shocks and SVAR matrices
%
[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;

if (abs((vy-vy0)/vy)>1e-4 | abs((vh-vh0)/vh)>1e-4);
  error('Error: original and constructed data do not match')
end;
disp('Variance Decomposition:')
disp('-----------------------')
disp('                    Output     Hours   ')
disp(sprintf(' Total variance    %7.3f   %7.3f',vy0,vh0))
disp(sprintf(' Sum of variances  %7.3f   %7.3f',vyd+vyz,vhd+vhz))
disp(' Fractions due to: ')
disp(sprintf('   Demand          %7.3f   %7.3f',vyd/(vyd+vyz),vhd/(vhd+vhz)))
disp(sprintf('   Technology      %7.3f   %7.3f',vyz/(vyd+vyz),vhz/(vhd+vhz)))
disp(' ')
disp('NOTE: total variance and sum of variance differ because of filtering')

