options double crt; name nbpan; ? ? Negative Binomial (I) Panel data models ? (fixed effects, random effects, etc.) ? Reproduces Limdep results for Table 9.1 of Cameron and Trivedi ? "Regression Analysis of Count Data", 1998, p.286. ? See files on Cameron's web page for rough corrections to the table, ? Limdep output with more/correct digits, and data file. ? http://www.econ.ucdavis.edu/~cameron ? Code adapted from hhgnb.tsp by Bronwyn H. Hall. ? by Clint Cummins 1/1999 ? ?File racd9d2.asc has 1730 lines = 346 firms times 5 years (1975-79) per firm ?For each firm in each year there is data on ? Time-invariant variables: OBSNO,YEAR,CUSIP,ARDSSIC,SCISECT,LOGK,SUMPAT ? Patents and up to 4 lags: PAT,PAT1,PAT2,PAT3,PAT4 ? LogR and up to 5 lags: LOGR,LOGR1,LOGR2,LOGR3,LOGR4,LOGR5 ?Each line is formatted as: ? 5F9.1,F12.5,6F8.1,6F12.5 const n,346 t,5; smpl 1,n; read(file='racd9d2.asc',format='(5F9.1,F12.5,6F8.1,6F12.5)') OBSNO,YEAR,CUSIP,ARDSSIC,SCISECT,LOGK,SUMPAT ?PAT,PAT1,PAT2,PAT3,PAT4 PAT75,no no no no ?LOGR,LOGR1,LOGR2,LOGR3,LOGR4,LOGR5 LOGR75,LOGR74 LOGR73 LOGR72 LOGR71 LOGR70 no no no no no no no PAT76,no no no no LOGR76,no no no no no no no no no no no no PAT77,no no no no LOGR77,no no no no no no no no no no no no PAT78,no no no no LOGR78,no no no no no no no no no no no no PAT79,no no no no LOGR79,no no no no no ; msd pat75-pat79 logr75-logr79 logk scisect ; ? Dependent variable is called n; compute constants for ? likelihoods. totaln = 0 ; dot 75-79 ; n.=pat. ; totaln = totaln+n. ; lconst. = lgamfn(n.+1) ; enddot ; sumconst = lconst75+lconst76+lconst77+lconst78+lconst79 ; lconst = lgamfn(totaln+1) ; msd totaln ; ? ? Construct the likelihood function for each model. ? list(prefix=Xb,first=75,last=79) Xbs; dot 75-79 ; ? ? Xb. is the Xb equation. Obviously one can ? add Z's (time-varying slopes) to this equation by including ? terms of the form d1.*z1, and so forth. ? Zd contains terms that do not vary over time (factor out of ? the lambda term). The full model for the expected number of ? patents is ? ? gamma(i,t) = exp(x(i,t)beta + z(i)delta) ? ? [In the version below the intercept varies over] ? [years but the slope coefficients do not. ] ? [frml Xb. a.+br*logr. ; ] ?frml Xb. a. + br*logr. ; param a.; frml lamb. exp(Xb.) ; frml gam. exp(Zd+Xb.) ; frml logl. lgamfn(gam. + n.) - lgamfn(gam.) ; eqsub logl. gam. ; frml loglc. lgamfn(lamb. + n.) - lgamfn(lamb.) ; eqsub loglc. lamb. ; enddot ; const a75,0; param br0 br1-br5; frml Xb75 a75 + br0*logr75 + br1*logr74 + br2*logr73 + br3*logr72 + br4*logr71 + br5*logr70; frml Xb76 a76 + br0*logr76 + br1*logr75 + br2*logr74 + br3*logr73 + br4*logr72 + br5*logr71; frml Xb77 a77 + br0*logr77 + br1*logr76 + br2*logr75 + br3*logr74 + br4*logr73 + br5*logr72; frml Xb78 a78 + br0*logr78 + br1*logr77 + br2*logr76 + br3*logr75 + br4*logr74 + br5*logr73; frml Xb79 a79 + br0*logr79 + br1*logr78 + br2*logr77 + br3*logr76 + br4*logr75 + br5*logr74; ? Sum of probabilities for denominator frml sumlam lamb75+lamb76+lamb77+lamb78+lamb79 ; eqsub sumlam lamb75-lamb79 ; frml sumlog logl75+logl76+logl77+logl78+logl79 ; eqsub sumlog logl75-logl79 ; frml sumlogc loglc75+loglc76+loglc77+loglc78+loglc79 ; eqsub sumlogc loglc75-loglc79 ; ? Model for non-time-varying variables. frml Zd b0 + bsize*logk + dsci*scisect ; param b0 bsize dsci ; frml expZd exp(Zd) ; frml sumgam expZd*sumlam; title 'NegBin no effects'; ? NegBin total - eq. (3.1) ? Note: it is much easier to reproduce these plain NegBin ? results by reading the data file as a cross section and ? using the NEGBIN(MODEL=1) command. frml leqt loglt = sumlog - (sumgam+totaln)*log(1+delta) + sumgam*log(delta) - sumconst ; eqsub leqt sumlog sumgam sumlam Xbs expZd Zd ; param delta .05; options nwidth=15,signif=8; ml(hiter=N,hcov=WN,tol=1e-7) leqt ; copy @logl ltotal ; copy @ncid ktotal ; mmake btotal a76-a79 br0-br5; ? save params for fixed effects starts title "NegBin Marginal (No Effects)" ; ? NegBin marginal - eq. (not in paper, same as 3.1 for the sum) ? Note that loglm is not exactly equal to loglt-loglc; full ? partioning is not possible in the negative binomial case. frml leqm loglm = lgamfn(sumgam+totaln) - lgamfn(sumgam) - lconst - (sumgam+totaln)*log(1+delta) + sumgam*log(delta) ; eqsub leqm sumgam sumlam Xbs expZd Zd ; ? year effects not identified in marginal const a76,0 a77,0 a78,0 a79,0; ml(hiter=N,hcov=WN,tol=1e-7) leqm ; copy @logl lmarg ; copy @ncid kmarg ; param a76-a79; title "NegBin Conditional -- Fixed Effects - Table 9.1" ; ? NegBin conditional - eq. (3.2) ? C-T numbers below are corrections from web page. ? C-T TSP ? lnR0 .3633 .363235562 ? lnR-1 .1555 .155526163 ? lnR-2 .1741 .174014901 ? lnR-3 .0147 .014852519 ? lnR-4 .0288 .028758446 ? lnR-5 .1360 .136057697 ? LogL -3391.35 -3391.35392 frml leqc loglc = sumlogc + lgamfn(sumlam) - lgamfn(sumlam+totaln) - sumconst + lconst ; eqsub leqc sumlogc sumlam Xbs ; unmake btotal a76-a79 br0-br5; ? use starts from total ml(hiter=N,hcov=WN,tol=1e-7) leqc ; copy @logl lcond ; copy @ncid kcond ; set testfe = 2*(lcond+lmarg-ltotal) ; set dfe = kcond+kmarg-ktotal ; title "Test for Presence of Firm Effects" ; cdf(df=dfe,chisq) testfe ; title "NegBin Total Model with Random Effects" ; ? C-T numbers below are from Limdep output file (not included in book). ? LIMDEP TSP ? B0 .8995721 .899561767 ? BSIZE .1619361 .161937013 ? DSCI .1176401 .117641997 ? BR0 .3503101 .350311917 ? BR1 -.003030878 -.00303172551 ? BR2 .1049859 .104987675 ? BR3 .01635281 .016352101 ? BR4 .03594328 .035942547 ? BR5 .07183289 .071832263 ? A76 -.04367290 -.043673546 ? A77 -.05565937 -.055659667 ? A78 -.1831053 -.183105479 ? A79 -.2300432 -.230043850 ? A 2.685216 2.68520989 ? B 2.015685 2.01568802 ? LogL -4948.494 -4948.49442 ? iters 23 9 ? ? NegBin total with R.E. - eq. (3.5) frml leqtre logltre = lgamfn(a+b) + lgamfn(a+sumgam) + lgamfn(b+totaln) - lgamfn(a) - lgamfn(b) - lgamfn(a+b+sumgam+totaln) + sumlog - sumconst ; eqsub leqtre sumgam sumlog sumlam Xbs expZd Zd ; param a,2 b,2; ml(hiter=N,hcov=WN,tol=1e-7) leqtre ; copy @logl ltotre ; copy @ncid ktotre ; title "NegBin Marginal Model with Random Effects" ; ? NegBin marginal with R.E. - eq. (4.3) frml leqmre loglmre = lgamfn(a+b) + lgamfn(a+sumgam) + lgamfn(b+totaln) - lgamfn(a) - lgamfn(b) - lgamfn(a+b+sumgam+totaln) + lgamfn(sumgam+totaln) - lgamfn(sumgam) - lconst ; eqsub leqmre sumgam sumlog sumlam Xbs expZd Zd ; ? year effects not identified in marginal const a76,0 a77,0 a78,0 a79,0; ml(hiter=N,hcov=WN,tol=1e-7) leqmre ; copy @logl lmargre ; copy @ncid kmargre ; ? ? The test below for random vs. fixed effects is probably not ? very good, because the likelihood does not partition exactly ? in negative binomial. It works fine for the Poisson model, ? though. ?set testfe = 2*(lcond+lmargre-ltotre) ; ?set dfe = kcond+kmargre-ktotre ; ?title "Test for Randomness of Firm Effects" ; ?cdf(df=dfe,chisq) testfe ;