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

From |
Charlotte Gary <charlottegary83@yahoo.com> |

To |
statalist@hsphsun2.harvard.edu |

Subject |
st: mvnp code Cappellari & Jenkins 2006 |

Date |
Tue, 7 Jul 2009 02:47:55 -0700 (PDT) |

Hi all, Does anyone have tried the code suggested by Cappellari and Jenkins Stata Journal (2006)in paragraph 3.4? (I have reported the code at the end) I am simply trying to run this code (without changes) but an error message after “ml maximize” impedes to obtain the results that the authors show on p. 170. The error message says: ^2 invalid name I check the code several times but I cannot see the error. I suspect that there is something that stata doesn’t like here: scalar `cf22' = sqrt( 1 - `c21'^2 ) Thank you! Charlotte here is the code: set seed 123456789 set obs 5000 matrix R = (1, .25, .5 \ .25, 1, .75 \ .5, .75, 1) drawnorm u1 u2 u3, corr(R) correlate u* gen x1 = uniform()-.5 gen x2 = uniform() + 1/3 gen x3 = 2*uniform() + .5 gen y1s = .5 + 4*x1 + u1 gen y2s = 3 + .5*x1 - 3*x2 + u2 gen y3s = 1 - 2*x1 + .4*x2 -.75*x3 + u3 gen y1 = y1s>0 gen y2 = y2s>0 gen y3 = y3s>0 program define myll args lnf xb1 xb2 xb3 c21 c31 c32 tempvar sp k1 k2 k3 quietly { gen double `k1' = 2*$ML_y1 - 1 gen double `k2' = 2*$ML_y2 - 1 gen double `k3' = 2*$ML_y3 - 1 tempname cf21 cf22 cf31 cf32 cf33 C su `c21', meanonly scalar `cf21' = r(mean) su `c31', meanonly scalar `cf31' = r(mean) su `c32', meanonly scalar `cf32' = r(mean) scalar `cf22' = sqrt( 1 - `c21'^2 ) scalar `cf33' = sqrt( 1 - `c31'^2 - `c32'^2 ) mat `C' = (1, 0, 0 \ `cf21', `cf22', 0 \ `cf31',`cf32', `cf33') egen `sp'=mvnp(`xb1'`xb2'`xb3'), chol(`C') draws($dr) prefix(z) signs(`k1'`k2'`k3') replace `lnf'= ln(`sp') } end quietly { probit y1 x1 mat b1 = e(b) mat coleq b1 = y1 probit y2 x1 x2 mat b2 = e(b) mat coleq b2 = y2 probit y3 x1 x2 x3 mat b3 = e(b) mat coleq b3 = y3 mat b0 = b1, b2, b3 } mdraws, dr(250) neq(3) prefix(z) random seed(123456789) antithetics replace global dr = r(n_draws) ml model lf myll (y1: y1=x1) (y2: y2=x1 x2) (y3: y3 = x1 x2 x3) ml init b0 ml maximize nlcom * * 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/

- Prev by Date:
**AW: st: xtmepoisson error variance** - Next by Date:
**Re: st: mvnp code Cappellari & Jenkins 2006** - Previous by thread:
**st: Frauke Ruether ist außer Haus (out of office)** - Next by thread:
**Re: st: mvnp code Cappellari & Jenkins 2006** - Index(es):

© Copyright 1996–2017 StataCorp LLC | Terms of use | Privacy | Contact us | What's new | Site index |