Thursday, July 31, 2008

Hessian via PathNLP

Same model as Hessian via Convert but now using PathNLP.


$ontext


Variance estimator for an MLE

Erwin Kalvelagen July 2008

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

$offtext

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 via pathnlp
*

$onecho > pathnlp.opt
hessian_gdx 1
$offecho

option nlp=pathnlp;
model estimation /loglik/;
estimation.optfile=1;
solve estimation minimizing L using nlp;
display "estimate",b.l;


*
* get hessian
*
scalar d2l__b_b 'hessian';
execute_load "hesData.gdx",d2l__b_b;
display d2l__b_b;

*
* inverting a 1x1 matrix is not very difficult
*
scalar v 'variance';
v = 1/d2l__b_b;
display v,var2greene;

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


PathNLP prints all kind of debug output. This seems not completely polished yet.