/* ----- PROGRAM c:\racd\racd7p2.gau July 1998 ------ This GAUSS program does Poisson PMLE for data in Chapter 7 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 racd7.1st The data is the same data as originally used in J. Kennan, "The Duration of Contract strikes in U.S. Manufacturing", Journal of Econometrics, 1985, Vol. 28, pp.5-28. The data in the ascii file racd7.asc and in gauss data set racd7.dat for 108 months, 4 variables per month. The data is monthly U.S. data from 1968(1) to 1976(12) 1. OUTPUT is level of economic activity (measured as cyclical departure of aggregate production from its trend level) 2. STRIKES is number of strikes (number of contract strikes in U.S. manufacturing beginning each month) 3. TIME is a time trend from 1 to 108 4. CONSTANT is unity This program performs Poisson PMLE regression of STRIKES on ten regressors including the intercept: intercept,OUTPUT */ /******* Initial setup: libraries used etc. and output file ********/ New; @ clear workspace @ outfile = "c:\\racd\\racd7p2.out"; output file = ^outfile reset; print "File racd7p2.out is output from Gauss program racd7p2.gau"; print "Date is:"; print date; /******* Read in the data ********/ data ="c:\\racd\\racd7"; print; call dstat(data,0); @ check that data okay @ print; let YVAR = STRIKES; let XVAR = CONSTANT,OUTPUT; /* 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 = 108; @ 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-rapson 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 7.2 p.232"; 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 provedures 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;