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

[Phil Clayton <philclayton@internode.on.net>]

Have you tried running it under version control? eg by putting:
version 8.2
at the top of your do-file.

Phil

On 05/03/2012, at 8:43 PM, 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/
> 
> -- 
> NEU: FreePhone 3-fach-Flat mit kostenlosem Smartphone!                                  
> Jetzt informieren: http://mobile.1und1.de/?ac=OM.PW.PW003K20328T7073a
> *
> *   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/


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