* * 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 * * * Get coefficients for a 5 term Gauss-Hermite integral. More terms makes * this more accurate, but considerably slows the estimation. * compute ngh=5 @GaussHermite(n=ngh,w=w,x=x,standard) * * 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,ngh compute prod=1.0 do j=bounds(t,1),bounds(t,2) compute z=zfrml(j)+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