options crt ;
name hhgfish "HHG Panel - Poisson Count Models" ;
?
?     ====================== POISSON ===========================
?     This TSP program sets up the likelihood function to
?     estimate the Poisson 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 equation 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 model.  
?
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.
?      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.l0 ;                      ]
  
  frml Xb. bt*year. + br*logr.l0 ;
  frml lamb. exp(Xb.) ;
  eqsub lamb. Xb. ;
  frml logl. n.*Xb. ;
  eqsub logl. Xb. ;
enddot ;

?    Sum of probabilities for denominator 

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 ;

?    Model for non-time-varying variables.

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

?    Likelihood functions.

?    Poisson total - eq. (1.2)

frml leqt loglt = -alpha*sumlam + sumlog + totaln*Zd - sumconst ;
eqsub leqt sumlam alpha sumlog Zd ;

?    Poisson conditional - eq. (2.5)

frml leqc loglc = sumlog - totaln*log(sumlam) - sumconst + lconst ;
eqsub leqc sumlog sumlam ;

?    Poisson marginal - eq. not in paper, but same as totals for sum(n)
?    Note that loglm = loglt-loglc, but it is more efficient to 
?    specify it separately. 

frml leqm loglm = -alpha*sumlam +totaln*(log(sumlam)+Zd) - lconst ;
eqsub leqm alpha sumlam Zd ; 

?    Poisson total with R.E. - eq. (2.3)

frml leqtre logltre = totaln*Zd + sumlog + delta*log(delta)
                   - (delta+totaln)*log(alpha*sumlam+delta) 
                   - lgamfn(delta) + lgamfn(delta+totaln) - sumconst ;
eqsub leqtre Zd sumlog alpha sumlam ;

?    Poisson marginal with R.E. - eq. (4.2)
?    Note that loglmre = logltre-loglc, also. 

frml leqmre loglmre = totaln*(log(sumlam)+Zd) + delta*log(delta)
                   - (delta+totaln)*log(alpha*sumlam+delta) 
                   - lgamfn(delta) + lgamfn(delta+totaln) - lconst ;
eqsub leqmre Zd sumlog alpha sumlam ;

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

title "Poisson 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 "Poisson 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 "Poisson Total Model with Random Effects" ;

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

title "Poisson 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 ;