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

[Phil1899@gmx.de]

Hi all,

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


-- 
Empfehlen Sie GMX DSL Ihren Freunden und Bekannten und wir
belohnen Sie mit bis zu 50,- Euro! https://freundschaftswerbung.gmx.de
*
*   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/