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

[Nick Cox <njcoxstata@gmail.com>]

-ml- was rewritten in Stata 11. A side-effect is that occasionally
some difficult-to-fit models are now more difficult to fit. I haven't
looked at your specifics.

By the way, your code makes assumptions about operating system and
directory structure which make it awkward for many other people to
test it without some rewriting.

Nick

On Mon, Mar 5, 2012 at 9:43 AM,  <Phil1899@gmx.de> wrote:
> Thanks for the hint Nick!
>
> I looked through the file, however, as far as I can understand none of the fixes / changes to ml can explain why a code working for Stata 8.2 does not work for Stata 11.2.  Does anyone have an idea why the code below no longer works?
>
> Best,
> Phil
>
> -------- Original-Nachricht --------
>> Datum: Mon, 5 Mar 2012 09:08:43 +0000
>> Von: Nick Cox <njcoxstata@gmail.com>
>> An: statalist@hsphsun2.harvard.edu
>> Betreff: Re: st: Changes in Stata\'s ml routine d0? Stata 8.2 vs. Stata 11.2
>
>> 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/