Re: RE: st: Random Effects Probit by Maximum Simulated Likelihood

[Phil1899@gmx.de]

Hi,

Thanks for the reply! I actually saw the post on statalist before. The conversation between Stephen Jenkins and Arne Uhlendorf led to an article in Stata Journal in 2006 by Arne Uhlendorf and Peter Haan with the title “Estimation of Multinomial Logit Models with Unobserved Heterogeneity Using Maximum Simulated Likelihood”. In this article they show how to estimate such a model using Stata’s ml routine. I used this paper as a guide to write the little program for the random effects probit model I posted in the first mail. They actually use the log-sum-exponentiate way to get the product which is not working in my case. That is why I asked the question on statalist. Would be great if anyone could provide me with an answer to that.  

Best,
Phil



-------- Original-Nachricht --------
> Datum: Sat, 25 Feb 2012 07:05:05 -0500
> Von: Cameron McIntosh <cnm100@hotmail.com>
> An: STATA LIST <statalist@hsphsun2.harvard.edu>
> Betreff: RE: st: Random Effects Probit by Maximum Simulated Likelihood

> You may find some clues here and in other similar archived threads:
> http://www.stata.com/statalist/archive/2005-08/msg00898.html ;
> 
> ----------------------------------------
> > Date: Fri, 24 Feb 2012 16:04:13 +0100
> > From: Phil1899@gmx.de
> > Subject: st: Random Effects Probit by Maximum Simulated Likelihood
> > To: statalist@hsphsun2.harvard.edu
> >
> > Dear all,
> >
> > I want to estimate a random effects probit model by maximum simulated
> likelihood. I am using Stata’s ml routine for this purpose and generate
> Halton sequences for the simulation using Cappellari and Jenkin’s (2006)
> mdraws program. My question relates to the last part of the estimation
> algorithm where the product over time T has to be taken. The program is very much
> affected by the way I obtain this product:
> >
> > 1) If I use the user written egen-function “prod”, I obtain results
> which are fairly close to the actual parameter estimates. (findit gprod)
> >
> > 2) If I use the alternative of taking logs, summing them up and then
> exponentiating, the model does not work.
> >
> > I cannot really understant why this is and would appreciate help in this
> matter.
> >
> > I estimate the model as suggested by Greene (Econometric Analysis, 5th
> edition, pp.693-694):
> >
> > lnL_simulated= SUM(1toN) ln (1/R SUM(1toR) (PRODUCT(1toT) F (k_it
> (x_it’beta + sigma_a * random_ir)))
> >
> > where R is the number of draws, F is the normal CDF, k_it=(2*y_it-1),
> x_it are the independent variables, sigma_a is the standard deviation of the
> individual effect a and random_id are draws from the standard normal
> distribution.
> >
> > I present the estimation code below. It starts with generating a data
> set and the Halton sequences. Then the program follows. Note that in the code
> below only 10 draws are taken; the same problems arises with more draws.
> >
> > Another question I have is about the estimated parameters obtained when
> using the gprod egen function; they are not that close to the real
> parameters. Is this because of a mistake in the program? Below I also compare
> actual and estimated parameters for 100 draws.
> >
> > Thank you very much!
> > Phil
> >
> >
> > _________Acutal Parameter__Esimated Parameter
> > x1______________0.5_______________0.56
> > x2______________-4.5______________-4.86
> > cons____________3_________________3.24
> > sigma_a_________1_________________1.16
> >
> > clear
> > clear matrix
> > clear mata
> > set mem 400m
> > set more off
> > cd C:\temp
> > set matsize 5000
> >
> > capture log close
> > cd "C:\"
> >
> > * create data set
> > set seed 123456789
> > set obs 500
> > gen pid = _n
> > global draws "10"
> >
> > * generate halton sequences:
> > matrix p = (3)
> > mdraws, neq(1) prime(p) dr($draws) prefix(c) burn(10)
> > local repl=${draws}
> > local r=1
> > sort pid
> > * convert sequnces into standard normal sequences
> > while `r' <= `repl' {
> > by pid: gen random_`r'=invnorm(c1_`r')
> > local r=`r'+1
> > }
> > sort pid
> > drop c*
> > save mdraws_re_${draws}, replace
> > drop random*
> > * generate individual specific effects
> > gen a = rnormal()
> > expand 5
> > bysort pid: gen t = _n
> > * explanatory variables
> > gen u = rnormal()
> > gen x1 = uniform()-.5
> > gen x2 = uniform() + 1/3
> > * data generating process
> > ge ys = 3 + .5*x1 - 4.5*x2 + a + u
> > sum ys
> > ge y = ys>0
> > sum y
> > tab y
> >
> > sort pid
> > merge m:1 pid using mdraws_re_${draws}.dta /* get Halton draws --> they
> vary on individual level, not time!*/
> > drop _merge
> > sort pid
> > * initial values
> > qui{
> > probit y x1 x2
> > mat b1 = e(b)
> > mat coleq b1 = y
> > }
> > program drop _all
> > program define reprob
> > set seed 1234567
> > args todo b lnf
> > tempvar lj k xb mu random_mu T last sum L1 L2 sp random xb_re xb_re_std
> xb_re_m atrho lam sig lnsig lnsp
> >
> > sort pid
> > mleval `xb' = `b' , eq(1)
> > mleval `lnsig' = `b' , eq(2) scalar /* rho_re*/
> > qui{
> > scalar `sig'=exp(`lnsig')
> > local y $ML_y1
> > gen double `sum'=0
> > gen double `L1'=0
> > gen double `L2'=0
> > gen double `sp'=0
> > gen double `random'=0
> > gen double `xb_re'=0
> > by pid: gen `last' = (_n==_N)
> > gen double `k' = (2*$ML_y - 1)
> > gen double `lnsp'=0
> > }
> > local r=1
> > forv r=1/${draws}{
> > qui{
> > replace `random' = `sig'*random_`r'
> > replace `xb_re' = `k'*(`xb' + `random')
> > replace `sp'=normal(`xb_re')
> > /*
> > * use ln, sum and exp to obtain product
> > replace `lnsp'=ln(`sp')
> > by pid: replace `sum'=sum(`lnsp')
> > by pid: replace `L2'=exp(`sum'[_N]) if _n==_N
> > by pid: replace `L1'=`L1'+`L2' if _n==_N
> > */
> > * use product command
> > capture drop `sum'
> > by pid: egen double `sum' = prod(`sp')
> > by pid: replace `L1'=`L1'+`sum' if _n==_N
> > }
> > }
> > qui: gen double `lj'=cond(!`last',0,ln(`L1' /$draws))
> > qui: mlsum `lnf'=`lj'
> > if (`todo'==0|`lnf'>=.) exit
> > end
> > ml model d0 reprob (y: y=x1 x2) /lnsig
> > ml check
> > ml init b1
> > ml maximize
> > nlcom (sig: exp([lnsig]_b[_cons]))
> >
> > --
> > 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/
>  		 	   		  
> *
> *   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/

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