%
% Estimate quarterly finanical shocks using quarterly gross output 
%

load shocks.out
load gypc.dat  
newphi    = shocks(4:43,2:5);  %1975-2014
gypc      = gypc(13:172,2:5);
nt        = 160;
n         = 2;
for i=1:4;
  qvar    = gypc(:,i);
  avar    = newphi(:,i);
  tqvar   = hptrend(log(qvar),1600);
  tavar   = hptrend(log(avar),100);
  Y       = zeros(nt,2);
  Y(:,1)  = log(qvar)-tqvar;
  ii      = 4:4:160;
  Y(ii,2) = log(avar)-tavar;
  a2qvar  = ann2qtr([[1975:2014]',log(avar)-tavar],4,1);
  a2qvar  = a2qvar(:,2);
  qY      = [Y(:,1),a2qvar];
  sumyx   = zeros(n,4*n+1);
  sumxx   = zeros(4*n+1,4*n+1);
  for t=5:160;
    Ylag  = [qY(t-1,:)';qY(t-2,:)';qY(t-3,:)';qY(t-4,:)';1];
    sumyx = sumyx+qY(t,:)'*Ylag';
    sumxx = sumxx+Ylag*Ylag';
  end;
  A       = sumyx*inv(sumxx);
  Sig     = zeros(n);
  for t=5:160;
    Ylag  = [qY(t-1,:)';qY(t-2,:)';qY(t-3,:)';qY(t-4,:)';1];
    eps   = qY(t,:)'-A*Ylag;
    Sig   = Sig+eps*eps'/156;
  end;
  A       = A(:,1:4*n);
  B       = chol(Sig)';
  Theta0  = [A(:);B(:)];

  [Theta1,f,g]  = uncmin(Theta0,'lestqtr',Y)
  if abs(sum(imag(Theta1)))>0;
    disp('Warning: imaginary numbers found in Theta')
    Y2    = estqtr(real(Theta1),Y);
  else
    Y2    = estqtr(Theta1,Y);
  end;
  qtrend  = ann2qtr([[1975:2014]',tavar],4,1);
  qtrend  = qtrend(:,2);
  qpred   = exp(Y2(:,n)+qtrend);
  Phi(:,i)= qpred;
end;

save qtrly
