Re: st: Changes in Stata's ml routine d0? Stata 8.2 vs. Stata 11.2

[Nick Cox <njcoxstata@gmail.com>]

Changes in -ml- (and much else) are documented in the whatsnew* help
files. Start with -help whatsnew- and scroll to the bottom of the
viewer for an index.

Nick

On Mon, Mar 5, 2012 at 8:20 AM,  <Phil1899@gmx.de> wrote:

> I would like to ask whether there have been any major changes in Stata’s ml routine during the last six years (in particular in case of d0)? I am wondering about that as I am trying to replicate results from a Stata Journal article from the year 2006 describing a code for Stata’s ml routine. I am using the same code and the same data as in the article, however, the program no longer converges (not concave). I am working with Stata 11.2, the authors of the Stata Journal article used Stata 8.2.
>
> The article I am referring to is the following: Peter Haan & Arne Uhlendorff, 2006. "Estimation of multinomial logit models with unobserved heterogeneity using maximum simulated likelihood," Stata Journal, vol. 6(2), pages 229-245.
>
> The data used in the article is available at: http://www.gllamm.org/jspmix.dat
>
> I present the code below. I would appreciate any advice in this matter.
>
> Best,
> Phil
>
>
> clear
> clear matrix
> clear mata
> set mem 400m
> set more off
> cd C:\temp
> set matsize 5000
>
> cd ""
>
> matrix p = (7, 11)
> global draws "50"
> infile scy3 id sex stag ravi fry3 tby using jspmix.dat, clear
> keep scy3
> sort scy3
>
> by scy3: keep if _n==1
> mdraws, neq(2) dr($draws) prefix(c) burn(10) prime (p)
>
> local repl=${draws}
> local r=1
> while `r' <= `repl' {
> by scy3: gen random_1`r'=invnorm(c1_`r')
> by scy3: gen random_2`r'=invnorm(c2_`r')
> local r=`r'+1
> }
> sort scy3
> save mdraws_${draws}, replace
>
> clear
> infile scy3 id sex stag ravi fry3 tby using jspmix.dat, clear
> sort scy3
> merge m:1 scy3 using mdraws_${draws}.dta, update
> drop _merge
> sort scy3
>
> * starting values
> mlogit tby sex, base(1)
> matrix init= e(b)
> matrix list init
> local c1 = init[1,3]
> local c2 = init[1,4]
> local c3 = init[1,5]
> local c4 = init[1,6]
> matrix C = (`c1', `c2', `c3', `c4')
> matrix list C
>
> gen a1=0
> gen a2=0
> gen a3=0
> replace a1=1 if tby==1
> replace a2=1 if tby==2
> replace a3=1 if tby==3
> sort scy3
>
> program drop _all
> program define mlogit_sim_d0
> args todo b lnf
> tempvar etha2 etha3 random1 random2 lj pi1 pi2 pi3 sum lnpi L1 L2 last
> tempname lnsig1 lnsig2 atrho12 sigma1 sigma2 cov12
> mleval `etha2' = `b', eq(1)
> mleval `etha3' = `b', eq(2)
> mleval `lnsig1' = `b', eq(3) scalar
> mleval `lnsig2' = `b', eq(4) scalar
> mleval `atrho12' = `b', eq(5) scalar
>
> qui {
> scalar `sigma1'=(exp(`lnsig1'))^2
> scalar `sigma2'=(exp(`lnsig2'))^2
> scalar `cov12'=[exp(2*`atrho12')-1]/[exp(2*`atrho12')+1]*(exp(`lnsig2'))*(exp(`lnsig1'))
> gen double `random1' = 0
> gen double `random2' = 0
> gen double `lnpi'=0
> gen double `sum'=0
> gen double `L1'=0
> gen double `L2'=0
> by scy3: gen byte `last'=(_n==_N)
> gen double `pi1'= 0
> gen double `pi2'= 0
> gen double `pi3'= 0
> }
> matrix W = ( `sigma1' , `cov12' \ `cov12' , `sigma2')
> capture matrix L=cholesky(W)
> if _rc != 0 {
> di "Warning: cannot do Cholesky factorization of rho matrix"
> }
> local l11=L[1,1]
> local l21=L[2,1]
> local l22=L[2,2]
>
> local repl=${draws}
> local r=1
> while `r' <= `repl' {
> qui {
> replace `random1' = random_1`r'*`l11'
> replace `random2' = random_2`r'*`l22' + random_1`r'*`l21'
>
> replace `pi1'= 1/(1 + exp(`etha2'+`random1')+exp(`etha3'+`random2'))
> replace `pi2'= exp(`etha2'+`random1')*`pi1'
> replace `pi3'= exp(`etha3'+`random2')*`pi1'
>
> replace `lnpi'=ln(`pi1'*a1+`pi2'*a2+`pi3'*a3)
>
> by scy3: replace `sum'=sum(`lnpi')
> by scy3: replace `L1' =exp(`sum'[_N]) if _n==_N
> by scy3: replace `L2'=`L2'+`L1' if _n==_N
>
> }
> local r=`r'+1
> }
> qui gen `lj'=cond(!`last',0, ln(`L2'/`repl'))
> qui mlsum `lnf'=`lj'
> if (`todo'==0|`lnf'>=.) exit
> end
>
>
> ml model d0 mlogit_sim_d0 ( tby = sex) ( tby = sex) /lnsig1 /lnsig2 /atsig12
> ml init C 0.5 0.5 0.5, copy
> ml maximize

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