%RUNLIKE   Use an example from CKM (with the baseline parameters) to
%          test the loglike.m program
%

%          Ellen McGrattan

nsim  = 1;
randn('seed',100462)
rand('seed',912)

%
% Theoretical VAR derived in RBC2BMS:
%
icase = 0;

if icase==1;
  alpha = .99;
  B1    = [0.01060538395901   0.03825012538435
           0.00924502646712  -0.04898282860064];
  M     = [0.96492795689229  -0.00725239676071
          -0.02182852871038   0.98369651108134];
else;
  alpha = 0;
  B1    = [0.01060538395901   0.03825012538435
           0.00924502646712   0.94101717139936];
  M     = [0.99440047322473  -0.04106164036378
           0.86689350016928  -0.03579633543679];
end;
Sige  = [0.06611001841241  -0.04808330083838
        -0.04808330083838   0.21893351828209]*1.0e-03 *100^2;
cSige = chol(Sige)';


%
% Simulate data
%

A     = B1+M;
B     = M;
Tdrop = 100;
L     = zeros(nsim);
Est   = zeros(11,nsim);
dL    = zeros(11,nsim);
imp   = zeros(13,nsim);

for n=1:nsim;
  T       = 400;
  eps     = randn(T,2);
  Y(1,:)  = zeros(1,2);
  e       = (cSige*eps')';
  for i=2:T;
    Y(i,:)= (A*Y(i-1,:)'+e(i,:)'-B*e(i-1,:)')';
  end;
  Y       = Y(Tdrop+1:T,:);

  %
  % Use same data as loglike.f90
  %
  load ydat.dat
  Y       = ydat;

  T       = length(Y);

  %
  % Call minimization routine with some initial guess for the parameters
  %

  A0      = A;
  B0      = B;
  cSige0  = cSige;
  Gam     = [vec([A0;B0]);cSige0(find(tril(ones(2))))];
  p       = 1;
  q       = 1;
  [T,k]   = size(Y);
  data    = [1;1;T;2;Y(:)];
  [x,L,dL,info,stat]= uncmin(Gam,'loglike',data)
  Est(:,n)= x;
  lnl(:,n)= L;
  dll(:,n)= dL;

  %
  % Use estimates in vector x to construct B1, M, Sige 
  %

  r       = max(p,q+1);
  AB      = zeros((p+q)*k,k);
  AB(:)   = x(1:(p+q)*k*k);
  Ass     = [[AB(1:p*k,:);zeros((r-p)*k,k)],[eye((r-1)*k);zeros(k,(r-1)*k)]];
  Bss     = [eye(k);-AB(p*k+1:(p+q)*k,:);zeros((r-1-q)*k,k)];
  Css     = [eye(k),zeros(k,(r-1)*k)];
  Mhat    = AB(p*k+1:(p+q)*k,:);
  B1hat   = AB(1:p*k,:)-Mhat;
  Sighat  = tril(ones(k));
  i       = find(Sighat);
  Sighat(i)= x(k*k*(p+q)+1:length(x));
  Sighat  = Sighat*Sighat';

  %
  % Calculate a B-Q impulse response:
  %
  imp(:,n)= calcimp(B1hat,Mhat,Sighat,alpha,13)

  %
  % Statistics for VARMA(2,1) on Z = [ log Output; log Hours]
  %
  A1z     = ([1,1;0,1]*[1,-1;0,alpha]+[1,1;0,1]*(B1+M)*[1,-1;0,1]);
  A2z     = -[1,1;0,1]*(B1+M)*[1,-1;0,alpha];
  Bz      = [1,1;0,1]*M;
  Sigz    = [1,1;0,1]*Sige*[1,1;1,0];

  % call program to generate population statistics,
  %   e.g., covariogram: cov(Z[t],Z[t+h]) = sum C_j Sigz C_j+h'
  %         variances conditioned on one shock off (look at thesis)
  % small sample:
  %   hpfiltered simulations of y,h 
  %   variance decomp 

end;
imp0    = calcimp(B1,M,Sige,alpha,13);
