options double crt;
name poisfe;
?
?  Poisson fixed effects model
?  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  hhgfish.tsp  by Bronwyn H. Hall.
?    by Clint Cummins  12/1998
?
?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 logl. n.*Xb. ;
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 ;

?    Model for non-time-varying variables.

frml Zd b0 + bsize*logk + dsci*scisect ;
param b0 bsize dsci ;
frml alpha exp(Zd) ; 

?    Likelihood functions.

?    Poisson total
title 'Poisson no effects - Table 9.1, columns 1 + 2';

?  Note: it is much easier to reproduce these plain Poisson
?  results by reading the data file as a cross section and
?  using the POISSON command.

?  Corrected results on web page (a couple were rounded wrong):
?          C-T      TSP
? lnR0    .1345    .1345
? lnR-1  -.0529   -.0529
? lnR-2   .0082    .0082
? lnR-3   .0660    .0661
? lnR-4   .0901    .0902
? lnR-5   .2395    .2395
? lnSIZE  .2528    .2529
? DSCI    .4543    .4543
? LogL  -17834.14  -17834.1382
frml leqt loglt = -alpha*sumlam + sumlog + totaln*Zd - sumconst ;
eqsub leqt sumlog sumlam Xbs alpha Zd ;
options nwidth=15,signif=8;
ml(hiter=N,hcov=WN,tol=1e-7) leqt ;     

?    Poisson conditional
title "Poisson Fixed Effects" ;
?   C-T numbers below are from Limdep output file.  The first 3 numbers
?  in the printed table are OK (2 digits), but the last 3 are wrong.
?          C-T       TSP
? lnR0    .3222     .3222
? lnR-1  -.08713   -.08713
? lnR-2   .007858   .007858
? lnR-3   .01060    .01060
? lnR-4  -.004642  -.004641
? lnR-5   .002607   .002607
? LogL  -3536.309  -3536.30861
frml leqc loglc = sumlog - totaln*log(sumlam) - sumconst + lconst ;
eqsub leqc sumlog sumlam Xbs ;
ml(hiter=N,hcov=WN,tol=1e-7) leqc ;