/* gauss\src\ghcount.src  
   Procedures for maxlik estimation of Poisson regression model.
   Based on pp1.src by Per Johansson for JAE 1997 paper.
   It provides ....

   polli(theta,f1)                  Log-likelihood function for Poisson MLE
   pogrd(theta,f1)                  Gradient for Poisson MLE
   pohess(theta,indy,indx,f1,nobs)  Hessian for Poisson MLE

   Note polli returns an nobs x 1 vector  (not a scalar)
        pogrd returns an nobs x k matrix  (not a k x 1 vector) 
        pohess returns a k x k matrix     (as usual) 
      
   The program assumes that data in a nobs x p matrix z are ordered as
      The first column of z is the dependent variable y.
      The remaining columns of z are regressors
*/

/* (1) Log-likelihood function for Poisson */
proc(1) = polli(theta,z);
local li,x,y,p;
p = cols(z);
    y = z[.,1];
    x = z[.,2:p];
li = y.*(x*theta) - exp(x*theta);
retp(li);
endp;

/* Gradient for Poisson */
proc(1) = pogrd(theta,z);
local gr,p,x,y,pdf;
p = cols(z);
    y = z[.,1];
    x = z[.,2:p];
gr = x.*(y-exp(x*theta));
retp(gr);
endp;

/* Analytical hessian for poisson */
proc(1) = pohess(theta,z);
local lambda,p,x,H,y;
p = cols(z);
    y = z[.,1];
    x = z[.,2:p];
lambda = exp(x*theta);
H = -(x.*lambda)'x;
retp(H);
endp;

/* The remaining code is for numerical derivatives if needed 
   It is a refinement of a Gauss program hessp.src
   For the Poisson this code is not necessary.
*/
/*
** hessperrc								419
** rewritten Gauss procedure
** Purpose:    Computes the matrix of second partial derivatives
**             (Hessian matrix) of a function defined by a procedure.
**
** Format:     h = HESSPER( &f, x0, z);
**
** Inputs:     &f -- pointer to a single-valued function f(x), defined
**                   as a procedure, taking a single Kx1 vector
**                   argument (f:Kx1 -> 1x1).  It is acceptable for
**                   f(x) to have been defined in terms of global
**                   arguments in addition to x:
**
**                      proc f(x);
**                          retp( exp(x'b) );
**                      endp;.
**
**           x0 -- Kx1 vector specifying the point at which the Hessian
**                 of f(x) is to be computed.
**           z  -- data matrix or file handle   (30/8/95)
** Output:   h -- KxK matrix of second derivatives of f with respect
**                to x at x0. This matrix will be symmetric.
**
** Remarks:    This procedure requires K(K+1)/2 function evaluations. Thus
**             if K is large it may take a long time to compute the
**             Hessian matrix.
**
**             No more than 3 - 4 digit accuracy should be expected from
**             this function, though it is possible for greater accuracy
**             to be achieved with some functions.
**
**             It is important that the function be properly scaled, in
**             order to obtain greatest possible accuracy. Specifically,
**             scale it so that the first derivatives are approximately
**             the same size.  If these derivatives differ by more than a
**             factor of 100 or so, the results can be meaningless.
**
*/

proc hessper(&f,x0,z);
    local k, hessian, grdd, ax0, dax0, dh, xdh, ee, f0, i, j, eps;
    local f:proc;

    /* check for complex input */
    if iscplx(x0);
        if hasimag(x0);
            errorlog "ERROR: Not implemented for complex matrices.";
            end;
        else;
            x0 = real(x0);
        endif;
    endif;

/* initializations */
    k = rows(x0);
    hessian = zeros(k,k);
    grdd = zeros(k,1);
    eps = 6.0554544523933429e-6;

/* Computation of stepsize (dh) */
    ax0 = abs(x0);
    if x0 /= 0;
        dax0 = x0./ax0;
    else;
        dax0 = 1;
    endif;
    dh = eps*maxc((ax0~(1e-2)*ones(rows(x0),1))').*dax0;

    xdh = x0+dh;
    dh = xdh-x0;    /* This increases precision slightly */
    ee = eye(k).*dh;

/* Computation of f0=f(x0) */
    f0 = f(x0,z);

/* Compute forward step */
    i = 1;
    do until i > k;
        grdd[i,1] = f(x0+ee[.,i],z);
        i = i+1;
    endo;

/* Compute "double" forward step */
    i = 1;
    do until i > k;
        j = i;
        do until j > k;
            hessian[i,j] = f(x0+(ee[.,i]+ee[.,j]),z);
            if i /= j;
                hessian[j,i] = hessian[i,j];
            endif;
            j = j+1;
        endo;
        i = i+1;
    endo;

    retp( ( ( (hessian - grdd) - grdd') + f0) ./ (dh.*dh') );
endp;

