$ontext


    Variance estimator for an MLE. Use Convert to retrieve hessian.

    Erwin Kalvelagen July 2008

    Reference:
       William H. Greene, "Econometric Analysis", 5th edition, 2003

$offtext

$version 228

set i 'cases' /i1*i20/;

table data(i,*)

                Y         E
     i1         20.5      12
     i2         31.5      16
     i3         47.7      18
     i4         26.2      16
     i5         44.0      12
     i6         8.28      12
     i7         30.8      16
     i8         17.2      12
     i9         19.9      10
    i10         9.96      12
    i11         55.8      16
    i12         25.2      20
    i13         29.0      12
    i14         85.5      16
    i15         15.1      10
    i16         28.5      18
    i17         21.4      16
    i18         17.7      20
    i19         6.42      12
    i20         84.9      16
;

parameter y(i),x(i);
y(i) = data(i,'y');
x(i) = data(i,'E');

scalar Var2Greene 'Variance estimate number two from Greene' / 46.16337 /;

variable b,L;
equation loglik;

loglik.. L =e= sum(i,log(b+x(i))) + sum(i,y(i)/(b+x(i)));


*
* estimation
*
model estimation /loglik/;
solve estimation minimizing L using nlp;
display "estimate",b.l;


*
* get hessian
*
option nlp=convert;
$onecho > convert.opt;
hessian
$offecho
estimation.optfile=1;
solve estimation minimizing L using nlp;

*
* gams cannot add elements at runtime so we declare the necessary elements here
*
set dummy /e1,x1/;

parameter h(*,*,*) '-hessian';
execute_load "hessian.gdx",h;
display h;

*
* inverting a 1x1 matrix is not very difficult
*
scalar v 'variance';
v = -1/h('e1','x1','x1');
display v,var2greene;

abort$(abs(v-var2greene)>1.0e-5) "Accuracy check failed";