The RATS Software Forum

Code: Select all

*
* Example of Gibbs sampling for an ordered probit model.
*
* Based upon Wooldridge, "Econometrics of Cross Section and Panel Data",
* 2nd ed. Example 16.2 from pp 657-658
*
open data pension.raw
data(format=free,org=columns) 1 226 id pyears prftshr choice female $
  married age educ finc25 finc35 finc50 finc75 finc100 finc101 $
  wealth89 black stckin89 irain89 pctstck
*
* Linear model.
*
linreg pctstck
# choice age educ female black married finc25 finc35 finc50 finc75 $
  finc100 finc101 wealth89 prftshr constant
*
* Ordered probit estimated by ML.
*
ddv(type=ordered,dist=probit,cuts=cuts) pctstck
# choice age educ female black married finc25 finc35 finc50 finc75 $
  finc100 finc101 wealth89 prftshr
*
compute gstart=%regstart(),gend=%regend()
*
* Flat prior on coefficients
*
dec symm hprior
dec vect bprior
compute hprior=%zeros(%nreg,%nreg)
compute bprior=%zeros(%nreg,1)
*
* YSTAR will be the latent variable which determines the choice. It's
* treated as a set of parameters in the sampler.
*
dec series ystar
*
* Create an equation with the regressors from the ML estimate of the
* ordered probit. Save the estimates into BDRAW, which will be used for
* the draws for the coefficients.
*
equation(lastreg) eqn ystar
compute bdraw=%beta
*
* Get the set of values of the dependent variable.
*
@uniquevalues(values=yvalues) pctstck
compute ngroup=%size(yvalues)
*
* For convenience, create a SERIES[INTEGER] which maps an entry to the
* position in the YVALUES vector based upon the dependent variable.
*
dec series[int] which
gset which = 0
do time=gstart,gend
   do i=1,ngroup
      if pctstck(time)==yvalues(i) {
         compute which(time)=i
         break
      }
   end do i
end do time
*
* This will be used inside the draw loop
*
dec vect lcat(ngroup) hcat(ngroup)
*
* These will be the parameters that we save on each draw: the regression
* parameters for the index function and the cut points.
*
nonlin(parmset=oprobit) bdraw cuts
*
compute ndraws=10000
compute nburn=10000
*
dec series[vect] bgibbs
*
gset bgibbs 1 ndraws = %parmspeek(oprobit)
*
infobox(action=define,progress,lower=-nburn,upper=ndraws) $
  "Ordered Probit Model-Gibbs Sampling"
do draw=-nburn,ndraws
   *
   * Draw y*'s given regression coefficients and cuts. This is a
   * truncated Normal with mean equal to X beta at the current set of
   * coefficients with lower and upper bounds determined by the cut
   * points.
   *
   set cmean gstart gend = %eqnvalue(eqn,t,bdraw)
   set lower gstart gend = %if(which==1,%na,cuts(which-1))
   set upper gstart gend = %if(which==ngroup,%na,cuts(which))
   *
   set ystar gstart gend = %rantruncate(cmean,1.0,lower,upper)
   *
   * Draw beta's given y*'s. Given the y*'s, this is just a
   * straightforward draw for a regression (with residual variance 1).
   *
   cmom(equation=eqn)
   compute bdraw=%ranmvpostcmom(%cmom,1.0,hprior,bprior)
   *
   * Draw cut points. These are drawn uniformly from the values in the
   * gaps between the y*'s for the different values of "which".
   *
   do i=1,ngroup
      ext(smpl=(which==i),noprint) ystar gstart gend
      compute lcat(i)=%minimum,hcat(i)=%maximum
   end do i
   *
   do i=1,ngroup-1
      compute cuts(i)=%uniform(hcat(i),lcat(i+1))
   end do i
   *
   infobox(current=draw)
   *
   * If we're out of the burn-in, save the parameters
   *
   if draw>0
      compute bgibbs(draw)=%parmspeek(oprobit)
end do draw
infobox(action=remove)
*
@mcmcpostproc(ndraws=ndraws,mean=bmean,stderrs=bstderrs,cd=bcd) bgibbs