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

[Phil1899@gmx.de]

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/