# st: Accept-reject procedure and -drawnorm- for simulating likelihood in -ml-

 From Terence Cheng (by way of Keith Dear ) To statalist@hsphsun2.harvard.edu Subject st: Accept-reject procedure and -drawnorm- for simulating likelihood in -ml- Date Sat, 10 Nov 2007 20:32:30 +1100

Dear Statalisters

I am writing a program that calculates the log likelihood value for use in -ml-. The density function that I am maximising does not have a closed form solution as it depends on three random coefficients. Simulating the density function requires drawing from the truncated trivariate normal distribution which I am attempting to do using the accept-reject method. This first involves generating a set of three random numbers using -drawnorm- and dropping the observations that are not within the accepted range, using the accepted observations to simulate the density.

The problem with the accept-reject approach as described above is that I would need to drop the data (using -drop _all-), which also drops the temporary variables that Stata generates for -ml- such as the lnf and coefficient vectors. I attempted to use scalars to capture the values of the coefficient vectors as I require them for the intermediate calculations. I have also tried to preserve the data prior to dropping them and for which Stata automatically restores after the program ends.

There must have been something wrong with my approach because while implementing -ml check-, Stata returns the following message:

"The initial values are not feasible. This may be because the initial values have been chosen poorly or because there is an error in progname and it always returns missing no matter what the parameter values."

I would appreciate if someone could point out where I had gone wrong. Is there a better way to generate draws from a truncated multivariate normal distribution? Thanks.

The beginning segment of my code is as follows:

****************
cap prog drop progname
program define poinorm_mlv3
version 9.0
args lnf theta1 theta2 theta3 theta4 theta5 theta6
/*Capture coefficient vectors as scalars prior to drop _all*/
tempname rho12 rho13 rho23 /*Correlation matrix*/
scalar `rho12'=[exp(2*`theta4')-1]/[exp(2*`theta4')+1]
scalar `rho13'=[exp(2*`theta5')-1]/[exp(2*`theta5')+1]
scalar `rho23'=[exp(2*`theta6')-1]/[exp(2*`theta6')+1]
if (`rho12'==. | `rho13'==. | `rho23'==.) {
qui replace `lnf'=.
exit
}

matrix Sig = (1 , `rho12' , `rho13' \ /* /*Construct corr mat*/
*/`rho12' , 1 , `rho23' \ /*
*/`rho13' , `rho23' , 1 )
matrix symeigen X V = Sig /*Test for positive df*/
if V[1,3] <= 0 {
qui replace `lnf' =.
exit
}

mat M = inv(Sig)
local repl = 5 /*Number of replications*/
local n = _N /*Sample size for density function*/
tempname xb1 xb2 xb3 y1 y2 y3
scalar `xb1' = `theta1' /*Capture xbs*/
scalar `xb2' = `theta2'
scalar `xb3' = `theta3'
scalar `y1' = \$ML_y1 /*Capture depvars*/
scalar `y2' = \$ML_y2
scalar `y3' = \$ML_y3
preserve
qui drop _all /*drop data and variables*/
/* Using -drawnorm- to generate random numbers, accept and reject. */

*
* For searches and help try:
* http://www.stata.com/support/faqs/res/findit.html
* http://www.stata.com/support/statalist/faq
* http://www.ats.ucla.edu/stat/stata/