%IRTheory   plots the theoretical impulse response generated by
%           truth.m (IR) along with population moments derived from
%           applying a Blanchard-Quah SVAR procedure with long
%           run restrictions imposed.

%           Ellen McGrattan, 2-13-05

nvar   = input('Number of variables? ');
nlag   = input('Number of VAR lags? ');
Sig    = D0*Q*Q'*D0';
Cbar   = inv(eye(nvar)-inv(eye(nvar)-(C2*inv(C1)-C1))*C1);
A0chol = inv(Cbar)*chol(Cbar*Sig*Cbar')';
M      = B2*inv(B1);

V0     = dlyap(C1+M,Sig+C1*Sig*C1'-(C1+M)*Sig*(C1+M)');
V1     = (C1+M)*V0-M*Sig;
VA     = zeros(nlag*nvar);
Vb     = zeros(nlag*nvar,nvar);
for i=2:nlag;
  eval(['V',int2str(i),'=(C1+M)^(',int2str(i-1),')*V1;'])
end;
for j=1:nlag
  for i=1:j-1;
    eval(['VA((i-1)*nvar+1:i*nvar,(j-1)*nvar+1:j*nvar) = V',int2str(j-i),';'])
  end;
  VA((j-1)*nvar+1:j*nvar,(j-1)*nvar+1:j*nvar) = V0;
  for i=j+1:nlag
    eval(['VA((i-1)*nvar+1:i*nvar,(j-1)*nvar+1:j*nvar) = V',int2str(i-j),''';'])
  end;
  eval(['Vb((j-1)*nvar+1:j*nvar,:) = V',int2str(j),''';'])
end;
Bols   = (inv(VA)*Vb)';
Sigols = V0-Vb'*inv(VA)*Vb;
sum1   = eye(nvar);
for i=1:nlag
  sum1 = sum1 - Bols(:,(i-1)*nvar+1:i*nvar);
end;  
Cbols  = sum1\eye(nvar);
A0ols  = inv(Cbols)*chol(Cbols*Sigols*Cbols')';

nt     = length(IR);
IRhld  = zeros(nt,nvar);
IRhld(1,:) = A0ols(:,1)';
Clag   = [eye(nvar),zeros(nvar,nvar*(nlag-1))];
for i=2:nt
  C          = zeros(nvar);
  for j=1:nlag;
    C        = C+Bols(:,(j-1)*nvar+1:nvar*j)*Clag(:,(j-1)*nvar+1:nvar*j);
  end;
  tem  = C*A0ols;
  eval(['A',int2str(i-1),'ols=C*A0ols;'])
  Clag = [C,Clag(:,1:nvar*(nlag-1))];
  IRhld(i,:) = tem(:,1)';
end;
for i=1:nt;
  IRsvar(i,2) = (alpha.^[i-1:-1:0])*IRhld(1:i,2)*100;
end;
IRsvar(:,1) = cumsum(IRhld(:,1))*100;

