Notice: On March 31, it was **announced** that Statalist is moving from an email list to a **forum**. The old list will shut down at the end of May, and its replacement, **statalist.org** is already up and running.

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

From |
Phil1899@gmx.de |

To |
statalist@hsphsun2.harvard.edu |

Subject |
st: Random Effects Probit by Maximum Simulated Likelihood |

Date |
Fri, 24 Feb 2012 16:04:13 +0100 |

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/

**Follow-Ups**:**RE: st: Random Effects Probit by Maximum Simulated Likelihood***From:*Cameron McIntosh <cnm100@hotmail.com>

- Prev by Date:
**st: mosaic plot with an intensity dimension** - Next by Date:
**st: RE: mosaic plot with an intensity dimension** - Previous by thread:
**st: mosaic plot with an intensity dimension** - Next by thread:
**RE: st: Random Effects Probit by Maximum Simulated Likelihood** - Index(es):