Dear all,
I am trying to program a choice model, in three stages with two types
in the population. The first is a binonial logit, the second a binomial
logit and the third a multinomial logit with 4 alternatives.
I have programed the maximum likelihood
function, but I have problems with convergence. I keep recieving a message
saying the function is not concave in all iterations. I am not sure if
something is wrong with my program (this the first one I have ever
wrote) or if the problem is that I need to set some appropiate initial
values. Any suggestions?
The second problem is that one of the parameters should be between zero
and one. Is there a way to constraint the domain, over which I wil
maximize?
Thanks,
Larissa
------------------------
program drop likefn
program define likefn
version 6
args lnf theta1 theta2 theta3 theta4 theta5 theta6 theta7 theta8 theta9
theta10
tempvar d1d0 d1d1 d1c0 d1c1 d2d0 d2d1 d2c0 d2c1 d3h0 d3h1 d3g0 d3g1 d3p0
d3p1 d3d0 d3d1
quietly gen double `d1d0'=1/(1+exp(`theta1'))
quietly gen double `d1d1'=1/(1+exp(`theta1'+1))
quietly gen double `d1c0'=exp(`theta1')/(1+exp(`theta1'))
quietly gen double `d1c1'=exp(`theta1'+1)/(1+exp(`theta1'+1))
quietly gen double `d2d0'=1/(1+exp(`theta2'))
quietly gen double `d2d1'=1/(1+exp(`theta2'+`theta6'))
quietly gen double `d2c0'=exp(`theta2')/(1+exp(`theta2'))
quietly gen double
`d2c1'=exp(`theta2'+`theta6')/(1+exp(`theta2'+`theta6'))
quietly gen double
`d3h0'=exp(`theta3')/(1+exp(`theta3')+exp(`theta4')+exp(`theta5'))
quietly gen double
`d3h1'=exp(`theta3'+`theta7')/(1+exp(`theta3'+`theta7')+exp(`theta4'+`theta8')+exp(`theta5'+`theta9'))
quietly gen double
`d3g0'=exp(`theta4')/(1+exp(`theta3')+exp(`theta4')+exp(`theta5'))
quietly gen double
`d3g1'=exp(`theta4'+`theta8')/(1+exp(`theta3'+`theta7')+exp(`theta4'+`theta8')+exp(`theta5'+`theta9'))
quietly gen double
`d3p0'=exp(`theta5')/(1+exp(`theta3')+exp(`theta4')+exp(`theta5'))
quietly gen double
`d3p1'=exp(`theta5'+`theta9')/(1+exp(`theta3'+`theta7')+exp(`theta4'+`theta8')+exp(`theta5'+`theta9'))
quietly gen double `d3d0'=1/(1+exp(`theta3')+exp(`theta4')+exp(`theta5'))
quietly gen double
`d3d1'=1/(1+exp(`theta3'+`theta7')+exp(`theta4'+`theta8')+exp(`theta5'+`theta9'))
quietly replace
`lnf'=ln((`d1d0'*`d2d0'*`d3h0'*(`theta10'))+(`d1d1'*`d2d1'*`d3h1'*(1-`theta10'))) if
$ML_y1==0 & $ML_y2==0 & $ML_y3==3
quietly replace
`lnf'=ln((`d1d0'*`d2d0'*`d3g0'*(`theta10'))+(`d1d1'*`d2d1'*`d3g1'*(1-`theta10'))) if
$ML_y1==0 & $ML_y2==0 & $ML_y3==2
quietly replace
`lnf'=ln((`d1d0'*`d2d0'*`d3p0'*(`theta10'))+(`d1d1'*`d2d1'*`d3p1'*(1-`theta10'))) if
$ML_y1==0 & $ML_y2==0 & $ML_y3==1
quietly replace
`lnf'=ln((`d1d0'*`d2d0'*`d3d0'*(`theta10'))+(`d1d1'*`d2d1'*`d3d1'*(1-`theta10'))) if
$ML_y1==0 & $ML_y2==0 & $ML_y3==0
quietly replace
`lnf'=ln((`d1c0'*`d2d0'*`d3h0'*(`theta10'))+(`d1c1'*`d2d1'*`d3h1'*(1-`theta10'))) if
$ML_y1==1 & $ML_y2==0 & $ML_y3==3
quietly replace
`lnf'=ln((`d1c0'*`d2d0'*`d3g0'*(`theta10'))+(`d1c1'*`d2d1'*`d3g1'*(1-`theta10'))) if
$ML_y1==1 & $ML_y2==0 & $ML_y3==2
quietly replace
`lnf'=ln((`d1c0'*`d2d0'*`d3p0'*(`theta10'))+(`d1c1'*`d2d1'*`d3p1'*(1-`theta10'))) if
$ML_y1==1 & $ML_y2==0 & $ML_y3==1
quietly replace
`lnf'=ln((`d1c0'*`d2d0'*`d3d0'*(`theta10'))+(`d1c1'*`d2d1'*`d3d1'*(1-`theta10'))) if
$ML_y1==1 & $ML_y2==0 & $ML_y3==0
quietly replace
`lnf'=ln((`d1d0'*`d2c0'*`d3h0'*(`theta10'))+(`d1d1'*`d2c1'*`d3h1'*(1-`theta10'))) if
$ML_y1==0 & $ML_y2==1 & $ML_y3==3
quietly replace
`lnf'=ln((`d1d0'*`d2c0'*`d3g0'*(`theta10'))+(`d1d1'*`d2c1'*`d3g1'*(1-`theta10'))) if
$ML_y1==0 & $ML_y2==1 & $ML_y3==2
quietly replace
`lnf'=ln((`d1d0'*`d2c0'*`d3p0'*(`theta10'))+(`d1d1'*`d2c1'*`d3p1'*(1-`theta10'))) if
$ML_y1==0 & $ML_y2==1 & $ML_y3==1
quietly replace
`lnf'=ln((`d1d0'*`d2c0'*`d3d0'*(`theta10'))+(`d1d1'*`d2c1'*`d3d1'*(1-`theta10'))) if
$ML_y1==0 & $ML_y2==1 & $ML_y3==0
quietly replace
`lnf'=ln((`d1c0'*`d2c0'*`d3h0'*(`theta10'))+(`d1c1'*`d2c1'*`d3h1'*(1-`theta10'))) if
$ML_y1==1 & $ML_y2==1 & $ML_y3==3
quietly replace
`lnf'=ln((`d1c0'*`d2c0'*`d3g0'*(`theta10'))+(`d1c1'*`d2c1'*`d3g1'*(1-`theta10'))) if
$ML_y1==1 & $ML_y2==1 & $ML_y3==2
quietly replace
`lnf'=ln((`d1c0'*`d2c0'*`d3p0'*(`theta10'))+(`d1c1'*`d2c1'*`d3p1'*(1-`theta10'))) if
$ML_y1==1 & $ML_y2==1 & $ML_y3==1
quietly replace
`lnf'=ln((`d1c0'*`d2c0'*`d3d0'*(`theta10'))+(`d1c1'*`d2c1'*`d3d1'*(1-`theta10'))) if
$ML_y1==1 & $ML_y2==1 & $ML_y3==0
end
ml model lf likefn (f1sch= byinc ed_mother ed_father nosib g8farm g8south
brokenf1 fem asian hisp black birthye)
(f2sch= f2inc ed_mother ed_father nosib g10farm g10south brokenf2 fem
asian hisp black birthye)
(f3sch= f2inc ed_mother ed_father nosib g10farm g10south brokenf2 fem
asian hisp black birthye)
(ged: f2inc ed_mother ed_father nosib g10farm g10south brokenf2 fem asian
hisp black birthye)
(purs: f2inc ed_mother ed_father nosib g10farm g10south brokenf2 fem
asian hisp black birthye)
(alpha2: ) (alphaH: ) (alphaG: ) (alphaP: ) (p: )
ml search p 0 1
ml maximize, difficult
OUTPUT from Stata:
initial: log likelihood = -7134.5371
rescale: log likelihood = -6703.5143
rescale eq: log likelihood = -5861.5134
numerical derivatives are approximate
flat or discontinuous region encountered
numerical derivatives are approximate
flat or discontinuous region encountered
Iteration 0: log likelihood = -5861.5134 (not concave)
numerical derivatives are approximate
flat or discontinuous region encountered
numerical derivatives are approximate
flat or discontinuous region encountered
ALL OTHER ITERATIONS HAVE SAME PROBLEM
Iteration 6: log likelihood = -5231.9467 (not concave)
numerical derivatives are approximate
flat or discontinuous region encountered
Iteration 7: log likelihood = -5231.9467 (not concave)
Number of obs =
8045
Wald chi2(12) =
233.62
Log likelihood = -5231.9467 Prob > chi2 =
0.0000
------------------------------------------------------------------------------
| Coef. Std. Err. z P>|z| [95%
Conf. Interval]
-------------+----------------------------------------------------------------
eq1 |
byinc | .1655285 .0292439 5.66 0.000 .1082115
.2228455
OTHER VARIABLES OMMITD, BUT THE OUTPUT IS FINE
_cons | 2.635513 7.461589 0.35 0.724 -11.98893
17.25996
-------------+----------------------------------------------------------------
eq2 |
f2inc | .2129379 .0412382 5.16 0.000 .1321124
.2937633
OTHER VARIABLES OMMITD, BUT THE OUTPUT IS FINE
_cons | -.3608168 9.079173 -0.04 0.968 -18.15567
17.43404
-------------+----------------------------------------------------------------
eq3 |
f2inc | .3751662 .0635452 5.90 0.000 .2506199
.4997126
OTHER VARIABLES OMMITED, BUT THE OUTPUT IS FINE
_cons | 1.189953 13.48772 0.09 0.930 -25.24549
27.6254
-------------+----------------------------------------------------------------
ged |
f2inc | .0868075 .036056 2.41 0.016 .016139
.1574759
OTHER VARIABLES OMMITED, BUT THE OUTPUT IS FINE
_cons | .6083109 8.647009 0.07 0.944 -16.33952
17.55614
-------------+----------------------------------------------------------------
purs |
f2inc | .0397004 .0371077 1.07 0.285 -.0330293
.1124301
OTHER VARIABLES OMMITED, BUT THE OUTPUT IS FINE
_cons | -.0012283 9.053507 -0.00 1.000 -17.74578
17.74332
-------------+----------------------------------------------------------------
alpha2 |
_cons | 41.03717
. . . . .
-------------+----------------------------------------------------------------
alphaH |
_cons | 26.21788
. . . . .
-------------+----------------------------------------------------------------
alphaG |
_cons | -31.41291
. . . . .
-------------+----------------------------------------------------------------
alphaP |
_cons | 15.89986
. . . . .
-------------+----------------------------------------------------------------
p |
_cons | .1606487 .0064533 24.89 0.000 .1480003
.173297
------------------------------------------------------------------------------
*
* 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/
