Dear Users,
here is my problem:
I am running a multivariate probit using Cappellari and Jenkins 2006 code.
At the end of the email I have enclosed their code.
I’m using individual level data (h) and my purpose is to obtain a
prediction of the probabilities of every possible outcomes for each h.
To do that I saved the variable called `sp' in the Cappellari and
Jenkins 2006 code, naming it prob.
In this way after the estimation I obtain a new column in the dataset, which - I guess - represents the predicted probability of the
observed combination of outcomes.
Plese see these few results as an example:
h z1 z2 z3 Prob
1 0 0 1 0.424
2 1 0 0 0.534
3 1 0 1 0.327
the last row contains the probability Pr(Z1=1, Z2=0, Z3=1)=0.327. This means that 1- Pr(Z1=1, Z2=0, Z3=1)=0.673 and represents the sum of the probabilities of all the other possible combinations for this individual. However, I cannot disentangle the specific probabilities of the other combinations, i.e. I don’t know what is the proabability of observing Z1=1, Z2=1, Z3=0 for guy #3.
In few words I don’t know how to get the probabilities of all the
possible combination of outcomes for each individual (i.e. seven new
columns in the dataset).
Does anyone have an idea on how to modify the code to get these
probabilities?
thanks
Charlotte
Here is the code for a multivarite probit with 3 equations:
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 - `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')
replace prob=`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(100) neq(3) prefix(z) random seed(123456789) antithetics replace
/*Code for doing this using ml evaluation method lf is set out below*/
gen prob=.
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/
