Stata The Stata listserver
[Date Prev][Date Next][Thread Prev][Thread Next][Date index][Thread index]

RE: st: RE: simulated maximum likelihood estimation


From   "Jun Xu" <mystata@hotmail.com>
To   statalist@hsphsun2.harvard.edu
Subject   RE: st: RE: simulated maximum likelihood estimation
Date   Mon, 09 May 2005 17:06:01 -0500

Nick,

Thanks a lot as always! What I am trying to do is to estimate a simulated maximum likelihood logit models. In logit models,

eta = x'b , where eta = ln(p/1-p), and based on soc theories, I suspect
eta = x'b + error, where error is a positively valued disturbance term (say expotentially distributed). Because there is probably no closed form solution, I need to turn to simulated maximum likelihood estimation. That is, in the codes I provided, I try to add an error (expotentially distributed with mean of 5) to x'b, repeat this process for 100 times, then averge them out and we get simulated likelihood (use expectation to approximate the integral, or to integrate out that disturbance term). Then I use ml to maximize the function.

I tried to dummy up mvprobit_ll.ado file, in which they kind of created # of replication times of variables like:

forval i = 1/$S_MLE_M {
gen double `d`i'' = 0
gen double `sp`i'' = 0
gen double `arg`i'' = 0
gen double `k`i'' = 2*``i''-1
}
...
...
...

replace `sp`i'' = normprob((`arg`i'')/`c`i'`i'')*`sp`j''
...


Jun Xu
Ph.D. Candidate
Department of Sociology
Indiana University at Bloomington
http://mypage.iu.edu/~junxu/home




From: "Nick Cox" <n.j.cox@durham.ac.uk>
Reply-To: statalist@hsphsun2.harvard.edu
To: <statalist@hsphsun2.harvard.edu>
Subject: RE: st: RE: simulated maximum likelihood estimation
Date: Mon, 9 May 2005 21:18:25 +0100

As I understand it, you don't need 200
temporary variables to do what you
want, as 2 will suffice. Thus your
code, at a hasty hack, seems to boil down
to this. I would still recommend building
up the log likelihood as a sum of logs, not
as a product to be logged later. I don't
understand what you are trying to do,
so can't help further.

set trace off
set more off
cap program drop my_ll
program my_ll
     version 8.2
     args lnf xb
     tempname b sp0
     sca `b' = 5

     tempvar ed sp
     gen double `sp' = 1
     gen double `ed' = 1

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

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

Nick
n.j.cox@durham.ac.uk

Jun Xu

> Sorry for that mistake. The following is a revised version.
> Still encounter
> difficulty calculating numerical derivatives. I don't think it's data
> problem. I tried to look into mvprob_ll.ado by Dr. Cappellari and Dr.
> Jenkins, and I am stuck there. No one is responsible for
> solving my problem,
> except myself; however, I do need some even slight hint to
> get me through.
>
> ********************************
> set trace off
> set more off
> cap program drop my_ll
> program my_ll
>     version 8.2
>     args lnf xb
>     tempname b sp0
>     sca `b' = 5
>
>     forval i = 1/100 {
>         tempname ed`i' like`j' sp`i'
>         gen double `sp`i'' = 0
>     }
>
>     gen double `sp0' = 1
>
>     qui replace `lnf' = 0
>     forval i = 1/100 {
>         loc j = `i' - 1
>         qui gen double `ed`i'' = -ln(uniform())*`b'
>         qui replace `sp`i''= `sp`j''*(invlogit(
> `xb'+`ed`i''))      if
> $ML_y1 == 1
>         qui replace `sp`i''=
> `sp`j''*(invlogit(-(`xb'+`ed`i'')))    if
> $ML_y1 == 0
>     }
>
>     qui replace `lnf' = ln(`sp100')
> end
>
> use binlfp2.dta, clear
> ml model lf my_ll (lfp = k5 k618 age), technique(nr bhhh dfp bfgs)
> ml search
> ml maximize
>

*
*   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/
_________________________________________________________________
Donít just search. Find. Check out the new MSN Search! http://search.msn.click-url.com/go/onm00200636ave/direct/01/

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




© Copyright 1996–2014 StataCorp LP   |   Terms of use   |   Privacy   |   Contact us   |   What's new   |   Site index