18 The Multinomial Probit Model

Part IV leaves the logit family. Every model so far — MNL, nested, latent class, mixed, hierarchical — has Gumbel errors somewhere in its bones, chosen in Chapter 7 for the closed form they buy. The multinomial probit (MNP) makes the other choice: multivariate normal errors with a freely estimated covariance matrix. Correlation between any pair of alternatives, unobserved-taste patterns of any shape — nothing pre-specified, everything estimated. The price is the one we have spent two parts learning to pay: the choice probabilities are integrals with no closed form, ever.

That price structure is why the MNP sits here, at the book’s culmination, rather than in Part II where its cousins live. Estimating it requires either the simulation-assisted likelihood machinery of Chapter 12 or the augmentation-and-Gibbs machinery of Chapter 16 — and it is the rare model where both routes are standard, which makes it the perfect capstone demonstration that model and estimation method are separate choices (Chapter 19 runs the head-to-head). This chapter does the groundwork: the model, its subtle identification, and its simulation.

18.1 Why Probit?

The nested logit of Chapter 10 repaired IIA by letting the researcher declare which alternatives share unobserved appeal. The mixed logit of Chapter 13 generated flexible substitution through random coefficients on observed attributes. The MNP attacks the unobservables directly:

\[ U_{ntj} = \mathbf{x}_{ntj}'\boldsymbol{\beta}+ \varepsilon_{ntj}, \qquad \boldsymbol{\varepsilon}_{nt} = (\varepsilon_{nt1}, \ldots, \varepsilon_{ntJ})' \sim N(\mathbf{0}, \Omega) \tag{18.1}\]

with $\Omega$ an unrestricted $J \times J$ covariance matrix. Return to the story that motivated nesting: Windows machines share unobserved OS-related utility. The nested logit encoded that as a nest with an estimated $\lambda$. The MNP simply lets the data put a positive covariance between the Acer and Dell errors — or not — along with any other pattern among the alternatives: no nests to specify, no researcher declaration of who competes with whom, no IIA anywhere.¹ For this part we return to the brand-labeled design of Chapter 10 (one Acer, one Dell, one Apple per task), because a covariance between “the Acer error” and “the Dell error” needs alternatives with stable identities to attach to.

18.2 The Choice Probability Is an Integral

The behavioral rule is unchanged — pick the maximum utility — so the probability that alternative $j$ is chosen is the probability that $J-1$ utility differences all fall on $j$’s side:

\[ P_{ntj} = \Pr\left( U_{ntj} > U_{ntk} \;\; \forall \, k \neq j \right) = \int_{\text{region where } j \text{ wins}} \phi_J(\boldsymbol{\varepsilon}\mid \mathbf{0}, \Omega) \, d\boldsymbol{\varepsilon} \tag{18.2}\]

— on its face a $J$-dimensional integral over the cone where $j$ wins; but only utility differences decide the comparison, and differences of normals are normal, so subtracting $U_{ntj}$ from the others collapses it to a $(J-1)$-dimensional normal orthant probability. That differenced form is the one every computation downstream actually uses. For $J = 2$ it collapses to the univariate $\Phi(\cdot)$ of the binary probit — closed form, essentially. For $J = 3$ it is a bivariate normal probability — computable by good numerical routines (pmvnorm), no longer a formula. For the $J$ of a real supermarket category, it is the reason this model waited until Part IV. Everything hard about the MNP is Equation 18.2; everything in Chapter 19 is a strategy for never quite computing it.

18.3 Identification: The Sternest Version Yet

Section 7.5 established two facts about what choice data can reveal — only utility differences matter, and the scale of utility is arbitrary — and the MNL absorbed both with quiet normalizations (no common intercept; error variance fixed by convention). The MNP has to pay the same two taxes, but on a much larger taxable estate, because $\Omega$’s $J(J+1)/2$ free-looking parameters are exactly where the unidentifiability hides. Walk it through for our $J = 3$ study; the logic is McCulloch and Rossi (1994)’s, and it governs everything in the next two chapters.

Differences. Since only differences matter, subtract a base alternative’s utility — say the Acer’s — from the others. The data can only ever speak about the differenced system: two equations, with differenced attributes and a differenced error vector whose covariance is

\[ \tilde\Omega = M \Omega M', \qquad M = \begin{pmatrix} -1 & 1 & 0 \\ -1 & 0 & 1 \end{pmatrix} \tag{18.3}\]

$\tilde\Omega$ is $2 \times 2$ symmetric: three numbers. The original $\Omega$ had six. Three whole dimensions of $\Omega$ are invisible to choice data — infinitely many $\Omega$’s produce identical behavior. (This is the multivariate version of “a constant added to all utilities cancels”; adding any common noise to all three alternatives changes nothing.)

Scale. Multiply all utilities by $\lambda > 0$: choices unchanged, so one more dimension dies. Convention: fix the first diagonal of the differenced covariance, $\tilde\omega_{11} = 1$. Two bookkeeping notes the next chapter leans on. First, from here on $\tilde\Omega$ means the differenced covariance after this normalization; its unnormalized cousin, raw $M \Omega M'$, will reappear in Chapter 19 under the name $W$ — the object the Gibbs sampler actually draws. Second, the scale tax hits $\boldsymbol{\beta}$ as well: dividing the differenced system by $\sqrt{(M\Omega M')_{11}}$ to enforce $\tilde\omega_{11}=1$ divides every utility, so what the data identify is $\boldsymbol{\beta}/ \sqrt{(M\Omega M')_{11}}$, not $\boldsymbol{\beta}$ itself. Coefficients scale like utility; covariances like utility squared. Count what survives: $3 - 1 = $ two free covariance parameters for our three-alternative study. From six down to two — the MNP’s flexibility is real but far smaller than $\Omega$’s dimensions suggest, and it shrinks further if you forget that each ASC also consumes Section 7.5’s level normalization (our brand-labeled design carries dell and apple ASCs against the Acer base, exactly as in Chapter 10).

Why belabor this? Because unidentified parameters are the classic MNP implementation disaster, in both estimation routes. Hand an optimizer a likelihood written in terms of raw $\Omega$ and it wanders flat directions forever, “converging” wherever it likes (Chapter 6’s warnings, in their natural habitat); run a Gibbs sampler on raw $\Omega$ and the chain drifts happily along the unidentified ridge — traces that never settle, $\widehat{R}$ that never approaches 1 — while the identified functions of the parameters mix beautifully. The two routes of Chapter 19 handle this differently (the MLE route parameterizes only identified quantities; the Bayesian route samples the unidentified space and reports identified ratios, post-processed), and knowing the identified count in advance — $(J-1)J/2 - 1$ covariance parameters — is how you audit either.

18.4 Warming Up: Binary Probit by MLE

For $J = 2$ everything is tractable, and it is worth thirty seconds to close a loop left open since Chapter 3, where the binary probit was one changed line in the simulator. Its choice probability is $\Phi(\mathbf{x}'\boldsymbol{\beta})$ (error scale fixed at 1 — the scale tax, paid), its log-likelihood mirrors Equation 5.5 with pnorm in place of plogis, and Part I’s machinery does the rest:

set.seed(180)
N <- 4000
price <- runif(N, 0.8, 2.4); ram16 <- rbinom(N, 1, 0.5)
X <- cbind(1, price, ram16)
beta_true <- c(1.0, -1.2, 0.6)
y <- as.integer(as.vector(X %*% beta_true) + rnorm(N) > 0)

loglik_bprobit <- function(beta, y, X) {
    P <- pnorm(as.vector(X %*% beta))
    sum(y * log(P) + (1 - y) * log(1 - P))
}
fit <- optim(rep(0, 3), loglik_bprobit, y = y, X = X, method = "BFGS",
             control = list(fnscale = -1), hessian = TRUE)
round(rbind(truth = beta_true, est = fit$par,
            se = sqrt(diag(solve(-fit$hessian)))), 3)

       [,1]   [,2]  [,3]
truth 1.000 -1.200 0.600
est   1.049 -1.220 0.623
se    0.081  0.052 0.045

Recovery, no drama — and compare with Chapter 16’s Albert–Chib sampler on data of exactly this form: two estimation philosophies already agree on the binary case. Everything difficult about the MNP is about $J > 2$: the moment a matrix of error covariances exists, Equation 18.2 stops being pnorm and identification stops being automatic.

18.5 Simulating MNP Data

Simulation, as usual, is the easy direction — the integral of Equation 18.2 never appears, because we run the mechanics (draw errors, add, argmax), not the probabilities. The errors are multivariate normal, so Section 2.4’s Cholesky transformation is the whole trick. We give the truth the OS-correlation story: Acer and Dell errors correlated at $0.6$, Apple independent, unit variances —

\[ \Omega^\ast = \begin{pmatrix} 1 & 0.6 & 0 \\ 0.6 & 1 & 0 \\ 0 & 0 & 1 \end{pmatrix} \tag{18.4}\]

— with the usual attribute coefficients and brand ASCs (Equation 7.6’s values, now interpreted against normal-scale errors²):

sim_mnp_data <- function(N, T, beta, Omega) {
    J <- 3
    n_rows <- N * T * J
    df <- data.frame(
        id    = rep(seq_len(N), each = T * J),
        task  = rep(rep(seq_len(T), each = J), times = N),
        alt   = factor(rep(c("Acer", "Dell", "Apple"), times = N * T),
                       levels = c("Acer", "Dell", "Apple")),
        ram   = factor(sample(c("8","16","32"), n_rows, replace = TRUE),
                       levels = c("8","16","32")),
        screen = factor(sample(c("13","15"), n_rows, replace = TRUE),
                        levels = c("13","15")),
        price = round(runif(n_rows, min = 0.8, max = 2.4), 2)
    )
    X <- model.matrix(~ alt + ram + screen + price, data = df)[, -1]
    V <- as.vector(X %*% beta)

    L   <- t(chol(Omega))                       # @sec-mvnorm
    eps <- as.vector(L %*% matrix(rnorm(n_rows), nrow = J))  # J x (N*T), stacked back
    U   <- V + eps
    task_id <- cumsum(!duplicated(df[, c("id", "task")]))
    df$choice <- as.integer(ave(U, task_id, FUN = function(u) u == max(u)))
    df
}

set.seed(181)
beta_true  <- c(dell = 0.5, apple = 1.0, ram16 = 0.6, ram32 = 0.9,
                screen15 = 0.3, price = -1.2)
Omega_true <- matrix(c(1, 0.6, 0,
                       0.6, 1, 0,
                       0,   0, 1), 3, 3)
mnp <- sim_mnp_data(N = 500, T = 8, beta = beta_true, Omega = Omega_true)

## save for the estimation chapters
saveRDS(list(data = mnp, beta_true = beta_true, Omega_true = Omega_true),
        file = "../data/mnp_study.rds")

One line rewards a careful look: the error draw reshapes the row vector into a $3 \times (NT)$ matrix so that each task’s three errors — one column — get the $L$ transformation as a unit, then stacks back to row order. Correlation lives within tasks, independence across them: the covariance structure of Equation 18.1, implemented in the shape of the data (Chapter 4’s sort contract doing silent work — alt varies fastest, so columns are tasks).

The usual battery, plus one check specific to this model. Marginals and directions first:

chosen <- mnp[mnp$choice == 1, ]
round(prop.table(table(chosen$alt)), 3)


 Acer  Dell Apple 
0.152 0.298 0.550

Apple leads, Acer trails — the ASC ordering, as always. Now the check that matters: is the error correlation visible in behavior? Under IIA-obeying errors, how often a shopper’s second-favorite brand is the other same-OS machine wouldn’t depend on error correlation; under Equation 18.4 it should. A sharper, simpler signature: compare the joint frequency that the two Windows machines both beat Apple against what independent errors would produce. Cheapest of all — lean on the model directly: recompute choice shares after deleting the Dell from every task, using fresh simulations, and watch where its share goes:

set.seed(182)
sim_shares <- function(Omega, drop_dell = FALSE) {
    ## one standardized task: identical specs, $1.5k each
    X_task <- rbind(acer = c(0,0,1,0,1,1.5), dell = c(1,0,1,0,1,1.5),
                    apple = c(0,1,1,0,1,1.5))
    keep <- if (drop_dell) c(1, 3) else 1:3
    V <- as.vector(X_task %*% beta_true)[keep]
    L <- t(chol(Omega[keep, keep, drop = FALSE]))
    R <- 200000
    eps <- L %*% matrix(rnorm(length(keep) * R), nrow = length(keep))
    tabulate(max.col(t(V + eps)), length(keep)) / R
}
before <- sim_shares(Omega_true)
after  <- sim_shares(Omega_true, drop_dell = TRUE)
list(before = round(setNames(before, c("acer","dell","apple")), 3),
     after  = round(setNames(after, c("acer","apple")), 3))

$before
 acer  dell apple 
0.106 0.299 0.595 

$after
 acer apple 
 0.24  0.76

Delete the Dell and its customers flow disproportionately to the Acer — the correlated-error sibling — not proportionally as MNL arithmetic would insist. The red-bus/blue-bus pathology of Chapter 7, cured by two off-diagonal numbers. (Run sim_shares(diag(3)) and the flow turns proportional: the cure really is the correlation, not the normality.)

18.6 The Estimation Problem Ahead

Take stock of the asymmetry we have manufactured. Simulating MNP choices: eight lines, all of them Part I material. Estimating MNP parameters: the likelihood needs Equation 18.2 — an integral we can only approximate — once per task per candidate parameter, inside an optimizer that will make thousands of candidate visits. Chapter 12’s playbook applies but needs an unbiased, smooth simulator of an orthant probability, which is a subtler object than the mixed logit’s average-of-closed-forms; building one (GHK, from Chapter 2’s truncated normals) is half of Chapter 19.

The other half never touches the integral. Chapter 16’s Albert–Chib logic generalizes: conditional on the latent utilities, the MNP is just seemingly-related regressions with normal errors — conjugate territory — and conditional on parameters and the observed choice, each task’s utilities are a multivariate normal truncated to the region where the chosen alternative wins. Draw utilities, draw parameters, repeat: the intractable integral is dissolved into truncated normal draws, at scale, with the identification post-processing that McCulloch and Rossi (1994) made standard. Two routes, one model, next chapter.

18.7 Key Learnings

The MNP replaces Gumbel errors with $N(\mathbf{0}, \Omega)$, $\Omega$ free (Equation 18.1): substitution patterns are estimated in the error covariance rather than specified in nests or proxied by random coefficients. No IIA anywhere — our deletion experiment showed Dell’s share flowing to its correlated sibling, disproportionately.
The choice probability (Equation 18.2) is a $(J-1)$-dimensional normal orthant integral: closed-form at $J=2$ (pnorm), numerical at $J=3$, and the central obstacle beyond. Everything in Chapter 19 is a strategy for living with it.
Identification is unforgiving: differencing kills $J$ of $\Omega$’s dimensions, scale kills one more, leaving $(J-1)J/2 - 1$ free covariance parameters (two, for our study). Unidentified MNP parameterizations produce optimizers that wander and samplers that never converge — count before you code.
Simulation needed none of it: Cholesky-correlated errors within tasks (shaped by the data’s sort contract), argmax, done — plus a saved answer key including $\Omega^\ast$ for the recovery tests ahead.
Two estimation routes await: simulate the integral inside MSL (GHK), or augment it away entirely (truncated normals + Bayesian regression). The capstone comparison of the book’s two estimation philosophies, on one model.

The mixed logit already offered estimated substitution flexibility, and McFadden and Train (2000)’s theorem says it can mimic the MNP arbitrarily well — so why bother? Three honest answers: the MNP puts the flexibility exactly where the story says it lives (the error covariance) rather than asking random coefficients to proxy for it; its Bayesian treatment (Chapter 19) is startlingly clean, with no tuning at all; and its identification lessons are the sternest in choice modeling — rite-of-passage material.↩︎
Alert readers of Section 7.5 will note the scale subtlety: coefficients against unit-normal errors are not numerically comparable to the same values against Gumbel errors. We reuse the familiar numbers for continuity; they simply mean slightly different signal-to-noise here. The recovery checks below compare against the correctly transformed identified quantities, which is the only comparison the data license.↩︎