options crt ;
name hhgnb "HHG Panel - Negative Binomial Count Models" ;
?
?     ====================== NEG. BIN. ===========================
?     This TSP program sets up the likelihood function to estimate
?     the Negative Binomial models for panel count data described 
?     in Hausman, Hall, and Griliches (Econometrica 1984). The 
?     equation numbers refer to that paper.
?
?     The example has 7 years of data. Note that all the data
?     for a single individual or firm are read into a single
?     observation (unlike the LOGIT procedure with COND option). 
?     The XB and Zd equations can easily be rewritten to add or 
?     subtract variables. The parameter for the variance of  
?     the random effects is delta.
?
?     TSP 4.3 (Sept. 95) 
?     Author:  Bronwyn H. Hall
?
smpl 1 57 ;
read (file="hhg84.dat") 
      cusip pat73 logr73l0-logr73l5 year73 logk scisect 
      skip pat74 logr74l0-logr74l5 year74 skip skip
      skip pat75 logr75l0-logr75l5 year75 skip skip
      skip pat76 logr76l0-logr76l5 year76 skip skip
      skip pat77 logr77l0-logr77l5 year77 skip skip
      skip pat78 logr78l0-logr78l5 year78 skip skip
      skip pat79 logr79l0-logr79l5 year79 skip skip ;
?
?   This is the data from Hall, Griliches, and Hausman (IER 86).
?   Only 399 obs = 57 firms x 7 years are used.
?
?   The variables are the following:
?
?   cusip        firm identifier.
?   patents      no. of successful patent applications in the year.
?   logr0        log of R&D spending in the year.
?   logr1-logr5  5 lags of log of R&D spending.
?   year         year of data.
?   logk         log of R&D stock at beginning of year.
?   scisect      dummy for firms in scientific sector (see HGH).
?
msd pat73-pat79 logr73l0-logr73l5 logk scisect ;

?    Dependent variable is called n; compute constants for
?    likelihoods.

totaln = 0 ;
dot 73-79 ;
  n.=pat. ;
  totaln = totaln+n. ; 
  lconst. = lgamfn(n.+1) ;
enddot ;
sumconst = lconst73+lconst74+lconst75+lconst76+lconst77+lconst78
         +lconst79 ;
lconst = lgamfn(totaln+1) ;
msd totaln ;
const t 7 ;                ? Number of years per firm
?
?   Construct the likelihood function for each of 5 models.  
?
?      Model for non-time-varying variables.
?      Zd contains terms that do not vary over time (factor out of
?      the lambda term).  

frml Zd b0 + bk*logk + bs*scisect ;
frml alpha exp(Zd) ; 
eqsub alpha Zd ;

dot 73-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.
?      The full model for the expected value 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.l0 ;                      ]
  
  frml Xb. bt*year. + br*logr.l0 ;
  frml lamb. exp(Xb.) ;
  eqsub lamb. Xb. ;
  frml gam. alpha*lamb. ;
  eqsub gam. alpha lamb. ;
  frml logl. lgamfn(gam. + n.) - lgamfn(gam.) ;
  eqsub logl. gam. ;
  frml loglc. lgamfn(lamb. + n.) - lgamfn(lamb.) ;
  eqsub loglc. lamb. ;
enddot ;

?    Sums of Terms. 

frml sumlam lamb73+lamb74+lamb75+lamb76+lamb77+lamb78+lamb79 ;
eqsub sumlam lamb73-lamb79 ;
frml sumlog logl73+logl74+logl75+logl76+logl77+logl78+logl79 ;
eqsub sumlog logl73-logl79 ;
frml sumlogc loglc73+loglc74+loglc75+loglc76+loglc77+loglc78
            +loglc79 ;
eqsub sumlogc loglc73-loglc79 ;

frml sumgam alpha*sumlam ;
eqsub sumgam sumlam alpha ;

?    Likelihood functions.

?    Neg. Bin. total - eq. (3.1) 

frml leqt loglt = sumlog - (sumgam+totaln)*log(1+delta) 
                  + sumgam*log(delta) - sumconst ;
eqsub leqt sumgam sumlog ;

?    Neg. Bin. conditional - eq. (3.2) 

frml leqc loglc = sumlogc + lgamfn(sumlam) - lgamfn(sumlam+totaln) 
                - sumconst + lconst ;
eqsub leqc sumlog sumlam sumlogc ;

?    Neg. Bin. marginal - eq. (not in paper, same as 3.1 for the sum) 
?    Note that loglm is not exactly equal to loglt-loglc; full 
?    partitioning is not possible in the negative binomial case.  

frml leqm loglm = -lgamfn(alpha*sumlam) + lgamfn(alpha*sumlam+totaln)  
                - (sumgam+totaln)*log(1+delta) + sumgam*log(delta) 
                - lconst ;
eqsub leqm alpha sumlam sumgam ; 

?    Neg. Bin. 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 sumlog sumgam ;

?    Poisson 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 ;

title "Negative Binomial Totals (No Effects)" ;
           
param br bt delta .05 ;
param b0 bk bs ;
ml (maxit=40,tol=.01,hcov=NBW,hiter=N) leqt ;     
copy @logl ltotal ; copy @ncid ktotal ;

title "Negative Binomial Marginal (No Effects)" ;
           
param br delta .05 ;
param b0 bk bs ;
const bt 0 ;              ? trend is not identified in marginal.
ml (maxit=40,tol=.01,hcov=NBW,hiter=N) leqm ;     
copy @logl lmarg ; copy @ncid kmarg ;

title "Negative Binomial Conditional" ;

param bt ;
ml (maxit=40,tol=.01,hcov=NBW,hiter=N) 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 "Negative Binomial Total Model with Random Effects" ;

param a 2 b 2 ;
ml (maxit=40,tol=.01,hcov=NBW,hiter=N) leqtre ;
copy @logl ltotre ; copy @ncid ktotre ;

title "Negative Binomial Marginal Model with Random Effects" ;

const bt 0 ;               ? trend is not identified in marginal
ml (maxit=40,tol=.01,hcov=NBW,hiter=N) leqmre ;   
copy @logl lmargre ; copy @ncid kmargre ;
set testfe = 2*(lcond+lmargre-ltotre) ;
set dfe = kcond+kmargre-ktotre ;
title "Test for Randomness of Firm Effects" ;
cdf(df=dfe,chisq) testfe ;

stop ; end ;