% Runmleq   Compute MLE estimates for the benchmark model using the
%           sample period 1959:1--2003:4.   The likelihood function
%           is in file mleq.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.12218441160373
   0.23456126785859
   0.21813418432846
  -1.71044434741722
   0.76890460108055
  -0.02328852467396
  -0.10128969659480
  -0.03063606505969
   0.04705520990984
   0.99352262916894
   0.02531186131043
   0.04233819795798
   0.43178571903678
   0.04171548961226
   1.18316739565328
  -0.00160108354708
  -0.04188834170829
  -0.00452014405740
  -0.01798869654933
   0.99730018412890
   0.01081927622854
  -0.00256942082981
   0.00370721649763
  -0.00500766071435
  -0.00623233345783
  -0.00088560323078
  -0.00505372754194
   0.00194528573172
   0.01479167606597
  -0.00000249657496];

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

ind    = [1:30]';
param  = 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(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(' ')

                                                                                
