I am interested in how the characteristics of policy i and the characteristics of insurance j influence the probability that insurance j is chosen as a participant on policy i. The standard multinomial choice model as McFadden’s choice model isn't the correct one because there aren't only x potential outcomes and one is chosen. Instead, here there are x potential outcomes and any number of them can be chosen.
I think a maximum likelihood probit estimation in a setting in which multiple outcomes can be simultaneously chosen would be the correct one.
One critical component of the analysis is the correlation structure of the error terms within a choice set. In particular, I would allow the errors to be freely correlated across all potential participants on any policy by a given policy holder. Some policy holder has more than one policy, and I allow errors to be correlated for all potential participants on any of the policies.

Code:
ssc install dataex
clear
input byte policy str6 Holder str9 insurance byte choice
1 "Jeff" "Health" 1
2 "Chris" "Mediocre" 1
3 "Jamie" "Health" 1
3 "Jamie" "AB" 1
3 "Jamie" "Medical" 1
4 "Mela" "Obamacare" 1
end
tempfile a
save `a'
expand 2
replace policy = -1 if _n >_N/2
fillin policy insurance
replace choice = 0 if choice ==.
joinby policy using `a', update
duplicates drop policy Holder insurance choice, force

Output:
policy Holder insurance choice
1 Jeff AB 0
1 Jeff Health 1
1 Jeff Medical 0
1 Jeff Mediocre 0
1 Jeff Obamacare 0
2 Chris AB 0
2 Chris Health 0
2 Chris Medical 0
2 Chris Mediocre 1
2 Chris Obamacare 0
3 Jamie AB 1
3 Jamie Health 1
3 Jamie Medical 1
3 Jamie Mediocre 0
3 Jamie Obamacare 0
4 Mela AB 0
4 Mela Health 0
4 Mela Medical 0
4 Mela Mediocre 0
4 Mela Obamacare 1




Probit-Estimation (first try):
meprobit choice, vce(cluster Holder) nogroup
margins, dydx(*)

--> But how can I compute my probit estimation with multiple simultaneous chosen insurances in Stata? Any thoughts are highly welcome.