RE: st: Random Effects Probit by Maximum Simulated Likelihood

[Cameron McIntosh <cnm100@hotmail.com>]

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/