A Practical Guide to Simulating Correlated Binary Outcomes

Abstract: Simulating dependent Bernoulli outcomes with prescribed means and pairwise Pearson correlations is a common task in risk modeling. A familiar approach is the Gaussian-threshold workflow for binary outcomes, often viewed as a Bernoulli analogue of the Gaussian copula construction. We show that setting latent Gaussian correlations equal to target Bernoulli correlations is generally incorrect after thresholding, and that pairwise tetrachoric calibration is exact only when the calibrated latent matrix is positive semidefinite. We therefore formulate the problem directly over the joint Bernoulli probability mass function. Given target means and pairwise correlations, we impose normalization, nonnegativity, mean constraints, and pairwise cross-moment constraints as a linear program over the $2^N$ atomic probabilities. The resulting PMF formulation either returns an exact law matching the requested first and second moments or certifies infeasibility. A convex-hull characterization further shows that every feasible target admits a law supported on at most $1+N+\binom{N}{2}$ states, while every infeasible target admits a separating quadratic certificate. We then develop a truncated-moment completion scheme that fits a reduced cross-moment table and generates samples by sequential conditioning, together with a sparse-support working-set refinement that can reduce memory usage on structured instances, although the worst-case complexity remains exponential. Together, these constructions provide an exact PMF-based framework for feasibility and simulation at moderate dimension and structured alternatives when the full atomic representation is impractical, while clarifying the limits of Gaussian-threshold constructions.

• The paper shows that calibrated Gaussian-threshold simulation is exact if and only if the resulting latent correlation matrix is positive semidefinite, while naive calibration can substantially miss target correlations.
• The PMF linear program directly tests whether means and pairwise correlations define a valid joint Bernoulli distribution, producing an exact sampler or a global infeasibility certificate.
• The paper presents truncated-moment and sparse-support methods for larger systems, but shows that truncation is not universally valid and exact high-dimensional simulation remains computationally difficult.

The paper "A Practical Guide to Simulating Correlated Binary Outcomes" (2607.16801) addresses a problem that arises routinely in credit risk, clinical development portfolios, and option modeling: given a mean vector $\vec{\mu}$ and a target Pearson correlation matrix $\mathbf{R}$ for $N$ Bernoulli variables, construct a simulator whose output matches those first two moments exactly. The paper's central argument is that the widely used latent-Gaussian threshold workflow does not answer this question correctly in general, and that the right object to work with is the joint probability mass function (PMF) itself. The authors formalize this as a linear feasibility program over the $2^N$ atomic probabilities — the PMF-LP — which either returns an exact law or certifies infeasibility.

The failure of Gaussian-threshold constructions

The paper first dissects the Gaussian-threshold approach, which it names CoMBO ("copula method for binary outcomes"). In its naive form, one draws $Z \sim \mathcal{N}(0, \mathbf{R}_{\text{lat}})$ with $\mathbf{R}_{\text{lat}} = \mathbf{R}$ and thresholds at probit cutoffs. Because thresholding is many-to-one rather than an invertible marginal transform, copula margin-freeness no longer applies: the induced Bernoulli Pearson correlation is

$$\operatorname{Corr}(X_i,X_j) = \frac{\Phi_2(t_i,t_j;\, r_{ij}) - p_i p_j}{\sqrt{p_i(1-p_i)\,p_j(1-p_j)}},$$

a nonlinear function of the latent correlation, so naive CoMBO systematically misses the target correlations. A three-variable simulation experiment makes the magnitude concrete: for targets $(\rho_{12},\rho_{13},\rho_{23}) = (0.1, 0.4, 0.8)$ with means $(0.2, 0.7, 0.6)$, uncalibrated CoMBO produced sample correlations of approximately $(0.052, 0.212, 0.567)$ over 1,000 replications of 10,000 draws — errors far larger than those of the analogous continuous logistic-marginal benchmark, which recovers targets almost exactly after calibration.

Pairwise tetrachoric calibration repairs each pair individually via the strictly increasing map $\mathbf{R}$0, yielding a unique candidate latent matrix $\mathbf{R}$1. The paper's key structural claim is that calibrated CoMBO is exact if and only if $\mathbf{R}$2 is positive semidefinite. For the illustrative example, the calibrated latent matrix has determinant approximately $\mathbf{R}$3, so no exact Gaussian-threshold representation exists even though the target is pairwise Bernoulli-feasible under the Chaganty–Joe bounds. Any PSD repair necessarily perturbs at least one induced Bernoulli correlation. This cleanly separates two distinct questions that are often conflated: whether the requested moments admit any joint Bernoulli law, and whether they admit one within the Gaussian-threshold family. A target can be feasible yet outside the latter family.

The PMF linear program

To answer the feasibility question directly, the paper works over the atomic probabilities $\mathbf{R}$4. Since means and pairwise cross-moments are linear in $\mathbf{R}$5 — $\mathbf{R}$6 and $\mathbf{R}$7 where $\mathbf{R}$8 — matching prescribed moments reduces to a pure feasibility LP with $\mathbf{R}$9 variables and only $N$0 equality constraints. A feasible solution yields an exact simulator by inverse transform sampling; infeasibility is a definitive global certificate, something pairwise bounds cannot provide.

A convex-hull characterization sharpens both directions. Feasibility is equivalent to membership of the target moment vector $N$1 in $N$2, where $N$3 collects $N$4. By Carathéodory's theorem, every feasible target admits a compatible law supported on at most $N$5 states; every infeasible target admits a separating quadratic certificate $N$6 that is nonnegative on all of $N$7 but has negative expectation under the requested moments. The support bound explains why sparse solutions exist; it does not, however, make finding them easy.

Truncated completion and sequential conditioning

For settings where the full state space is impractical, the paper develops a truncated-moment completion scheme indexed by subsets of $N$8. Fixing a truncation order $N$9, the pair moments are pinned to the targets, moments of order $2^N$0 through $2^N$1 become decision variables, and moments above order $2^N$2 are set to zero. By Möbius inversion on the Boolean lattice, the implied atoms satisfy $2^N$3, so validity reduces to a linear feasibility problem over $2^N$4 variables. Sampling then proceeds by sequential conditioning: cylinder probabilities $2^N$5 give conditional ratios computable in at most $2^N$6 terms per step, polynomial in $2^N$7 for fixed $2^N$8.

Two propositions bound what truncation can achieve. First, there exist pairwise-consistent specifications (three fair variables, all pairwise cross-moments zero) that no choice of third-order moments can repair — pairwise feasibility is genuinely insufficient for joint feasibility. Second, and more restrictively, any order-$2^N$9 truncation forces support on states with at most $Z \sim \mathcal{N}(0, \mathbf{R}_{\text{lat}})$0 ones, hence requires $Z \sim \mathcal{N}(0, \mathbf{R}_{\text{lat}})$1. Whenever the mean sum is $Z \sim \mathcal{N}(0, \mathbf{R}_{\text{lat}})$2 — which holds with high probability for i.i.d. means with positive mean, by Hoeffding's inequality — feasibility within the truncated family demands $Z \sim \mathcal{N}(0, \mathbf{R}_{\text{lat}})$3, so $Z \sim \mathcal{N}(0, \mathbf{R}_{\text{lat}})$4 cannot suffice in general. The truncated method is therefore a model class, not a universal feasibility test: infeasibility may reflect too small a truncation order rather than true impossibility.

Sparse-support refinement

The gap between $Z \sim \mathcal{N}(0, \mathbf{R}_{\text{lat}})$5 variables and $Z \sim \mathcal{N}(0, \mathbf{R}_{\text{lat}})$6 constraints motivates a working-set scheme: solve a restricted master LP with slacks over a small candidate support, use dual multipliers to define a quadratic pricing score $Z \sim \mathcal{N}(0, \mathbf{R}_{\text{lat}})$7, add improving states via binary quadratic heuristics (local search, simulated annealing, QUBO-style methods), and iterate until slack vanishes. Combined with the Carathéodory bound, this guarantees that some exact law uses at most $Z \sim \mathcal{N}(0, \mathbf{R}_{\text{lat}})$8 states if the procedure converges to slack zero, though identifying that support offers no worst-case escape from exponential complexity. The paper also notes that naive variable-partitioning decompositions are not exact unless conditional moment constraints per prefix are enforced, which reintroduces exponential cost.

Limitations and open questions

The limitations are stated plainly by the authors. The exact PMF-LP scales exponentially in $Z \sim \mathcal{N}(0, \mathbf{R}_{\text{lat}})$9 because the input itself grows as $\mathbf{R}_{\text{lat}} = \mathbf{R}$0; it is practical only at moderate dimension. The truncated method sacrifices universality, and the sparse-support method retains worst-case exponential behavior despite often reducing memory and runtime on structured instances. Open questions left explicitly include symmetry reduction, cut generation, and column-generation schemes that avoid materializing all atoms; controlled relaxations enforcing critical moments exactly while bounding error elsewhere; adaptive rules for selecting the truncation order from data; and a systematic empirical comparison of truncated versus sparse-support representations.

Conclusion

The paper contributes a clean decomposition of a familiar applied problem into three logically distinct issues: Bernoulli-feasibility of the requested moments, representability within the Gaussian-threshold family, and computational cost of constructing an exact law. Its main results — the exactness condition for tetrachorically calibrated CoMBO, the PMF-LP with its convex-hull characterization and separating certificates, and the polynomial-per-step sequential sampler under fixed-order truncation — together clarify when standard approximations are safe and when they silently fail. The framework is exact and transparent at moderate dimension, while high-dimensional exact moment matching remains open.