# st: maximum simulated likelihood

 From Jenkins S P To statalist@hsphsun2.harvard.edu Subject st: maximum simulated likelihood Date Sun, 28 Aug 2005 14:10:31 +0100 (BST)

```Date: Sat, 27 Aug 2005 12:15:36 -0500
From: "Jun Xu" <mystata@hotmail.com>
Subject: st: simulate maximum likelihood estimation--where to
construct the random variables

A while ago, Professor Stephen P. Jenkins and Arne Uhlenforff
on this:

Arne Uhlenforff <auhlendorff@diw.de> wants to estimate a MNL
model with a
random intercept using simulated ML (rather than using -gllamm-).

**********
I haven't looked at your program in detail, but one thing
that strikes me is
that I think you are not treating the pseudo-random uniform
draws correctly.
You want 50 replications/draws per ML iteration, but are
creating them anew
every iteration -- you shouldn't. They should be created just
once. (See
e.g. the gory details of the code underlying -mvprobit- or
associated SJ
article.) I don't know if this is having knock-on effects elsewhere.

good luck
Stephen
***********

I am also trying to write some codes on simulated ml as
practice. I did it
in two ways, a first method is to generate 100 random
variables outside the
loop as below; a second method is to create `ed`i' (the
random variable)
anew in each iteration.

Method 1:
forval i = 1/100   {
tempname ed`i'
gen double `ed`i'' = ...whatever ways to generate the random
variable
}

qui forval i = 1/100 {
replace `sp'= `sp' * (invlogit(`xb' + `ed`i'')) if
\$ML_y1 == 1
replace `sp'= `sp' * (invlogit(-(`xb' + `ed`i'')))
if \$ML_y1 == 0
}

qui replace `lnf' = ln(`sp')

Method 2:
forval i = 1/100   {
tempname ed`i'

}

qui forval i = 1/100 {
gen double `ed`i'' = ...whatever ways to generate the random
variable
replace `sp'= `sp' * (invlogit(`xb' + `ed`i'')) if
\$ML_y1 == 1
replace `sp'= `sp' * (invlogit(-(`xb' + `ed`i'')))
if \$ML_y1 == 0
}

qui replace `lnf' = ln(`sp')

Strangely enough Method 1 will converge and Method 2 will
never converge. I
am not sure if this is what Professor Jenkins talked about
and why it is so?
I will appreciate any thoughts on this.
```
There were 2 points raised by Arne Uhlendorff and myself in the earlier discussion:
1. the need to create the pseudo-random variables "outside" the -ml- iteration loop, doing it once and for all ab initio. [Your method 1.]
2. (Arne's follow-up point) the need to average the simulated
likelihoods (the average is over the simulations) per -ml- iteration.

Remember that the method is really "maximum simulated likelihood" (not "simulated maximum likelihood"). The draws, outside the loop as you phrase it, are used to calculate the simulated likehood at each iteration. That calculation also requires averaging. Your code appears not to average. If you look inside "mvprob_ll.ado", the subroutine associated with the -mvprobit- program, you will see a line like
replace `lnf' = `lnf' + `sp\$S_MLE_M'/\$S_MLE_D This is related to the averaging process.

A third point, implicit in Arne's discussion, was the use of another program, -gllamm- in his case, to check (or more generally, "certify") that his program was giving the results expected. [I conjecture that Arne was writing his own specialist program because -gllamm- was slow-running for his particular problem.] An alternative is to generate some simulated data and equations with 'known' coefficients, and then to estimate this with your model. [There is an example of this in the Stata Journal SJ 3-3 article about -mvprobit-. The program itself was updated in SJ 5-2.]
When you have got your program running, I suggest some of this sort of checking.

Stephen
=============================================
Professor Stephen P. Jenkins <stephenj@essex.ac.uk>
Institute for Social and Economic Research (ISER)
University of Essex, Colchester CO4 3SQ, UK
Phone: +44 1206 873374. Fax: +44 1206 873151.
http://www.iser.essex.ac.uk
Survival Analysis Using Stata: http://www.iser.essex.ac.uk/teaching/degree/stephenj/ec968/index.php

*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/