st: Bivariate Random Effects Probit by Simulated ML

[Phil1899@gmx.de]

Dear all,

I want to estimate a bivariate random effects probit model by maximum simulated likelihood. I follow Haan and Uhlendorf (2006) and Cappellari and Jenkins (2006) as far as possible in writing the program. Unfortunately, the estimator does not converge. In particular, there seems to be a problem with the Cholesky factorization. I am really stuck here as I generate a data set to perform the estimations. Perhaps anyone on the list has an idea about what I am doing wrong? I would greatly appreciate help in this matter! Please find below the basic model setup and the program code.   

Best
Phil

Cappellari & Jenkins (2006), Calculation of Multivariate Normal Probabilities by Simulation, with Applications to Maximum Simulated Likelihood Estimation. IZA Discussion Paper Series.

Haan & Uhlendorff (2006), Estimation of Multinomial Logit Models with Unobserved Heterogeneity Using Maximum Simulated Likelihood. DIW Discussion Papers.

The basic model looks as follows:

Y1_it = b1x1_it + u1_i + v1_it
Y2_it = b2x1_it + b3x2_it + u2_i + v2_it

With  
Cov u_i = (sig_u1 , rho1* sig_u1* sig_u2 \ rho1* sig_u1* sig_u2 , sig_u2 ),  where sig_u1 stands for the variance of u1 and rho1 is the correlation of unobserved heterogeneity and Cov v_it = (1 , rho2 \ rho2, 1 )  

* create data set
set seed 123456789
set obs 500
gen pid = _n
global draws "50"

* generate halton sequences:
matrix p = (7, 11)
mdraws, neq(2) dr($draws) prefix(c) burn(10) prime(p)

local repl=${draws}
local r=1
sort pid
while `r' <= `repl' {
by pid: gen random_1`r'=invnorm(c1_`r')
by pid: gen random_2`r'=invnorm(c2_`r')
local r=`r'+1
}
sort pid
drop c*
save mdraws_${draws}, replace

drop random*

* generate individual specific effects
matrix R = (1, .5 \ .5, 1)
drawnorm u1 u2 , corr(R)
expand 5
bysort pid: gen t = _n
xtset pid t
xtsum
matrix R = (1, .3 \ .3, 1)
drawnorm v1 v2 , corr(R)

ge x1 = uniform()-.5
ge x2 = uniform() + 1/3

* Equations
ge y1s = .5 + 4*x1 + u1 + v1
ge y2s = 3 + .5*x1 - 3*x2 + u2 + v2
ge y1 = y1s>0
ge y2 = y2s>0
tab y1 y2
sort pid
merge m:1  pid using mdraws_${draws}.dta
drop _merge
sort pid


program drop _all
program myll_d0
args todo b lnf
tempvar xb1 xb2 random1 random2 lj atrho1 atrho2 lnsig1 lnsig2 
tempname  sig1 sig2 rho1 rho2 rho3 xb1_re xb2_re sp sum L1 last k1 k2
mleval `xb1' = `b' , eq(1) 				
mleval `xb2' = `b' , eq(2) 				
mleval `lnsig1' = `b' , eq(3) scalar	/*sigma_u1*/
mleval `lnsig2' = `b' , eq(4) scalar	/*sigma_u2*/
mleval `atrho1' = `b' , eq(5) scalar	/* rho_u*/
mleval `atrho2' = `b' , eq(6) scalar	/* rho_v*/

qui {
scalar `sig1' = (exp(`lnsig1'))^2 
scalar `sig2' = (exp(`lnsig2'))^2
scalar `rho1'=tanh(`atrho1')*(exp(`lnsig2'))*(exp(`lnsig1'))
scalar `rho2'=tanh(`atrho2')

gen double `rho3'=0
gen double `random1' = 0
gen double `random2' = 0
gen double `sp' = 0
gen double `xb1_re' = 0
gen double `xb2_re' = 0
gen double `sum'=0
gen double `L1'=0
by pid: gen `last' = (_n==_N)

gen double `k1' = 2*$ML_y1 - 1
gen double `k2' = 2*$ML_y2 - 1
}


matrix W = ( `sig1' , `rho1' \ `rho1' , `sig2' )
capture matrix L = cholesky(W)
if _rc != 0 {
di "Warning: cannot do Cholesky factorization of rho matrix"
}
local l11=L[1,1]
local l21=L[2,1]
local l22=L[2,2]

local r=1
while `r' <= $draws  {
qui {
replace `random1' = random_1`r'*`l11'
replace `random2' = random_2`r'*`l22' + random_1`r'*`l21'

replace `xb1_re' = `k1'*(`xb1' + `random1')
replace `xb2_re' = `k2'*(`xb2' + `random2')
replace `rho3' = `k1' * `k2' * `rho2'

replace `sp' = binormal(`xb1_re',`xb2_re',`rho3')

by pid: replace `sum'=sum(`sp'[_N]) if _n==_N
by pid: replace `L1'=`L1'+`sum' if _n==_N
}
local r = `r' + 1
}

qui: gen double `lj'=cond(`last',ln(`L1' / 50),0)
qui: mlsum `lnf'=`lj'
if (`todo'==0|`lnf'>=.) exit
end

qui {
xtprobit y1 x1,re
mat b1 = e(b)
mat B1 = (b1[1,1],b1[1,2])
mat coleq B1 = y1
xtprobit y2 x1 x2,
mat b2 = e(b)
mat B2 = (b2[1,1],b2[1,2],b2[1,3])
mat coleq B2 = y2
mat b0 = (B1, B2)
mat list b0
}
ml model d0 myll_d0 (y1: y1=x1) (y2: y2=x1 x2) /lnsig1 /lnsig2 /atrho1 /atrho2  
ml check
ml init b0, skip
ml maximize

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