Does any one have any idea about how to estimate a sequential ordered
probit model using Stata?
I am aware of bioprpobit program written by Zurab Sajaia wihich allows
for simultaneous specification. I am trying to modify this program to
estimate a sequential ordered probit model. In my model the first
categorical variable is observed for the entire sample. The values it
takes are 0,1,2,3,4 and 5.
The second categorical variable is observed only for those individuals
for whom the first variable is strictly greater than zero. The second
categorical variable takes only two values, 0 and 1.
For starting values, I used the estimates from two regressions, an
ordered probit regression for the first ordered variable and a probit
regression for the second variable. Despite providing the initial
values, my ml program keeps saying unable to find feasible values.
I would highly appreciate if you could provide some hints as to how to
write a code for implementing a sequential ordered probit model using
stata and overcome this problem of feasible values.
Thanks a ton,
Sridhar Telidevara
Link to Zurab Sajaia's bioprobit:
http://fmwww.bc.edu/repec/bocode/b/bioprobit_d2.ado
Here is the modified program:
capture program drop bi_prob_d0
version 11
use stata_1.dta,replace
program define bi_prob_d0, eclass
quietly{
forvalues i = 1/10 {
local grad "`grad' g`i'"
}
args todo b lnf g negH `grad'
tempvar p1_b p2_b qq c12
tempname r rho gamma d dtanh_dr drho_dr drho_dgamma rr drr1_dr d2tanh_dr2
mleval `p1_b' = `b' , eq(1)
mleval `p2_b' = `b' , eq(2)
mleval `r' = `b' , eq(3) scalar
mleval `gamma' = `b' , eq(4) scalar
forvalues k = 1/5 {
tempvar c`k'1
mleval `c`k'1' = `b' , eq(`=4+`k'') scalar
}
mleval `c12' = `b' , eq(10) scalar
scalar `rr' = 1/sqrt(1+(`gamma')^2+2*`gamma'*tanh(`r'))
scalar `rho' = (tanh(`r')+`gamma')*`rr'
tempvar cut1 cut1_1 cut2 cut2_1
tempname CUT1 CUT2
local neginf = minfloat()
local posinf = maxfloat()
matrix `CUT1' = `neginf' \ `c11' \ J(4,1,.) \ `posinf'
matrix `CUT2' = `neginf' \ `c12' \ `posinf'
forvalues k = 2/5 {
matrix `CUT1'[`k'+1,1] = `CUT1'[`k',1]+`c`k'1'^2
}
generate double `cut1' = `CUT1'[$ML_y1+1,1]
generate double `cut1_1' = `CUT1'[$ML_y1 ,1]
generate double `cut2' = `CUT2'[$ML_y2+1,1]
generate double `cut2_1' = `CUT2'[$ML_y2 ,1]
local arg11 "(`cut1' -`p1_b')"
local arg21 "(`rr'*(`cut2' -`p2_b'-`gamma'*`p1_b'))"
local cond11 "$ML_y1"
local cond21 "$ML_y2"
local arg12 "(`cut1_1'-`p1_b')"
local arg22 "(`rr'*(`cut2' -`p2_b'-`gamma'*`p1_b'))"
local cond12 "$ML_y1-1"
local cond22 "$ML_y2"
local arg13 "(`cut1' -`p1_b')"
local arg23 "(`rr'*(`cut2_1'-`p2_b'-`gamma'*`p1_b'))"
local cond13 "$ML_y1"
local cond23 "$ML_y2-1"
local arg14 "(`cut1_1'-`p1_b')"
local arg24 "(`rr'*(`cut2_1'-`p2_b'-`gamma'*`p1_b'))"
local cond14 "$ML_y1-1"
local cond24 "$ML_y2-1"
forvalues i = 1/4 {
tempvar q`i'
generate double `q`i''=0
replace `q`i''= binormal(`arg1`i'',`arg2`i'',`rho') if ($ML_y1>0 )
}
replace `q1' = normal(`arg11') if $ML_y1==0
generate double `qq' = `q1'-`q2'-`q3'+`q4'
replace `qq' = 0.1D-60 if(`qq'<=0)
mlsum `lnf' = ln(`qq')
if (`todo'==0 | `lnf'==.) exit
// calculate first derivatives
scalar `d' = 1/sqrt(1-(`rho')^2)
scalar `dtanh_dr' = 4*exp(2*`r')/((1+exp(2*`r'))^2)
scalar `drho_dr' = (`rr'-`gamma'*`rho'*(`rr')^2)*`dtanh_dr'
scalar `drr1_dr' = -`gamma'*(`rr')^2*`dtanh_dr'
scalar `drho_dgamma' = `rr'*(1-(`rho')^2)
forvalues i = 1/4 {
tempvar dF1`i' dF2`i' dF3`i' d1q`i' d3q`i' d10q`i'
gen double `dF1`i''=0
gen double `dF2`i''=0
gen double `dF3`i''=0
replace `dF1`i'' = normalden(`arg1`i'')*normal(`d'*(`arg2`i'' -
`rho'*`arg1`i'')) if $ML_y1>0
replace `dF2`i'' = normalden(`arg2`i'')*normal(`d'*(`arg1`i'' -
`rho'*`arg2`i'')) if $ML_y1>0
replace `dF3`i'' =
exp(-0.5*((`arg1`i'')^2+(`arg2`i'')^2-2*`rho'*`arg1`i''*`arg2`i'')*`d'^2)/(2*_pi*sqrt(1-(`rho')^2))
if $ML_y1>0
replace `dF11' = normalden(`arg11') if $ML_y1 ==0
gen double `d1q`i'' = - `dF1`i'' - `gamma'*`rr'*`dF2`i''
local `d2q`i'' = "(-`dF2`i''*`rr')"
gen double `d3q`i'' = `dF2`i''*`arg2`i''*`drr1_dr'+`dF3`i''*`drho_dr'
gen double `d4q`i'' =
-(`rho'*`arg2`i''+`p1_b')*`rr'*`dF2`i''+`dF3`i''*`drho_dgamma'
local `d5q`i'' = "`dF1`i''"
generate double `d10q`i'' = `rr'*`dF2`i''
forvalues k = 2/5 {
local eqn = 4+`k'
tempvar d`eqn'q`i'
generate double `d`eqn'q`i'' = cond(`cond1`i''>=`k',
`dF1`i''*2*`c`k'1', 0)
}
}
capture matrix drop `g'
forvalues i = 1/10{
tempvar g`i'i
generate double `g`i'i' = `d`i'q1'-`d`i'q2'-`d`i'q3'+`d`i'q4'
replace `g`i'' = `g`i'i'/`qq'
mlvecsum `lnf' `g`i'' = `g`i'', eq(`i')
matrix `g' = (nullmat(`g'), `g`i'')
}
}
end
*
* For searches and help try:
* http://www.stata.com/help.cgi?search
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
