Dear Nick,
really thank you now it works perfectly!
May I ask you once again another suggestion. I'm trying to compare the results obtained using ml (paragraph 3.4 cappellari and jenkins ) with those obtained using mvprobit instead of ml.
I just would like to be sure that I’m estimating the predicted probabilities of multiple outcome variables after mvprobit correctly.
hence my question is: Does anyone see something wrong in the following of proceeding?
I followed what Cappellari and Jenkins SJ 2006 suggest at page 166, that is:
• Step 1→ estimate the multivariate probit model using mvprobit
• Step 2→ save the upper integration point variables – i.e. the 3
(M=3) linear index variables Im =bM’xM which can be derived by
using mvppred xb after mvprobit – In this way I obtained xb1, xb2,
xb3.
• Step 3→ save the V matrix. I used mvppred stdp after mvprobit. I
obtained stdp1, stdp2, stdp3. Then I used correlate stdp1 stdp2 stdp3
and I named this matrix V.
• Step 4→ create the Cholesky matrix using: matrix c=cholesky(V)
• Step 5→ create the signs variable k1=2*y1-1, k2=2*y2-1, k3=2*y3-1
• Step 6→ use mdraws to generate the draws variables – e.g.:mdraws,
dr(250) neq(3) prefix(z) random seed(123456789) antithetics replace
• Step 7→ calculate the probabilities by using the egen command with
the Im variables as arguments and referring to a matrix equal to cholesky(V ) in the
chol() option. i.e.:
egen pr_joint=mvnp( xb1 xb2 xb3), dr(100) chol(V) prefix(ff) signs(k1 k2 k3)
--- On Wed, 7/8/09, Nick Cox <[email protected]> wrote:
> From: Nick Cox <[email protected]>
> Subject: st: RE: RE: AGAIN on mvnp code Cappellari & Jenkins 2006
> To: [email protected]
> Date: Wednesday, July 8, 2009, 1:00 PM
> You definitely should not change the
> code
>
> egen `sp' =
>
> to
>
> egen sp =
>
> as that produces the error you report.
>
> From my understanding, what I would suggest is
>
> 0. Before you call this code,
>
> gen safecopy = .
>
> 1. Within your -myll- code add as last line
>
> replace safecopy = `sp'
>
> 2. Access the variable -safecopy- afterwards.
>
> Nick
> [email protected]
>
>
> Charlotte Gary
>
> My purpose is to get the probabilities of every possible
> multiple outcome, therefore sp in the code of Cappellari and
> Jenkins 2006 (p.168) is the variable of my interest.
> However, in the code (see again below) sp is defined
> as a tempvar hence it is not saved as a new variable in the
> dataset. I tried to modify the program omitting sp from the
> tempvar list and changing this line:
> egen `sp'= mvnp(`xb1' `xb2' `xb3'), chol(`C') draws($dr)
> prefix(z) signs(`k1' `k2' `k3' )
>
> into this:
> egen sp= mvnp(`xb1' `xb2' `xb3'), chol(`C') draws($dr)
> prefix(z) signs(`k1' `k2' `k3' )
>
> However, this was of proceeding is not correct as after ml
> maximize an error message appears: “sp is already
> defined”.
> Hence my question is: how is it possible to obtain the
> probabilities of every possible multiple outcome at unit of
> analysis level?
> Really thank you in advance!
>
>
> program define myll
> args lnf xb1 xb2 xb3 c21 c31 c32
> tempvar sp k1 k2 k3
> quietly {
> gen double `k1' = 2*$ML_dummy_C - 1
> gen double `k2' = 2*$ML_dummy_I - 1
> gen double `k3' = 2*$ML_dummy_S - 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 - `cf21'^2 )
> scalar `cf33' = sqrt( 1 - `cf31'^2 - `cf32'^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)/c21 /c31 /c32, title("MV Probit by MSL, $dr pseudorandom
> draws")
> ml init b0
> ml maximize
>
> *
> * 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/
>
*
* 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/
