/* jagmle1.prg does MLE with application to Poisson for Jaggia data.
   In practice we need more complicated data read in than this,
   and need to associate variables names with the output.
*/

/* The original data are from Sanjiv Jaggia and Satish Thosar, 1993,
 "Multiple Bids as a Consequence of Target Management Resistance"
  Review of Quantitative Finance and Accounting, 447-457.

  The data in the ascii file jaggia.asc are
  on 126 firms, 12 variables per firm:

  1. DOCNO     Doc No.
  2. WEEKS     Weeks
  3. NUMBIDS   Count (Dependent Variable)
  4. TAKEOVER  Delta (1 if taken over)
  5. BIDPREM   Bid Premium
  6. INSTHOLD  Institutional Holdings
  7. SIZE      Size measured in billions
  8. LEGLREST  Legal Restructuring
  9. REALREST  Real Restructuring
  10. FINREST  Financial Restructuring
  11. REGULATN Regulation
  12. WHTKNGHT White Knight

  In addition we need to create
  13. SIZESQ   Size Squared

  This program performs MLE regression of
     NUMBIDS
  on ten regressors including the intercept:
  intercept,LEGLREST,REALREST,FINREST,WHTKNGHT,BIDPREM,
  INSTHOLD,SIZE,SIZESQ,REGULATN
*/

/******* Initial setup: libraries used etc. and output file *******
*/

New;                              @ clear workspace @
outfile = "jagmle1.out";
output file = ^outfile reset;
print "File jagmle1.out is output from Gauss program jagmle.src";
print "Date is:";
print date;

/******* Read in the data and create variables              *******
 Here for simplicity we read the data in as a flat ascii file.
 The data are ordered by observation, with four variables
 per line so that the first observation (on 12 variables)
 takes up the first three lines etc.
 In practice we would use a Gauss data set with names provided.
*/

load data[126,12] = jaggia.asc;

DOCNO = data[.,1];
WEEKS = data[.,2];
NUMBIDS = data[.,3];
TAKEOVER = data[.,4];
BIDPREM = data[.,5];
INSTHOLD = data[.,6];
SIZE = data[.,7];
LEGLREST = data[.,8];
REALREST = data[.,9];
FINREST = data[.,10];
REGULATN = data[.,11];
WHTKNGHT = data[.,12];
SIZESQ = SIZE.*SIZE;
ONE = ones(126,1);

X = ONE~LEGLREST~REALREST~FINREST~WHTKNGHT~BIDPREM~
        INSTHOLD~SIZE~SIZESQ~REGULATN;
y = NUMBIDS;

/* Here X is an Nxk matrix and y is an Nx1 vector */

/* Check data by getting means */
n = rows(X);
xmean = X'one/n;
ymean = y'one/n;

print;
print "Sample means of the regressors:";
print xmean;
print "Sample mean of the dependent variable";
print ymean;

/******** Do the statistical analysis    *********/

/* Obtain starting values from ols regression of lny on x */
   b = invpd(X'X)*X'(ln(y+(y.==0)*0.1));

/* Estimate parameters using Newton-Raphson algorithm */

/* For the Poisson log-density
     ln L = Sum_i {-mu_i + y_i*ln(mu_i) - y_i!}
   where
     mu_i = exp(x_i'b)

   The first-order conditions are
     Sum_i (y_i - exp(x_i'b))x_i = 0

   and the second-derivative matrix is
     Sum_i -exp(x_i'b))x_i*x_i' = 0          */


/* To understand the code below you need to understand Gauss matrix operators.
   1.   * is the usual matrix multiplication.
   2.  .* is something new. It is element by element multiplication.
       For example if a is n x 1 and b is n x 1
       then a.*b is n x 1 with i-th entry a_i x b_i
       For example if A is n x k and b is n x 1
       then A.*b is n x k with i-th row a_i x b_i (a_i is i-th row of a)
   3.  .* is similarly define except element by element division.
       For example if A is n x k and b is n x 1
       then A./b is n x k with i-th row a_i / b_i (a_i is i-th row of a).


   In the code below sco is the first derivative vector
                                  divided by n (for numerical stability)
   So  sco = (X'(y-mu))/n
           = [ Sum_i (y_i - exp(x_i'b))x_i ] / n
   and hes is the second-derivative matrix
                                  divided by n (for numerical stability)
   So  hes = -X'(X.*mu)/n
           = [ Sum_i (-exp(x_i'b))x_i*x_i' ] / n                       */

    cha = 1;
    do until cha<1e-16;
        mu = exp(X*b);
        sco = (X'(y-mu))/n;
        hes = -X'(X.*mu)/n;
        bo = b;
        b = bo + invpd(-hes)*sco;
        cha = (bo-b)'(bo-b)/b'b;
    endo;

/* Obtain standard errors and t-statistics. 
   Following uses sandwich form with degrees of freedom adjustment */
bmle = b;
gradmatrix = X.*(y-mu);
covmle = invpd(-n*hes)*(gradmatrix'gradmatrix)*invpd(-n*hes)*(n/(n-cols(x)));
/* Do covmle = invpd(-n*hes); if instead want hessian form */
varmle = diag(covmle);
semle = sqrt(varmle);
tmle = bmle./semle;

/* print results */
print;
print "Poisson ML regression of NUMBIDS on intercept, LEGLREST, REALREST,
       FINREST,";
print "WHTKNGHT,BIDPREM,INSTHOLD,SIZE,SIZESQ and REGULATN";
print;
print "ML Poisson estimates, standard errors and t-statistics";
print "These are Huber-White sandwich standard errors with dof adjustment";
print bmle~semle~tmle;