Thank you Maarten for your suggestion.
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
--- On Tue, 7/7/09, Nick Cox <[email protected]> wrote:
> From: Nick Cox <[email protected]>
> Subject: RE: st: mvnp code Cappellari & Jenkins 2006
> To: [email protected]
> Date: Tuesday, July 7, 2009, 1:57 PM
> The Stata Journal is perfectly happy
> to include substantial do-files with the electronic media
> supporting each article. There is no policy to exclude or
> discourage such files. At the same time, any principle that
> every segment of code in each article should be matched by a
> separate do-file would lead to much more work to no obvious
> good end, and so we don't support it. The Stata Journal also
> associates itself with common-sense advice to users to think
> about their code and check it carefully.
>
> Nick
> [email protected]
>
> Editor, SJ
>
> Martin Weiss
>
> I am a long-time advocate of do-files to accompany the SJ,
> so far w/o any success. That would help spot such problems
> and also make them easy to rectify by just changing the
> do-file...
>
> Nick Cox
>
> A good way to spot that something was wrong was to notice
> that `c21' etc are variables, and thus that the code as
> originally given implied stuffing variables into scalars.
> That in itself is not illegal in Stata, as Stata would use
> `c21'[1] etc., but as a general rule that's unlikely to be
> what you want.
>
> Maarten buis
>
> --- On Tue, 7/7/09, Charlotte Gary wrote:
> > 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.
>
> You'll need the change the lines:
> scalar `cf22' = sqrt( 1 - `c21'^2 )
> scalar `cf33' = sqrt( 1 - `c31'^2 - `c32'^2 )
>
> into:
> scalar `cf22' = sqrt( 1 - `cf21'^2 )
> scalar `cf33' = sqrt( 1 - `cf31'^2 - `cf32'^2 )
>
> and change the line:
> egen `sp'=mvnp(`xb1'`xb2'`xb3'), chol(`C') ///
> draws($dr) prefix(z) signs(`k1'`k2'`k3')
>
> into:
> egen `sp'=mvnp(`xb1' `xb2' `xb3'), chol(`C') ///
> draws($dr) prefix(z) signs(`k1' `k2'
> `k3')
>
> (notice the spaces between variable names)
>
> *
> * 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/
