/* ----- PROGRAM c:\racd\racd5p2.gau July 1998 ------ This GAUSS program does Poisson PMLE for data in Chapter 5 of A.C. Cameron and Pravin K. Trivedi (1998) REGRESSION ANALYSIS OF COUNT DATA, Econometric Society Monograph No.30, Cambridge University Press. A title such as Table 3.3, col.3, p.69 means Table 3.2 column 3 page 69 For more details on the data see file racd5.1st 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 racd5.asc and in gauss data set racd5.dat on 126 firms, 14 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 the data set racd5.dat includes 13. SIZESQ Size Squared 14. CONSTANT Unity This program performs Poisson PMLE 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 = "c:\\racd\\racd5p2.out"; output file = ^outfile reset; print "File racd5p2.out is output from Gauss program racd5p2.gau"; print "Date is:"; print date; /******* Read in the data ********/ data ="c:\\racd\\racd5"; print; call dstat(data,0); @ check that data okay @ print; let YVAR = NUMBIDS; let XVAR = CONSTANT,LEGLREST,REALREST,FINREST,WHTKNGHT,BIDPREM, INSTHOLD,SIZE,SIZESQ,REGULATN; /* Set data size dimensions and some other variables Will need to change nobs in different examples */ 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); /******** 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 */ /* 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 */ /* Begin newton-raphson iterations */ cha = 1; n = nobs; 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 mle standard errors and t-statistics Several ways are given: (1) Usual mle using minus inverse of hessian (2) Robust sandwich or Eicker-White or Huber (3) Corrected mle assuming variance multiple of mean This program then uses (3) */ covusual = invpd(-n*hes); @ (1) usual ml estimate of covariance @ opg = (x.*(y-mu))'(x.*(y-mu)); @ outer product of the gradient @ covrobus = covusual*opg*covusual; @ (2) robust sandwich/ Huber/ Eicker-White @ alphanb1 = (1/(n-k))*(y-mu)'((y-mu)./mu); covnb1 = alphanb1*covusual; @ (3) nb1 estimate using nb1 variance function @ bmle = b; covmle = covnb1; @ Change this line for alternative variance estimators @ varmle = diag(covmle); semle = sqrt(varmle); tmle = bmle./semle; /* print results */ print; print "Cameron & Trivedi REGRESSION ANALYSIS OF COUNT DATA Table 5.3 p.148"; print; print "Poisson PML regression"; print "Standard errors and t-statistics assume NB1 variance function"; print; print $yvar "is the dependent variable"; print; print "Variable Poisson PMLE Standard errors t-statistics"; results=ones(k,4); results[.,1]=xvar; results[.,2:4]=bmle~semle~tmle; let resmask[1,4]= 0 1 1 1; @ character numeric numeric numeric @ call printfmt(results,resmask); /* End of program with procedures following below */ /* Procedures opendat, rad, readdat ********************************** Following procedures written 4/16/97 based on Per's fsa2.src. In turn I think this is based on Frank Windmeijer's code. Basically given a Gauss data set with names Calling opendat, rad, readdat in seldata.src 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 ----------------------------------------------------------------------- Open the dataset and get the indices of the explanatory and dependent variables { f1,n,indy,indx } = opendat(data,yvar,xvar,&indices,m); INPUT: data - String data set to be used yvar - Character vecor of names of dependent variables xvar - Charcter vector of names of explantory variables &indices - Proc Indices.src, should be included in the main program m - The Number of dependent variables OUTPUT: f1 - File handle for the open data set n - Number of observations in the data set indy - Index for the dependent variables indx - Index for the exogenous variables ------------------------------------------------------------------- { nror, rest } = rad(nobs,n); INPUT: n - number of observations in the data set nobs - number of observations to read OUTPUT: nror - number of times data set needs to be read to be complete y - number of obervations in the time data is read ----------------------------------------------------------------------- Reads the nobs observations of the data set { x,y } = Readdat(f1,r,nobs,indy,indx); INPUT: f1 - File handle for the open data set r - the row where to start the reading nobs - number of observations to read indy - Index for the dependent variables indx - Index for the exogenous variables program OUTPUT: x - nobs*k matrix of explanatory variables y - nobs*m matrix of dependent variables ----------------------------------------------------------------------- */ proc (4) = opendat(data,yvar,xvar,&index,m); local f1,n,p,var; local depvar,indy,ind,indx; local index:proc; p = rows(yvar|xvar); { var,ind } = index(data,yvar|xvar); indy = ind[1:m]; indx = ind[m+1:p]; open f1 = ^data; n = rowsf(f1); retp(f1,n,indy,indx); endp; proc (2) = rad(nobs,n); local a,nror,rest; if n > nobs; a = n/nobs; nror = ceil(a); rest = nobs*(a- (nror-1)); else; nror = 1; rest = n; endif; retp(nror,rest); endp; proc (2) = Readdat(f1,r,nobs,indy,indx); local y,z,x; y = seekr(f1,r); z = readr(f1,nobs); x = z[.,indx]; y = z[.,indy]; retp(y,x); endp;