/* jagmle2.prg does maximum likelihood using Gauss procedure MAXLIK
using Jaggia data. It is more advanced than jagmle1.prg
TO RUN THIS PROGRAM
(1) You need an implementation of Gauss that has Maxlik (a Gauss add-on)
(2) You need my program SELDATA.SRC in the same directory (to select data)
(3) You need my program GHCOUNT.SRC in the same directory (for gradient etc.)
You can get seldata.src and ghpoiss.src from my web site.
Program jagmle2.prg
1. Reads in a flat ascii data set.
2. Creates new data from transformation of original data
3. Converts the data to a Gauss data set, with names.
4. Selects a subset of the data to be used in the OLS regression
(This uses my separate procedure seldata.src)
5. Calls my procedure ghpoiss.src to get logl, gradient and hessian
6. Calls the Gauss add-on maxlik to perform maximum likelihood.
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:
For description of the variables see jagols1.prg
In addition we need to create a constant and Size Squared
*/
/******* Initial setup: libraries used etc. and output file ********
*/
New; @ clear workspace @
library maxlik; @ So Gauss knows we use procedures in maxlik @
maxset; @ A maxlik command @
disable; @ Allows missing values in calculation @
outfile = "jagmle2.out"; @ all output goes to this file @
output file = ^outfile reset; @ delete earlier contents @
print "File jagmle2.out is output from Gauss program jagmle2.src";
print "Date is:";
print date;
print;
/******* Include additional Gauss code in separate files **********
seldata.src: Gauss program to select the subset of data
and variable names used in analysis.
This includes procedures opendat and readdat used below.
ghcount.src: Gauss program with gradient and hessian for poisson
This includes procedures polli, pogrd, pohess used below.
I assume you have these files in the same directory as this program.
*/
#include seldata.src;
#include ghcount.src;
@ or if e.g. in directory gauss then include c:\gauss\seldata.src @
/******* Read in the data and create variables **********
Here we first 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.
Then we create a Gauss data set with names provided.
*/
load data1[126,12] = jaggia.asc; @ Data read in to matrix @
one = ones(126,1); @ A column vector of ones @
sizesq = data1[.,7].*data1[.,7]; @ A column vector = seventh column squared @
data2 = data1~one~sizesq; @ New variables added to existing data @
/******* Create Gauss data set with names *********
The created Gauss data set will be called jaggiatr.dat
and the associated names will be in file jaggiatr.dht
Note that an alternative is to use the Gauss utility program ATOG386.EXE
(where ATOG means ascii to Gauss). See Gauss Systems and Graphics manual.
*/
dataset = "jaggiatr"; @ or e.g. "c:\\gauss\\myex\\jaggiatr" @
let vnames = DOCNO WEEKS NUMBIDS TAKEOVER BIDPREM INSTHOLD SIZE LEGLREST
REALREST FINREST REGULATN WHTKNGHT ONE SIZESQ;
call saved(data2,dataset,vnames);
/******** Extract the subset of the Gauss data set to be used *********
*/
data ="jaggiatr"; @ or e.g. "c:\\gauss\\myex\\jaggiatr" @
print;
call dstat(data,0); @ check data okay @
print;
let YVAR = NUMBIDS;
let XVAR = ONE LEGLREST REALREST FINREST WHTKNGHT BIDPREM INSTHOLD
SIZE SIZESQ REGULATN;
/* Set data size dimensions and some other variables
*/
m = 1; @ Number of dependent variables @
k = rows(xvar); @ Number of regressors @
nobs = 126; @ Maximum number of rows to read in dataset @
/* Call opendat, rad, readdat in seldata.src
This gives
x - nobs*k matrix of explanatory variables
y - nobs*m matrix of dependent variables
indy - Index for the dependent variables
indx - Index for the exogenous variables
f1 - File handle for the open data set
n - Number of observations in the data set
nror - number of times data set needs to be read to be complete
rest - number of obervations in the time data is read
The data are read into memory one time.
Nobs must be equal to number of observations in the data set
*/
{f1,n,indy,indx} = opendat(data,yvar,xvar,&indices,m);
{nror, rest } = rad(nobs,n);
{ y,x } = Readdat(f1,1,nobs,indy,indx);
/******** Estimation by maxlik version 4 *********
The log-likelihood, gradient and hessian procedures
used are polli, pogrd and pohess.
These appear in ghcount which was included earlier
*/
@ Some setting up not needed for maxlik @
z = y~x;
clear y,x;
parnames = XVAR;
output on;
" POISSON MAXIMUM LIKELIHOOD using Maxlik and ghmnl1.src";
" Dataset ";;Data;
" Nr of observations ";;n;
" Nr of times data is read ";;nror;
" Nr of obs. read last chew ";;rest;
" Dependent variable ";;$yvar;
" Parameters (Variables) names ";
$parnames';
"";
output off;
/* Analytical derivatives and hessian.
Can delete if instead use numerical derivatives
*/
_max_HessProc = &pohess; @ analytical hessian provided in hgcount.src @
_max_GradProc= &pogrd; @ analytical gradient provided in hgcount.src @
/* Setting various parameters for maxlik.
These are generally ones that over-ride the default values in maxlik.
For default values and descriptions see the maxlik manual
or the Gauss code in maxlik.src
*/
_max_Parnames = xvar; @ Parameter names which here = variable names @
_mlditer = 40; @ Not sure @
_max_Algorithm = 4; @ Newton Raphson = 4 BFGS = 2 (def) @
_max_MaxIters = 250; @ maximum number of iterations @
_mlmtime = 1e+5; @ Not sure 1e+5 = default @
_max_CovPar = 3; @ Method to compute covariance matrix @
@ 3 = heteroskedastic-consistent estimate @
_max_GradTol = 0.00001; @ Convergence tolerance for gradient of coeffs @
__output = 0; @ Minimal output. Use for Unix without X-windows @
start = (0.01*rndn(k,1)); @ start values for coeffs here randomly chosen @
/* Call maxlik.
Arguments are
data = the name of the Gauss data set (or matrix if data in matrix)
yvar|xvar = the variables to be selected from the Gauss data set
polli = procedure to compute log-likelihood (in ghcount.src)
start = starting values of the parameters
*/
output on;
{bp,phi0,grdp,covp,retcode} = maxlik(data,yvar|xvar,&polli,start);
/* Print the results from maxlik */
call maxprt(bp,phi0,grdp,covp,retcode);
output off;
f1 = close(f1);