"Fixed Effects" Logit

[TomDoan]

This does a "fixed effects" logit in a panel data set. (See Wooldridge (2002), section 15.8.3 for a description). This particular data set is unbalanced, with the individuals identified by the <<idcode>> series. With a balanced panel, set up with CAL with the PANELOBS option, you could replace the SSTATS that computes nfloat with

Code: Select all

compute n=NUMBEROFINDIVIDUALS

and replace the do i=1,26200 loop with

Code: Select all

do i=1,n
   compute bounds(1,i)=1+(i-1)*TIMEDIMENSION
   compute bounds(2,i)=i*TIMEDIMENSION
   sstats bounds(1,i) bounds(2,i) DEPENDENTVARIABLE>>count(i)
end do i

NUMBEROFINDIVIDUALS, TIMEDIMENSION and DEPENDENTVARIABLE are whatever is appropriate for your data.

There's also one occurrence of <<union>>, which is the dependent variable in this example, in the FELogitElement function, which you'll have to replace with your dependent variable.

Code: Select all

*
* "Fixed Effects" logit model
*
open data union.dat
data(format=prn,org=columns) 1 26200 idcode year age grade not_smsa south union t0 southxt black
*
* Data set is unbalanced. This next line counts the number of distinct
* idcodes (assuming that an individual's data is in consecutive entries)
*
sstats 1 26200 (idcode.ne.idcode{1})>>nfloat
compute n=fix(nfloat)
*
* This locates the boundaries for each individual since the likelihood
* function has to be integrated across all data points for an individual.
*
dec rect[int] bounds(n,2)
dec vect[int] count(n)
compute lastcode=0.0,fill=0
do i=1,26200
   if idcode(i)==lastcode {
      compute bounds(fill,2)=i
      compute count(fill)=count(fill)+fix(union(i)==1)
   }
   else {
      compute fill=fill+1
      compute bounds(fill,1)=i
      compute bounds(fill,2)=i
      compute count(fill)=fix(union(i)==1)
      compute lastcode=idcode(i)
   }
end do i
*
* This does the standard logit
*
ddv(dist=logit) union
# age grade not_smsa south southxt constant
*
* Set up the formula. The constant can't be estimated
*
frml(regressors,vec=b) zfrml
# age grade not_smsa south southxt
nonlin b
*
* This computes the partial log likelihood for individual "t". This is the
* relative probability of the individual having a particular set of "1"'s among
* all possible ways of having that many.
*

function FELogitElement t
type real FELogitElement
type integer t
*
local vector zvalues
local integer nt nones i j done
local vect[integer] ones
local real sumcombos sumz actual
*
compute nt=bounds(t,2)-bounds(t,1)+1
dim zvalues(nt)
*
* If all values are 0 or 1, the ML value is to make the individual effect either
* -infinity or +infinity, so there is no information about the slope
* coefficients.
*
compute nOnes=count(t)
if nOnes==0.or.nOnes==nt {
   compute FELogitElement=%na
   return
}
*
* Compute the index values given the current setting for the coefficients
*
ewise zvalues(i)=zfrml(bounds(t,1)+i-1)
*
* Set up the list of "ones"
*
dim ones(nOnes)
ewise ones(i)=i
*
compute sumcombos=0.0
compute done=0
while .not.done {
   compute sumz=0.0
   do i=1,nOnes
      compute sumz=sumz+zvalues(ones(i))
   end do i
   compute sumcombos=sumcombos+exp(sumz)
   *
   * Update the set of ones
   *
   compute done=1
   do i=nOnes,1,-1
      compute ones(i)=ones(i)+1
      if ones(i)>nt-nOnes+i
         next
      do j=i+1,nOnes
         compute ones(j)=ones(j-1)+1
      end do j
      compute done=0
      break
   end do i
}
compute actual=0.0
do j=bounds(t,1),bounds(t,2)
   if union(j)==1
      compute actual=actual+zvalues(j-bounds(t,1)+1)
end do j
compute FELogitElement=exp(actual)/sumcombos
end
*
* The maximize runs over count of individuals, not the actual entries.
*
frml felogit = log(FELogitElement(t))
maximize(method=bfgs,title="Fixed Effects Logit",trace) felogit 1 n