% Runmleq   Compute MLE estimates for the benchmark model using the
%           sample period 1959:1--2003:4.   In this case, restrict
%           shocks to be orthogonal.  The likelihood function
%           is in file mleq_orth.m
%

%           Ellen McGrattan, 2-16-02
%           Revised, ERM, 3-16-04 

%------------------------------------------------------------------------------

load uszvarq.dat
global ZVAR
t    = uszvarq(:,1);
ZVAR = uszvarq(:,2:5);

x0=[ ...
   0.00423387
   0.30547259
   0.48495409
  -1.59023879
   0.99102181
   0.98909856
   0.99102472
   0.99102472
   0.00971733
   0.00691816
   0.00757662
   0.01637490];

[x1,f,g,code,status]   = uncmin(x0,'mleq_orth',0);
[L,Sbar,P0,P,Q]        = mleq_orth(x1);

ind   = [1,1,1,1, ...
         1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1, ...
         1,0,0,0,   1,0,0,     1,0,       1]';
param = zeros(30,1);
param(find(ind)) = x1;

[L,Lt] = mleseq(param,ind);

del = diag(max(abs(param)*1e-4,1e-8));
for i=1:length(param);
  if ind(i);
    [f1,f2] = mleseq(param+del(:,i),ind);
    [m1,m2] = mleseq(param-del(:,i),ind);
    dL(i,1)  = (f1-m1)/(2*del(i,i));
    dLt(i,:) = (f2-m2)'/(2*del(i,i));
  end;
end;
dLt =dLt(find(ind),:);
sum1=0;
[n,m]=size(dLt);
for t=1:m;
  sum1=sum1+dLt(:,t)*dLt(:,t)';
end;
se=diag(sqrt(inv(sum1)));

disp('Results from Maximum Likelihood Estimation')
disp('------------------------------------------')
disp(' ')
disp('  Parameters  Standard Errors ')
disp(sprintf('   %8.5f      %8.5f\n',[x1,se]'))
disp(' ')
fprintf('  Log-Likelihood = %8.2f ',L)
disp(' ')

                                                                                
