Hi Sridhar,
First of all, I'm sorry for not being able to provide any specific comments on your code, I usually have hard time reading and understanding even my own ML programs (including -bioprobit- that is :$ ) with all the 'smart' ways to loop over N equations, cut-off points, excessive use of macro substitution (powerful tool but very difficult to read), modifying even slightly an ML program might become an nightmare of continuous debugging ... but here are few general ideas:
If you want to learn how to write ML estimators in Stata, the very first thing to do would be buying this book: http://www.stata.com/bookstore/mle.html
The model you're describing, although sounds similar as bioprobit, will have different likelihood function, so I'd suggest first starting with writing it down. You have to program the (log) likelihood first and see how it performs on artificial data (you can plan your DGP and equations and see hew well the parameters are recovered). First and second derivatives will come much later, so for starts it might be easier to work with (Stata's official) -biprobit- code and modify it for your model. Since I have arbitrary number of cut-offs for both equations in the -bioprobit-, plus to make derivative calculations easier I create some constructs in the beginning that perhaps complicate reading the ll calculation part.
So I'd start from -biprobit-, modify the ll to accommodate the sequential model, and support only fixed number of cutoffs (6 as in your example) for now, and see how it goes. Then once you see it works, start adding derivatives etc.
Hope this helps,
zurab
----------------------------------------
> Date: Thu, 21 Jan 2010 10:23:54 -0500
> Subject: st: Sequential ordered probit model
> From: [email protected]
> To: [email protected]
>
> 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/
*
* 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/
