/*
**  max4.e - A Poisson Model
**
**  In this example the heteroskedastic-consistent covariance matrix is
**  computed, and since the hessian and cross-product of the gradient are
**  calculated in the process of this computation they are also printed
**  out for comparison.  The poisson model estimated here is misspecified
**  because a random disturbance was added along with the exogenous
**  variables (see PSN.SIM for more discussion).  The correct standard
**  errors, even though the model is misspecified, may be computed from the
**  heteroskedastic-consistent covariance matrix.  When you run the model
**  note that the covariance matrix computed from the hessian is different
**  from that computed from the cross-product of the gradient (if the model
**  were correctly specified then they should be approximately the same).
**  Also note that the standard errors computed from the heteroskedastic-
**  consistent covariance matrix are larger than the standard errors from
**  either the hessian or the cross-product matrix (here called the
**  information matrix).  This reflects the fact that the true sampling
**  variance is larger because of the loss of efficiency by using a
**  misspecified model.  Note also that the standard errors from the
**  hessian (i.e., the usual standard errors) are under-estimated as a
**  result of the misspecification.
**
**    This example also illustrates the use of MAXLIK with a user-supplied
**  gradient procedure.
*/

library maxlik;
#include maxlik.ext;
maxset;

proc lpsn(b,z);
    local m;
    m = z[.,2:4]*b;
    retp(z[.,1].*m-exp(m));
endp;

proc lgd(b,z);
    retp((z[.,1]-exp(z[.,2:4]*b)).*z[.,2:4]);
endp;

x0 = { .5, .5, .5 };

_mlgdprc = &lgd;
_mlcovp = 3;

output file = max4.out reset;

    { x,f0,g,h,retcode } = maxlik("psn",0,&lpsn,x0);
    call maxprt(x,f0,g,h,retcode);
    print;
    print "heteroskedastic-consistent v-c matrix of par's";
    print h;
    print;
    print "hessian method v-c matrix of par's";
    print _mlhsvcp;
    print;
    print "cross product v-c matrix of par's";
    print _mlcpvcp;

output off;