Dear all,

I have been using CMP, a wonderful package for models whose disturbance terms are multivariate normally distributed, to simultaneously estimate a model with a probit, a multinomial probit and (in some cases) linear models. I'm currently using Stata 16 under Windows 10 and a household survey with around 5000 observations for context.

CMP can now generate predicted probabilities for the multinomial probit part of the model, but they don't appear to sum to 1. Or more precisely, when asking for the probabilities of the modeled outcomes and calculating the probability of the reference outcome as 1 minus the sum of the modeled outcomes, I regularly get a negative predicted probabilty for the reference outcome (and for some observations, the sum of the predicted outcomes of my 4-outcome mprobit is 3). Here is an example of the code used:
Code:
constraint 1 [_outcome_1_2]: ExtraRooms Livestock
constraint 2 [_outcome_1_3]: Certified Livestock
constraint 3 [_outcome_1_4]: Certified ExtraRooms

gen inmprob=Employed*$cmp_mprobit
count
local ndraws=round(sqrt(r(N)*2))
cmp (JobType: part_status = `JobTypeVars') /*
*/ (WageEmp: WageEmp = `RedFormVars')  /*
*/ (Select: Active = `SelectVars') [pw=${weightvar}], /*
*/ ind(inmprob $cmp_probit $cmp_probit) struc ghkdraws(`ndraws', type(halton)) nolr nonrtol constraints(1/3)

levelsof part_status, local(ylevs)
local ny: word count `ylevs'
gen double Pr_1=1
forvalues j=2/`ny' {
     predict Pr_`j', eq(#`j') pr
     replace Pr_1=Pr_1-Pr_`j'
}
gen double sumps = Pr_2 + Pr_3 + Pr_4
sum sumps, d
which gives the following output:
Code:
                            sumps
-------------------------------------------------------------
      Percentiles      Smallest
 1%     .3317354       .2309818
 5%     .4679055       .2455179
10%     .5658572       .2640299       Obs               4,742
25%      .784298       .2640347       Sum of Wgt.       4,742

50%     .9918313                      Mean           1.223272
                        Largest       Std. Dev.      .8293418
75%     1.004787              3
90%            3              3       Variance       .6878078
95%            3              3       Skewness       1.545598
99%            3              3       Kurtosis       3.743591
Perhaps I am using the predict command wrong, but the sum of 3 out of the 4 predicted probabilities shouldn't be greater than 1 for any observations, whereas it's the case for at least 75% of them here (see Q3). Any help would be greatly appreciated.

Thanks.

David