Error message with a log-likelihood program

[Olivier Supplisson]

Dear Stata community,

I'm trying to estimate a semi-structural model using log-likelihood estimator. Given my model, I have to code my own program.

According to stata coma "ml check" my program has an accurate syntax and, a priori, my (theoretical) log-likelihood is ok. However, so far, when I've launched my program stata have sent the following error message : "could not calculate numerical derivatives -- flat or discontinuous region encountered" (with a simple or more complete model). (I have the same message using any other altorithm and the "difficult" option).

I'd like to know if my coded log-likelihood matches the theoretical one (e.g. that I coded what I wanted to code or if I made a mistake) and if that error message comes from my data and not from my program.


*Theoretical LL :

Ln(L)=
\sum_{h_{i}<\bar{H}_{i}}  ln(S_{1i})
+\sum_{h_{i}>\bar{H}_{i}}  ln(S_{2i})
+\sum_{h_{i}=\bar{H}_{i}}  ln(S_{3i})

With :

S_{1i}=\frac{1}{\sigma} \phi(\frac{h_{i}-x'_{i}\beta-z'\psi}{\sigma})
S_{2i}=\frac{1}{\sigma} \phi(\frac{h_{i}-x'_{i}\gamma-z'\psi}{\sigma})
S_{3i}=\Phi(\frac{\bar{h]_{i}-x'_{i}\gamma-z'\psi}{\sigma})-\Phi(\frac{\bar{h]_{i}-x'_{i}\beta-z'\psi}{\sigma})


Where \phi stands for the normal (0,1) density and \Phi the normal (0,1) cumulative distribution function

*Stata code :

program define ModeleNormatifEchantillonTotal
version 1.0 args lnf xbeta xgamma sigma
tempvar f1 f2 F1 F2

quietly {
*Définition des parties de la ML
*P1
gen double `f1' = (1/(`sigma'))*normalden((1/(`sigma'))*($ML_y1-`xbeta'))
*P2
gen double `f2' = (1/(`sigma'))*normalden((1/(`sigma'))*($ML_y2-`xgamma'))
*P3
gen double `F1' = normal(($ML_y3-`xbeta')/`sigma')
*P4
gen double `F2' = normal(($ML_y4-`xgamma')/`sigma')

* Fonction de vraisemblance selon les 3 cas possibles
replace `lnf' = ln(`f1') if $ML_y1 < $ML_y3
replace `lnf' = ln(`f2') if $ML_y1 > $ML_y3
replace `lnf' = ln(`F2'-`F1') if $ML_y1 == $ML_y3

}
end

constraint drop _all
constraint 1 [P1]e=[P2]e
constraint 2 [P1]e=[P3]e
constraint 3 [P1]a=[P3]a
constraint 4 [P2]a=[P4]a
constraint 5 [P1]a=[P3]a
constraint 6 [P2]b=[P4]b
constraint 7 [P1]b=[P3]b
constraint 8 [P2]b=[P4]b

noi xi:ml model lf ModeleNormatifEchantillonTotal (P1: Y=a b e) (P2:Y=a b c e) (P3:Z=a b e)(P4:Z=a b c e)(sigma: ), constraint (1-8)
ml check
ml search
noi ml maximize, difficult

Do I need to put $ML_y3 instead of $ML_y4 in F2 ? etc.

I hope my message was clear.