*
* REPROBIT.RPF
* Random Effects probit 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<>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)
compute lastcode=0.0,fill=0
do i=1,26200
   if idcode(i)==lastcode
      compute bounds(fill,2)=i
   else
      compute fill=fill+1,bounds(fill,1)=i,lastcode=idcode(i)
end do i
*
* This does the standard probit
*
ddv union
# age grade not_smsa south southxt constant
*
* Set up the formula. The standard error of the RE component also needs
* to be estimated.
*
frml(lastreg,vec=b) zfrml
nonlin b sigma
compute sigma=1.0
*
* These are the coefficients for a five term Hermite integral
*
dec vect w(5)
dec vect x(5)
input w
 .9453087204829
 .3936193231522
 .3936193231522
 .01995324205905
 .01995324205905
input x
   0
  .958572464613819
 -.958572464613819
 2.020182870456086
-2.020182870456086
*
* This does the integral over the data for individual "t"
*
function REIntegral t
type real REIntegral
type integer t
*
local integer i j
local real prod z
*
compute REIntegral=0.0
do i=1,5
   compute prod=1.0
   do j=bounds(t,1),bounds(t,2)
      compute z=zfrml(j)+sqrt(2)*sigma*x(i)
      compute prod=prod*%if(union(j)==1,%cdf(z),%cdf(-z))
   end do j
   compute REIntegral=REIntegral+w(i)*prod
end do i
return
end REIntegral
*
* The maximize runs over count of individuals, not the actual entries.
*
frml reprobit = log(REIntegral(t))
maximize(method=bfgs,title="Random Effect Probit") reprobit 1 n