% Runmleq_ur   Compute MLE estimates for the model with random walk 
%              technology using the sample period 1959:1--2003:4.   
%              The likelihood function is in file mleq_ur.m
%
                                                                                
%              Ellen McGrattan, 2-16-02
%              Revised, ERM, 3-16-04
                                                                                
%------------------------------------------------------------------------------
load uszvarq_ur.dat
global ZVAR
t    = uszvarq_ur(:,1);
ZVAR = uszvarq_ur(:,2:5);

x0 =[ ...
   3.9347e-03
   3.0121e-01
   9.2394e-01
   5.3643e-02
   1.8795e-01
   5.1438e-02
   9.7705e-01
  -8.4147e-02
  -2.6835e-02
   1.2658e-02
   1.0471e+00
  -6.5525e-03
  -1.9874e-03
  -1.9399e-02
  -1.0983e-02
  -1.2456e-02
   6.3435e-04
   1.7157e-05];

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

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

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

del = diag(max(abs(param)*1e-4,1e-8));
for i=1:length(param);
  if ind(i);
    [f1,f2] = mleseq_ur(param+del(:,i),ind);
    [m1,m2] = mleseq_ur(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(' ')

                                                                                
