st: A code for Maximum Simulated Likelihood for an ordered heckman with panel data

[antonio martuscelli <totissoit@yahoo.it>]

Dear all,
I am trying to write a code using the ml command for the following
model: the selection equation is an ordered probit with three classes
and for each class there is a linear outcome equation. The dataset is
a panel so I want to include individual time invariant random effects
in both selection and the three outcome equations. I am assuming for
the moment that these individual effects are uncorrelated with each
other and thus each is distributed as an univariate normal. I am using
MSL to approximate the likelihood function by drawings random numbers
from an univariate normal for each of the four random effects. I use
the d0 method. The code doesn't work and I am failing in identifying
the problem(s). I report the code below and the message I get from
stata. Given that the code works for the simpler pooled case without
random effects I think there is some mistake in the way I am
introducing them or in the simulation but I have tried different
alternatives and none works. Another option could be that the model is
not identified and I need to impose some restrictions. If you have any
suggestions on what I am doing wrong that would be very helpful.

+++++++++++++++++++++++++++++++++++++++++++++++++++code++++++++++++++++++++++++++++++++++++
matrix p = (7, 11, 13, 17)
global draws "100"
by id2: keep if _n==1
mdraws, neq(4) dr($draws) prefix(c) burn(15) prime(p)
forvalues r=1/$draws{
gen random_1`r'=invnormal(c1_`r')
gen random_2`r'=invnormal(c2_`r')
gen random_3`r'=invnormal(c3_`r')
gen random_4`r'=invnormal(c4_`r')
 }
save "C:\mdraws.dta”, replace
 use "C:\simulation_panel_1.dta", clear
merge id2 using "C:\mdraws.dta"

program define mysim_d0
args todo b lnf
tempvar etha1 etha2 etha3 etha4 random1 random2 random3 random4 lj pi1
pi2 pi3 sum lnpi L1 L2 last z1 z2 z3
tempname lnsig1 lnsig2 lnsig3 lnsig4 sigma1 sigma2 sigma3 sigma4 lns1
lns2 lns3 mu0 lndelta athrho1 athrho2/// athrho3 rho1 rho2 rho3 delta
mu1

mleval `etha1' = `b', eq(1)***coefficients
mleval `etha2' = `b', eq(2)
mleval `etha3' = `b', eq(3)
mleval `etha4' = `b', eq(4)
mleval `mu0' = `b', eq(5) scalar*****cut-points
mleval `lndelta' = `b', eq(6) scalar
mleval `lnsig1' = `b', eq(7) scalar****log sd time invariant effects**
mleval `lnsig2' = `b', eq(8) scalar
mleval `lnsig3' = `b', eq(9) scalar
mleval `lnsig4' = `b', eq(10) scalar
mleval `lns1' = `b', eq(11) scalar***log sd errors, restricted to 1
for the oprobit***
mleval `lns2' = `b', eq(12) scalar
mleval `lns3' = `b', eq(13) scalar
mleval `athrho1' = `b', eq(14) scalar**correlation coeff errors*
mleval `athrho2' = `b', eq(15) scalar
mleval `athrho3' = `b', eq(16) scalar

qui {
scalar `sigma1'=(exp(`lnsig1'))^2
scalar `sigma2'=(exp(`lnsig2'))^2
scalar `sigma3'=(exp(`lnsig3'))^2
scalar `sigma4'=(exp(`lnsig4'))^2
gen double `random1' = 0
gen double `random2' = 0
gen double `random3' = 0
gen double `random4' = 0
gen double `lnpi'=0
gen double `sum'=0
gen double `L1'=0
gen double `L2'=0
by id2: gen byte `last'=(_n==_N)
gen double `pi1'=0
gen double `pi2'=0
gen double `pi3'=0
scalar `rho1' = tanh(`athrho1')
scalar `rho2' = tanh(`athrho2')
scalar `rho3' = tanh(`athrho3')
gen double `z1'=0
gen double `z2'=0
gen double `z3'=0
scalar `delta'=exp(`lndelta')
scalar `mu1'=`mu0'+`delta'
}
matrix W = ( `sigma1' , 0 , 0, 0 \ 0, `sigma2' , 0, 0 \ 0, 0, `sigma3'
, 0 \  0, 0, 0, `sigma4')
capture matrix L=cholesky(W)
if _rc != 0 {
di "Warning: cannot do Cholesky factorization of rho matrix"
}
local l11=L[1,1]
local l22=L[2,2]
local l33=L[3,3]
local l44=L[4,4]

forvalues r=1/$draws{
qui {
replace `random1' = random_1`r'*`l11'
replace `random2' = random_2`r'*`l22'
replace `random3' = random_3`r'*`l33'
replace `random4' = random_4`r'*`l44'

*****log likelihood function******
replace `pi1' = ln(normalden(($ML_y1 - `etha1'-`random1') /
exp(`lns1')))-`lns1' if a1==1
replace `z1' = `etha4'+`random4' +`rho1'*($ML_y1 -
`etha1'-`random1')/exp(`lns1') if a1==1
replace `pi1' = `pi1' + ln(normal((`mu0'-`z1')/sqrt(1-`rho1'^2))) if a1==1

replace `pi2' = ln(normalden(($ML_y2 - `etha2'-`random2') /
exp(`lns2')))-`lns2' if a2==1
replace `z2' = `etha4'+`random4' +`rho2'*($ML_y2 -
`etha2'-`random2')/exp(`lns2') if a2==1
replace `pi2' = `pi2' + ln(normal((`z2'-`mu0')/sqrt(1-`rho2'^2))-
normal((`z2'-`mu1')/sqrt(1-`rho2'^2))) if a2==1

replace `pi3' = ln(normalden(($ML_y3 - `etha3'-`random3') /
exp(`lns3')))-`lns3' if a3==1
replace `z3' = `etha4'+`random4' +`rho3'*($ML_y3 -
`etha3'-`random3')/exp(`lns3') if a3==1
replace `pi3' = `pi3' + ln(normal((`z3'-`mu1')/(sqrt(1-`rho3'^2)))) if a3==1

replace `lnpi'=`pi1'*a1+`pi2'*a2+`pi3'*a3

by id2: replace `sum'=sum(`lnpi')
by id2: replace `L1' =exp(`sum'[_N]) if _n==_N
by id2: replace `L2'=`L2'+`L1' if _n==_N
}
}
qui gen `lj'=cond(!`last',0, ln(`L2'/$draws))
qui mlsum `lnf'=`lj'
if (`todo'==0|`lnf'>=.) exit
end
ml model d0 mysim_d0 (y = `x1' `x1bar' if mktpos==0, nocons) (y = `x2'
`x2bar' if mktpos==1, nocons) (y = `x3' `x3bar' if mktpos==2, nocons)
(mktpos = `x4' `x4bar', nocons) /mu0 /lndelta /lnsig1 /lnsig2 /lnsig3
/lnsig4 /lns1 /lns2 /lns3 /athrho1 /athrho2 /athrho3

ml check
ml search
ml maximize

+++++++++++++++++++++++++++++++++++++++++++++++++++++stata message
++++++++++++++++++++++++

.ml search
initial:       log likelihood =         0
rescale:       log likelihood =      0
rescale eq:    log likelihood =     0

.ml maximize

initial:       log likelihood =          0
rescale:       log likelihood =       0
rescale eq:    log likelihood =      0
could not calculate numerical derivatives -- discontinuous region with
missing values encountered
r(430);

*
*   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/