Notice: On April 23, 2014, Statalist moved from an email list to a forum, based at statalist.org.
From | Nick Cox <njcoxstata@gmail.com> |
To | statalist@hsphsun2.harvard.edu |
Subject | Re: st: Changes in Stata's ml routine d0? Stata 8.2 vs. Stata 11.2 |
Date | Mon, 5 Mar 2012 09:08:43 +0000 |
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/