/* 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;