st: Random Effects Probit by Maximum Simulated Likelihood

[Phil1899@gmx.de]

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/