GIBBSORDEREDPROBIT.RPF

GIBBSORDEREDPROBIT.RPF is an example of the use of Gibbs sampling for an ordered probit model. The model estimated is from Example 16.2 of Wooldridge, "Econometrics of Cross Section and Panel Data", 2nd ed.

Gibbs sampling for probit models involves adding parameters that represent the latent "Y" value that govern the choice (here among a set of ordered choices)—in the program that is the series YSTAR.

Full Program

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,nse=bnse) bgibbs

report(action=define)

report(atrow=1,atcol=1,align=center) "Variable" "Coeff" $

"Std Error" "NSE" "CD"

do i=1,%nreg

report(row=new,atcol=1) %eqnreglabels(0)(i) bmean(i) $

bstderrs(i) bnse(i) bcd(i)

end do i

report(action=format,atcol=2,tocol=3,picture="*.###")

report(action=format,atcol=4,picture="*.##")

report(action=show)

Output

Ordered Probit - Estimation by Newton-Raphson

Convergence in 22 Iterations. Final criterion was 0.0000100 <= 0.0000100

Dependent Variable PCTSTCK

Usable Observations 194

Degrees of Freedom 180

Skipped/Missing (from 226) 32

Log Likelihood -201.9865

Average Likelihood 0.3530422

Pseudo-R^2 0.1039465

Log Likelihood(Base) -212.3703

LR Test of Coefficients(14) 20.7676

Significance Level of LR 0.1077377

Variable Coeff Std Error T-Stat Signif

************************************************************************************

1. CHOICE 0.371172029 0.199647148 1.85914 0.06300727

2. AGE -0.050051595 0.021274872 -2.35262 0.01864190

3. EDUC 0.026138084 0.037030562 0.70585 0.48028038

4. FEMALE 0.045564135 0.216217360 0.21073 0.83309561

5. BLACK 0.093391158 0.319116817 0.29266 0.76978581

6. MARRIED 0.093599058 0.235621843 0.39724 0.69118850

7. FINC25 -0.578425724 0.381360200 -1.51674 0.12933138

8. FINC35 -0.134667802 0.399848979 -0.33680 0.73627019

9. FINC50 -0.262035808 0.409480343 -0.63992 0.52222279

10. FINC75 -0.566227297 0.467497878 -1.21119 0.22582379

11. FINC100 -0.227892331 0.452959300 -0.50312 0.61488076

12. FINC101 -0.864109775 0.635012819 -1.36078 0.17358470

13. WEALTH89 -0.000095572 0.000400055 -0.23890 0.81118438

14. PRFTSHR 0.481716027 0.204519202 2.35536 0.01850485

15. Cut(1) -3.087369386 1.567420776 -1.96971 0.04887125

16. Cut(2) -2.053549700 1.564412582 -1.31267 0.18929584

Variable Coeff Std Error NSE CD

CHOICE 0.388 0.183 0.00 3.18

AGE -0.040 0.011 0.00 18.64

EDUC 0.032 0.033 0.00 3.89

FEMALE 0.097 0.184 0.00 6.85

BLACK 0.121 0.280 0.00 -0.09

MARRIED 0.100 0.232 0.00 -0.55

FINC25 -0.526 0.399 0.01 3.05

FINC35 -0.072 0.411 0.01 2.71

FINC50 -0.185 0.393 0.01 3.54

FINC75 -0.483 0.439 0.01 4.79

FINC100 -0.160 0.439 0.01 3.46

FINC101 -0.800 0.499 0.01 2.41

WEALTH89 -0.000 0.000 0.00 -2.97

PRFTSHR 0.503 0.219 0.00 0.92