/*
**     max5.e
**
**  A Poisson Model with constraints on a parameter and an
**  illustration of the delta method used to compute its standard
**  error
**
**  In this example the two regression coefficients are constrained
**  to be equal. This is done because we know that they are equal in
**  the population (see psn.sim).
**
**  A procedure to compute the gradient could be included here as was done
**  in max3.e but it is simpler to use numerically computed gradients
**  since we will not have to take into account the constraint.
*/

library maxlik;
#include maxlik.ext;
maxset;

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

let a = { 1  0,
          0  1,
          0  1 };

proc cnstrn(b);
    retp(a*b);
endp;

let b0 = 1 .5;
_mlcovp = 3;

output file = max5.out reset;

    { b,f0,g,h,retcode } = maxlik("psn",0,&lpsn,b0);

    format /rd 1,0;
    print "algorithm: " _mlalgr;
    format /ro 10,5;
    print;
    print "function   " f0;
    print;
    print "unconstrained parameters " b';
    print;
    print "constrained parameters " cnstrn(b)';
    print;
    print "gradient   " g';
    print;

/*  To get the v-c matrix of the constrained parameters it is necessary
**  to compute the Jacobian of the constrained parameters with respect
**  to the unconstrained parameters.  In general the Jacobian can be
**  produced by calling GRADP with the constraint procedure and the
**  unconstrained parameters as arguments, for example,
**
**           j = gradp(&cnstrn,b);
**
**  However, in the present example the Jacobian is a constant and
**  is equal to "a" the matrix of zeros and ones above.
*/
    j = a;

    print "heteroskedastic-consistent v-c matrix of constrained par's";
    print j*h*j';
    print;
    print "hessian method v-c matrix of constrained par's";
    print j*_mlhsvcp*j';
    print;
    print "cross product v-c matrix of constrained par's";
    print j*_mlcpvcp*j';

output off;