Dear Stata users,
I am trying to estimate a hurdle model using ml. The first part of the
likelihood function is a Probit, the second part is a Weibull duration model
(PH) with right censoring.
I think my problems are due to the duration and the censoring variables: the
variable _MLw1 created by Stata during the maximisation has missing values,
corresponding to the values for which the duration and the censoring
variable are not observed (and the binary choice variable in the probit
takes value 0).
The hurdle models code in Stata, those that I know, are for counts and
require a 0 value in the count variable. So I am replacing missings values
in the time variable with 0.
capture program drop hurdleSQ
program define hurdleSQ
version 9.2
args lnf xb1 xb2 alfa
tempvar normp lambda df sf ldf lsf weib k1
local s "$ML_y1" /*choice variable; if s=1 then t is observed,
otherwise t is missing*/
local t "$ML_y2" /*time variable*/
local q "$ML_y3" /*censoring variable*/
gen byte `k1'= 2*`s' - 1
quietly{
gen double `normp'= ln(norm(`k1'*`xb1'))
gen double `lambda1' = exp(`xb2')
gen double `df' =
`alfa'*`t'^(`alfa'-1)*`lambda'*exp(-`lambda'*`t'^(`alfa'))
gen double `sf' = exp(-`lambda'*`t'^(`alfa'))
gen double `ldf'= ln(`cdf')
gen double `lsf'= ln(`sf')
gen double `weib'= (`q'*`lcdf'+(1-`q')*`lsf' )
replace `lnf'= cond(`t'==0,`normp',`weib')
}
end
replace smy=0 if smy==.
ml model lf hurdleSQ (probit: s = $rhs1) (duration: t q = $rhs2 ) /alfa
ml init initmx, skip copy /*initmx are initial values from the separate
probit and weibull model*/
ml max, difficult
and then I get...
initial: log likelihood = -<inf> (could not be evaluated)
feasible: log likelihood = -4034.2574
rescale: log likelihood = -1785.9959
rescale eq: log likelihood = -1438.3167
could not calculate numerical derivatives
flat or discontinuous region encountered
I don't think that replacing missing values in the censoring variable can be
of any help. Probably there is something to do inside the program. I have
tried to use the option "missing" after ml model, but still I get the same
error message.
Do you see a way to solve this problem or do you spot any mistake in the
likelihood?
Thanks in advance,
Silvia
*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/
