% Runmleq   Compute MLE estimates for the benchmark model using Mexican
%           data over the sample period 1980:1--2003:4.   The likelihood 
%           function is in file mleq.m
%

%           Ellen McGrattan, 12-20-04
%           Revised, ERM, 1-4-05 

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

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

x0=[ ...
   0.43293923031769
   0.31117554608763
   0.03608516713915
  -2.44524344398936
   0.94403244804755
   0.95659677758810
   0.99544682509922
   0.87847354429560
   0.02280871328480
   0.01089889837824
  -0.01021730705886
  -0.06275278881480
  -0.02153949191098
   0.03815262712519
  -0.20695945468782
  -0.00719892655545
  -0.08660130991552
  -0.12364447778832];


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

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

[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',[x2,se]'))
disp(' ')
fprintf('  Log-Likelihood = %8.2f ',L)
disp(' ')

                                                                                
