The math behind choicer: nested logit

This document provides a detailed mathematical description of the nested logit (NL) model as implemented in src/nestlogit.cpp.

Table of Contents

1. Notation
2. Model Definition
3. Log-Likelihood Function
4. Gradient Computation
5. Hessian Computation
   • 4.1 Notation and Scope
   • 4.2 V-Space Hessian Q_i
   • 4.3 Beta-Beta Block
   • 4.4 Beta-Delta and Delta-Delta Blocks
   • 4.5 Lambda-V Cross Derivative
   • 4.6 Lambda-Beta and Lambda-Delta Blocks
   • 4.7 Lambda-Lambda Block
   • 4.8 Single-Pass Implementation
6. Implementation Details
7. Post-Estimation Quantities
8. Choice-Based Sampling and WESML Weighting

---

Notation

Symbol	Description
$i = 1, \ldots, N$	Index for individuals (choice situations)
$j = 1, \ldots, J_i$	Index for inside alternatives available to individual $i$
$C_i$	Choice set of inside alternatives available in situation $i$
$j_i$	The alternative chosen by individual $i$
$w_i$	Weight for individual $i$
$X_{ij}$	Row vector of covariates for individual $i$ and alternative $j$ ($1 \times K$)
$\beta$	Coefficient vector ($K \times 1$)
$\delta_j$	Alternative-specific constant (ASC) for alternative $j$
$V_{ij}$	Systematic (deterministic) utility
$B_k$	Global set of inside alternatives belonging to nest $k$; situation $i$ contributes $B_k\cap C_i$
$k(j)$	The nest containing alternative $j$
$\lambda_k$	Inclusive value (dissimilarity) coefficient for nest $k$
$I_k$	Inclusive value for nest $k$
$P(j \mid k)$	Probability of choosing $j$ conditional on choosing nest $k$
$P_k$	Marginal probability of choosing nest $k$
$P_{ij}$	Probability that individual $i$ chooses alternative $j$

---

1. Model Definition

1.1 Utility Specification

The utility that individual $i$ derives from alternative $j$ is:

$$U_{ij} = V_{ij} + \varepsilon_{ij}$$

where the systematic utility is:

$$V_{ij} = X_{ij}\beta + \delta_j$$

The error term $\varepsilon_{ij}$ follows a Generalized Extreme Value (GEV) distribution that induces within-nest correlation.

1.2 Nest Structure

The $J$ inside alternatives are partitioned into $K$ mutually exclusive and exhaustive nests $B_1, \ldots, B_K$. An outside option (alternative 0) with $V_{i0} = 0$ may optionally be included, and it is always assigned its own singleton nest with $\lambda_0 = 1$.

1.3 Inclusive Value Coefficients

Each nest $k$ has a dissimilarity parameter $\lambda_k > 0$ that governs within-nest substitution:

• $\lambda_k = 1$ corresponds to independent (MNL-like) alternatives within the nest.
• $0 < \lambda_k < 1$ indicates positive within-nest correlation (the typical case).
• $\lambda_k > 1$ is mathematically evaluable in the implemented probability formula but lies outside the standard globally RUM-consistent nested-GEV restriction $0<\lambda_k\leq1$. It should not automatically be interpreted as a valid negative within-nest correlation parameter.

The optimizer in run_nestlogit() currently imposes the lower bound $\lambda_k > 0$ but does not impose $\lambda_k \leq 1$. An estimate above one should therefore be read as a diagnostic about the nesting tree or maintained assumptions, not as evidence of especially strong within-nest substitution.

Singleton nests (nests with exactly one alternative) have $\lambda_k$ fixed to 1 and are not estimated. Only non-singleton nests contribute $\lambda$ parameters to $\theta$. The low-level parser requires at least one globally non-singleton nest; a specification in which every nest is a singleton is rejected rather than estimated as an MNL through this entry point.

Code reference: src/nestlogit.cpp

1.4 Alternative-Specific Constants (ASCs)

The ASC identification follows the same convention as MNL:

• Without outside option (include_outside_option = FALSE): The first inside alternative’s ASC is fixed to zero ($\delta_1 = 0$). The parameter vector contains $J-1$ free ASC parameters.
• With outside option (include_outside_option = TRUE): The outside option serves as the reference with $V_{i0} = 0$. All $J$ inside alternatives have free ASC parameters.

Code reference: src/nestlogit.cpp

---

2. Log-Likelihood Function

2.1 Probability Decomposition

The nested logit choice probability decomposes as:

$$P_{ij} = P(j \mid k(j)) \cdot P_{k(j)}$$

The two components are derived from the GEV generating function.

2.2 Inclusive Value

The inclusive value for nest $k$ is the log-sum of scaled utilities within the nest:

$$I_{ik} = \sum_{j \in B_k\cap C_i} \exp\!\left(\frac{V_{ij}}{\lambda_k}\right)$$

$$\log I_{ik} = \log\!\left(\sum_{j \in B_k\cap C_i} \exp\!\left(\frac{V_{ij}}{\lambda_k}\right)\right)$$

The rest of the article suppresses the $i$ subscript and writes $I_k$. Nests with no available alternative in $C_i$ have log inclusive value $-\infty$ and contribute zero to the upper-level denominator.

2.3 Conditional Choice Probability

The probability of choosing alternative $j$ given that nest $k(j)$ is chosen:

$$P(j \mid k) = \frac{\exp(V_{ij}/\lambda_k)}{I_k}$$

In log form:

$$\log P(j \mid k) = \frac{V_{ij}}{\lambda_k} - \log I_k$$

2.4 Nest (Marginal) Probability

The probability of choosing nest $k$:

$$P_k = \frac{\exp(\lambda_k \log I_k)}{\sum_{l} \exp(\lambda_l \log I_l)}$$

If an outside option is included, its contribution enters the denominator as $\exp(\lambda_0 \log I_0) = \exp(1 \cdot \log(\exp(0/1))) = \exp(0) = 1$.

In log form:

$$\log P_k = \lambda_k \log I_k - \log\!\left(\sum_l \exp(\lambda_l \log I_l)\right)$$

2.5 Joint Choice Probability

Taking logs of the joint probability:

$$\log P_{ij} = \log P(j \mid k(j)) + \log P_{k(j)}$$

$$\log P_{ij} = \frac{V_{ij}}{\lambda_{k(j)}} - \log I_{k(j)} + \lambda_{k(j)} \log I_{k(j)} - \log D_i$$

where $D_i = \sum_l \exp(\lambda_l \log I_l)$ (plus 1 if an outside option is present).

Code reference: src/nestlogit.cpp

2.6 Log-Likelihood

The log-likelihood is the weighted sum of log-probabilities of the observed choices:

$$\ell(\theta) = \sum_{i=1}^{N} w_i \log P_{ij_i}$$

where $\theta = (\beta, \lambda, \delta)$ and $j_i$ is the alternative chosen by individual $i$.

Code reference: src/nestlogit.cpp

2.7 Log-Sum-Exp Trick for Numerical Stability

Two log-sum-exp stabilizations are applied.

Within-nest inclusive value. To compute $\log I_k$ without overflow, the implementation subtracts the maximum scaled utility in each nest:

$$\log I_k = V_{\max,k} + \log\!\left(\sum_{j \in B_k} \exp\!\left(\frac{V_{ij}}{\lambda_k} - V_{\max,k}\right)\right)$$

where $V_{\max,k} = \max_{j \in B_k} V_{ij}/\lambda_k$.

Across-nest denominator. The denominator $D_i = \sum_l \exp(\lambda_l \log I_l)$ is similarly stabilized by subtracting $\max_l \lambda_l \log I_l$ before exponentiation.

Code reference: src/nestlogit.cpp

2.8 Relation to Other Nested-Logit Parameterizations

The probabilities in §§2.3–2.4 are the utility-normalized GEV form used by choicer (McFadden 1978; Train 2009, Ch. 4): utility is scaled by $1/\lambda_k$ inside the conditional probability, and the inclusive value is multiplied by $\lambda_k$ in the nest probability. Some software uses a normalization in which its reported systematic index enters the within-nest logit unscaled. These conventions can represent different definitions of utility coefficients and error scales; they are not generally numerically identical under a common-coefficient restriction, and their reported $\beta$ and $\lambda$ cannot be compared without deriving the normalization mapping. It is too strong to label every unscaled convention non-RUM; the precise point is which utility, inclusive-value, and scale normalization is being estimated (Koppelman & Wen 1998; Hensher & Greene 2002; Heiss 2002).

---

3. Gradient Computation

3.1 Gradient with Respect to $V_{ia}$ (Utility Scores)

To derive the gradients with respect to $\beta$ and $\delta$, we first compute the derivative of the log-probability of the chosen alternative with respect to the systematic utility $V_{ia}$ of each alternative $a$.

Let $k_i$ denote the chosen nest (i.e., $k(j_i)$). The log-probability of the chosen alternative is:

$$\log P_{ij_i} = \log P(j_i \mid k_i) + \log P_{k_i}$$

Derivative of the conditional term. Since $\log P(j \mid k) = V_{ij}/\lambda_k - \log I_k$ and $\frac{\partial \log I_k}{\partial V_{ia}} = \frac{1}{\lambda_k} \mathbf{1}_{a \in B_k} P(a \mid k)$:

$$\frac{\partial \log P(j_i \mid k_i)}{\partial V_{ia}} = \frac{1}{\lambda_{k_i}}\left(\mathbf{1}_{a = j_i} - \mathbf{1}_{a \in B_{k_i}} P(a \mid k_i)\right)$$

Derivative of the nest term. Using $\frac{\partial \log D_i}{\partial V_{ia}} = P_{k(a)} P(a \mid k(a)) = P_{ia}$:

$$\frac{\partial \log P_{k_i}}{\partial V_{ia}} = \mathbf{1}_{a \in B_{k_i}} P(a \mid k_i) - P_{ia}$$

Combined. Adding both terms:

$$\frac{\partial \log P_{ij_i}}{\partial V_{ia}} = \frac{\mathbf{1}_{a = j_i}}{\lambda_{k_i}} + \left(1 - \frac{1}{\lambda_{k_i}}\right)\mathbf{1}_{a \in B_{k_i}} P(a \mid k_i) - P_{ia}$$

This is implemented as the vector grad_vec of length $J_i$:

$$\text{grad\_vec}[a] = -P_{ia} + \left(1 - \frac{1}{\lambda_{k_i}}\right) \mathbf{1}_{a \in B_{k_i}} P(a \mid k_i) + \frac{1}{\lambda_{k_i}} \mathbf{1}_{a = j_i}$$

Code reference: src/nestlogit.cpp

3.2 Gradient with Respect to $\beta$

Since $\frac{\partial V_{ij}}{\partial \beta} = X_{ij}^T$, the contribution from individual $i$ is:

$$\frac{\partial \ell_i}{\partial \beta} = w_i \sum_{j=1}^{J_i} \text{grad\_vec}[j] \cdot X_{ij}^T = w_i \, X_i^T \, \text{grad\_vec}$$

The full gradient is accumulated as a single BLAS matrix-vector product, identically to MNL.

Code reference: src/nestlogit.cpp

3.3 Gradient with Respect to $\delta$

Since $\frac{\partial V_{ij}}{\partial \delta_a} = \mathbf{1}_{j = a}$, the gradient entries are the corresponding elements of grad_vec, scattered to the appropriate positions in $\theta$:

$$\frac{\partial \ell_i}{\partial \delta_a} = w_i \cdot \text{grad\_vec}[a]$$

The same free-ASC convention as MNL applies (first inside alternative’s ASC excluded when no outside option).

Code reference: src/nestlogit.cpp

3.4 Gradient with Respect to $\lambda_k$

Define the within-nest probability-weighted mean utility:

$$\bar{V}_k = \sum_{j \in B_k} P(j \mid k) \cdot V_{ij}$$

For the chosen nest ($k = k_i$). The log-probability decomposes as:

$$\log P_{ij_i} = \underbrace{\frac{V_{ij_i}}{\lambda_{k_i}} - \log I_{k_i}}_{\log P(j_i \mid k_i)} + \underbrace{\lambda_{k_i} \log I_{k_i} - \log D_i}_{\log P_{k_i}}$$

Using $\frac{\partial \log I_k}{\partial \lambda_k} = -\bar{V}_k / \lambda_k^2$:

$$\frac{\partial \log P(j_i \mid k_i)}{\partial \lambda_{k_i}} = \frac{\bar{V}_{k_i} - V_{ij_i}}{\lambda_{k_i}^2}$$

$$\frac{\partial \log P_{k_i}}{\partial \lambda_{k_i}} = (1 - P_{k_i})\!\left(\log I_{k_i} - \frac{\bar{V}_{k_i}}{\lambda_{k_i}}\right)$$

Combining:

$$\frac{\partial \ell_i}{\partial \lambda_{k_i}} = w_i \left[(1 - P_{k_i})\!\left(\log I_{k_i} - \frac{\bar{V}_{k_i}}{\lambda_{k_i}}\right) + \frac{\bar{V}_{k_i} - V_{ij_i}}{\lambda_{k_i}^2}\right]$$

For a non-chosen nest ($k \neq k_i$). The conditional probability $P(j_i \mid k_i)$ is independent of $\lambda_k$, so only the nest probability contributes:

$$\frac{\partial \log P_{k_i}}{\partial \lambda_k} = -P_k \!\left(\log I_k - \frac{\bar{V}_k}{\lambda_k}\right)$$

$$\frac{\partial \ell_i}{\partial \lambda_k} = -w_i \, P_k \!\left(\log I_k - \frac{\bar{V}_k}{\lambda_k}\right)$$

The common factor $\log I_k - \bar{V}_k / \lambda_k$ is labeled term_in_brackets in the code.

Code reference: src/nestlogit.cpp

---

4. Hessian Computation

This section documents the Hessian of the negated log-likelihood. The function nl_loglik_hessian_parallel returns $H(\theta) = -\nabla^2 \ell(\theta)$, the same quantity that invert_hessian() expects (positive semi-definite at the MLE). The sign convention matches the numeric oracle nl_loglik_numeric_hessian (which differences the negated gradient), so both functions are interchangeable as inputs to invert_hessian().

Two implementations are available:

• nl_loglik_hessian_parallel — default (se_method = "hessian"): single-pass analytical Hessian derived below. Validated against the numeric oracle at a maximum absolute difference of $\approx 2.6 \times 10^{-8}$ across the package’s test fixtures, well below the $10^{-6}$ tolerance those tests enforce. At extreme $\lambda_k < 0.1$, the oracle’s own finite-difference accuracy degrades; the analytical Hessian remains finite and symmetric throughout.
• nl_loglik_numeric_hessian — retained oracle (se_method = "numeric"): central finite differences on the analytical gradient with step $\epsilon_i = 10^{-6} \cdot \max(|\theta_i|, 1)$. Use for diagnostic comparison or whenever a direct verification of the analytical path is required.

---

4.1 Notation and Scope

All notation follows §§0–3. Subscripts $i$ (individual), $j/a/b$ (alternative local indices within $i$’s choice set, 0-based), and $k$/$l$ (nest indices, 0-based) are used throughout. Additional symbols specific to this section:

Symbol	Description
$k_i$	0-based index of the chosen nest ($-1$ when outside option is chosen)
$c_a$	0-based local index of the chosen inside alternative
$\bar{V}_k$	$\sum_{j \in B_k} P(j \mid k) \cdot V_{ij}$ — $P(\cdot \mid k)$-weighted mean utility in nest $k$
$T_k$	$\log I_k - \bar{V}_k / \lambda_k$ — the “bracket term” from §3.4
$\widetilde{V}^2_k$	$\sum_{j \in B_k} P(j \mid k) \cdot V_{ij}^2$ — conditional second moment
$\mathrm{VarV}_k$	$\widetilde{V}^2_k - \bar{V}_k^2$ — $P(\cdot\mid k)$-weighted variance of $V$ in nest $k$
$\mathbf{m}_k$	$\sum_{j \in B_k} P_{ij} \cdot \mathbf{x}_{ij}$ — $P_{ij}$-weighted mean of $X$ in nest $k$ ($K$-vector)
$\mathbf{c}_k$	$\sum_{j \in B_k} P(j \mid k) \cdot \mathbf{x}_{ij}$ — $P(\cdot\mid k)$-weighted mean of $X$ in nest $k$ ($K$-vector)
$\mathbf{S}^{(2)}_k$	$\sum_{j \in B_k} P_{ij} \cdot \mathbf{x}_{ij} \mathbf{x}_{ij}^T$ — $P_{ij}$-weighted second moment of $X$ ($K \times K$)
$\mathbf{Q}_k$	$\sum_{j \in B_k} P(j \mid k) \cdot \mathbf{x}_{ij} \mathbf{x}_{ij}^T$ — $P(\cdot\mid k)$-weighted second moment of $X$ ($K \times K$)
$\bar{\mathbf{x}}_i$	$\sum_j P_{ij} \cdot \mathbf{x}_{ij}$ — $P_{ij}$-weighted global mean of $X$ for individual $i$ ($K$-vector)

The negated per-individual log-likelihood contribution is $f_i = -w_i \log P_{i,c_a}$. The full Hessian is the sum $H = \sum_{i=1}^N H^{(i)}$ where $H^{(i)} = -w_i \nabla^2 \log P_{i,c_a}$.

---

4.2 V-Space Hessian $Q_i$

Because $V_{ij} = \mathbf{x}_{ij} \beta + \delta_j$ is linear in $\beta$ and $\delta$, the second derivatives with respect to these parameters factor through the $m_i \times m_i$ V-space Hessian. Define:

$$Q_i[a, b] = -\frac{\partial^2 \log P_{i,c_a}}{\partial V_{ia} \, \partial V_{ib}}$$

Starting from the gradient formula (§3.1):

$$g[a] = \frac{\partial \log P_{i,c_a}}{\partial V_{ia}} = -P_{ia} + \left(1 - \frac{1}{\lambda_{k_i}}\right) \mathbf{1}_{a \in B_{k_i}} P(a \mid k_i) + \frac{\mathbf{1}_{a = c_a}}{\lambda_{k_i}}$$

differentiating with respect to $V_{ib}$ yields:

$$\frac{\partial g[a]}{\partial V_{ib}} = -\frac{\partial P_{ia}}{\partial V_{ib}} + \left(1 - \frac{1}{\lambda_{k_i}}\right) \mathbf{1}_{a \in B_{k_i}} \frac{\partial P(a \mid k_i)}{\partial V_{ib}}$$

Partial derivative of $P_{ia}$ with respect to $V_{ib}$. Write $k_a = k(a)$, $k_b = k(b)$. Expanding $P_{ia} = P(a \mid k_a) \cdot P_{k_a}$ via the product rule and using $\partial P_{k_a}/\partial V_{ib} = P_{k_a}[\mathbf{1}_{k_b=k_a} P(b\mid k_a) - P_{ib}]$:

$$\frac{\partial P_{ia}}{\partial V_{ib}} = P_{ia} \left\{ \mathbf{1}_{k_a = k_b}\left[ \frac{\mathbf{1}_{a=b}}{\lambda_{k_a}} + P(b \mid k_a)\left(1-\frac{1}{\lambda_{k_a}}\right) \right] - P_{ib} \right\}$$

Partial derivative of $P(a \mid k_i)$ with respect to $V_{ib}$. Non-zero only when $a \in B_{k_i}$:

$$\frac{\partial P(a \mid k_i)}{\partial V_{ib}} = \frac{1}{\lambda_{k_i}} \mathbf{1}_{k_b = k_i} P(b \mid k_i) \left(\mathbf{1}_{a=b} - P(a \mid k_i)\right)$$

Combined V-space Hessian (correct formula). Negating and combining:

$$\boxed{ Q_i[a, b] = P_{ia} \left\{ \mathbf{1}_{k_a = k_b}\left[ \frac{\mathbf{1}_{a=b}}{\lambda_{k_a}} +P(b \mid k_a)\left(1-\frac{1}{\lambda_{k_a}}\right) \right] - P_{ib} \right\} - \underbrace{\mathbf{1}_{k_a = k_i} \mathbf{1}_{k_b = k_i} \cdot \frac{1 - 1/\lambda_{k_i}}{\lambda_{k_i}} \cdot P(a \mid k_i)\bigl(\mathbf{1}_{a=b} - P(b \mid k_i)\bigr)}_{C_{a,b}: \text{ chosen-nest correction}} }$$

The chosen-nest correction $C_{a,b}$ is nonzero only when both $a$ and $b$ are in the chosen nest $k_i$ (and $k_i \geq 0$). It arises from the $(1 - 1/\lambda_{k_i})$ factor in $g[a]$ coupling the conditional probability gradient back to $V_{ib}$.

Structural notes and MNL limit:

• For $k_a \neq k_b$: $Q_i[a,b] = -P_{ia} P_{ib}$ (cross-nest entries are the negative product of marginal probabilities, the familiar multinomial covariance term).
• MNL limit ($\lambda_k = 1$ for all $k$): $Q_i[a,b] = P_{ia}(\mathbf{1}_{a=b} - P_{ib})$, the standard multinomial Hessian of $-\log P_{c_a}$. Own diagonal $= P_{ia}(1-P_{ia})$; off-diagonal $= -P_{ia}P_{ib}$. The chosen-nest correction vanishes because $1 - 1/\lambda_{k_i} = 0$.
• When the outside option is chosen ($k_i = -1$) or when $\lambda_{k_i} = 1$ (singleton chosen nest), $C_{a,b} = 0$ for all pairs.

Implementation note: The implementation materializes $Q_i$ as an explicit $m_i \times m_i$ matrix, which is efficient for typical choice set sizes ($J \leq 20$). The outer-product decomposition used for the beta-beta block (§4.3) avoids forming this matrix explicitly for the $K \times K$ accumulation.

Code reference: src/nestlogit.cpp

---

4.3 Beta-Beta Block

Since $\partial V_{ij}/\partial \beta_r = x_{ij,r}$:

$$H^{(i)}_{\beta_r, \beta_s} = w_i \sum_{a} \sum_{b} x_{ia,r} \cdot Q_i[a,b] \cdot x_{ib,s} = w_i \, \mathbf{X}_i^T \mathbf{Q}_i \mathbf{X}_i$$

The matrix product is decomposed into three terms that can be accumulated in a single nest-level pass, exploiting the outer-product structure of $Q_i$:

Term A — universal outer product (subtracted): $$-w_i \, \bar{\mathbf{x}}_i \, \bar{\mathbf{x}}_i^T$$

This arises from the $-P_{ia} P_{ib}$ term universal to all pairs: $\sum_{a,b} x_{ia,r} P_{ia} P_{ib} x_{ib,s} = \bar{x}_{i,r} \bar{x}_{i,s}$.

Term B — per-nest within-nest correction (added), summed over all nests $k$: $$+w_i \sum_{k} \left[ \left(1 - \frac{1}{\lambda_k}\right) \mathbf{m}_k \mathbf{c}_k^T + \frac{1}{\lambda_k} \mathbf{S}^{(2)}_k \right]$$

For singleton nests ($\lambda_k = 1$): $(1 - 1/\lambda_k) = 0$, so the contribution reduces to $\mathbf{S}^{(2)}_k$.

Term C — chosen-nest correction (subtracted), only when $k_i \geq 0$: $$-w_i \cdot \frac{1 - 1/\lambda_{k_i}}{\lambda_{k_i}} \left( \mathbf{Q}_{k_i} - \mathbf{c}_{k_i} \mathbf{c}_{k_i}^T \right)$$

Here $\mathbf{Q}_{k_i} - \mathbf{c}_{k_i}\mathbf{c}_{k_i}^T$ is the $P(\cdot\mid k_i)$-weighted covariance of $X$ within the chosen nest.

Full beta-beta block (Hessian of the negated log-likelihood; positive semi-definite at the MLE): $$H^{(i)}_{\beta\beta} = w_i \left[ -\bar{\mathbf{x}}_i \bar{\mathbf{x}}_i^T + \sum_k \!\left(\!\left(1-\frac{1}{\lambda_k}\right)\mathbf{m}_k \mathbf{c}_k^T + \frac{1}{\lambda_k}\mathbf{S}^{(2)}_k\right) - \mathbf{1}_{k_i \geq 0} \cdot \frac{1-1/\lambda_{k_i}}{\lambda_{k_i}} \!\left(\mathbf{Q}_{k_i} - \mathbf{c}_{k_i}\mathbf{c}_{k_i}^T\right) \right]$$

MNL limit check. At $\lambda_k = 1$ for all $k$: Term B $= \sum_k \mathbf{S}^{(2)}_k = \mathbf{X}_i^T \mathrm{diag}(P_i) \mathbf{X}_i$; Term C vanishes; the formula yields the standard MNL expression $w_i(\mathbf{X}_i^T \mathrm{diag}(P_i) \mathbf{X}_i - \bar{\mathbf{x}}_i \bar{\mathbf{x}}_i^T)$.

Code reference: src/nestlogit.cpp

---

4.4 Beta-Delta and Delta-Delta Blocks

Since $\partial V_{ij}/\partial \delta_a = \mathbf{1}_{j=a}$, these blocks are direct projections of $Q_i$ onto free-ASC alternatives.

Delta-delta block (for free-ASC alternatives $a'$ and $b'$): $$H^{(i)}_{\delta_{a'}, \delta_{b'}} = w_i \cdot Q_i[a', b']$$

Beta-delta block (for covariate index $r$ and free-ASC alternative $b'$): $$H^{(i)}_{\beta_r, \delta_{b'}} = w_i \sum_{a} x_{ia,r} \cdot Q_i[a, b'] = w_i \, \bigl(\mathbf{X}_i^T \mathbf{Q}_i[:, b']\bigr)_r$$

Both blocks reuse the same $Q_i$ matrix computed in §4.2. No additional auxiliary quantities are required.

ASC normalization. Without an outside option, the first inside alternative’s ASC is fixed to zero: free-ASC alternatives satisfy alt_idx0[a] > 0, mapped to theta positions delta_start_idx + alt_idx0[a] - 1.

Code reference: src/nestlogit.cpp

---

4.5 Lambda-V Cross Derivative

The lambda parameters enter nonlinearly through both $\log I_k$ and the nest probabilities, so cross-derivatives between $\lambda_k$ and $\beta/\delta$ must be computed via the chain rule. Define:

$$d^{(k)}_b = \frac{\partial g_{\lambda_k}}{\partial V_{ib}}, \qquad g_{\lambda_k} = \frac{\partial \log P_{i,c_a}}{\partial \lambda_k}$$

Key intermediates (derived from §3.4):

$$\frac{\partial T_k}{\partial V_{ib}} = -\mathbf{1}_{k(b)=k} P(b \mid k) \frac{V_{ib} - \bar{V}_k}{\lambda_k^2}, \qquad \frac{\partial P_k}{\partial V_{ib}} = P_k \bigl[\mathbf{1}_{k(b)=k} P(b \mid k) - P_{ib}\bigr]$$

For non-chosen nests ($k \neq k_i$), using $g_{\lambda_k} = -P_k T_k$: $$d^{(k)}_b = P_k T_k P_{ib} + \mathbf{1}_{k(b)=k} P_k P(b \mid k) \left[\frac{V_{ib} - \bar{V}_k}{\lambda_k^2} - T_k\right]$$

The first term is a universal “shift” proportional to $P_{ib}$; the second is the within-nest correction active only for $b \in B_k$.

For the chosen nest ($k = k_i$, only when $k_i \geq 0$), using $g_{\lambda_{k_i}} = (1-P_{k_i})T_{k_i} + (\bar{V}_{k_i} - V_{c_a})/\lambda_{k_i}^2$: $$d^{(k_i)}_b = P_{k_i} T_{k_i} P_{ib} - \frac{\mathbf{1}_{b = c_a}}{\lambda_{k_i}^2} + \mathbf{1}_{k(b)=k_i} P(b \mid k_i) \left[ -P_{k_i} T_{k_i} + \frac{V_{ib} - \bar{V}_{k_i}}{\lambda_{k_i}^2} \left(P_{k_i} - 1 + \frac{1}{\lambda_{k_i}}\right) + \frac{1}{\lambda_{k_i}^2} \right]$$

The extra $-1/\lambda_{k_i}^2$ term for $b = c_a$ arises from differentiating the $(\bar{V}_{k_i} - V_{c_a})/\lambda_{k_i}^2$ component of the chosen-nest gradient.

---

4.6 Lambda-Beta and Lambda-Delta Blocks

Using $\partial V_{ib}/\partial \beta_r = x_{ib,r}$ and $\partial V_{ib}/\partial \delta_{b'} = \mathbf{1}_{b=b'}$:

$$H^{(i)}_{\lambda_k, \beta_r} = -w_i \sum_{b} d^{(k)}_b \cdot x_{ib,r} = -w_i \, \bigl(\mathbf{X}_i^T \mathbf{d}^{(k)}\bigr)_r$$

$$H^{(i)}_{\lambda_k, \delta_{b'}} = -w_i \, d^{(k)}_{b'}$$

Both are filled symmetrically. Singleton nests (theta index $= -1$) contribute no rows or columns to these blocks.

Code reference: src/nestlogit.cpp

---

4.7 Lambda-Lambda Block

Differentiating $g_{\lambda_k}$ with respect to $\lambda_l$, using: $$\frac{\partial T_k}{\partial \lambda_k} = \frac{\mathrm{VarV}_k}{\lambda_k^3} \geq 0, \qquad \frac{\partial P_k}{\partial \lambda_k} = P_k (1 - P_k) T_k, \qquad \frac{\partial P_k}{\partial \lambda_l} = -P_k P_l T_l \;(l \neq k), \qquad \frac{\partial T_k}{\partial \lambda_l} = 0 \;(l \neq k)$$

Off-diagonal $k \neq l$ (any combination of chosen/non-chosen nests): $$H^{(i)}_{\lambda_k, \lambda_l} = -w_i \, P_k P_l T_k T_l$$

This formula holds uniformly regardless of whether $k$ or $l$ equals $k_i$. It can be negative (whenever $T_k$ and $T_l$ share the same sign), reflecting correlation between nest-level inclusive values.

Diagonal $k \neq k_i$ (non-chosen nest): $$H^{(i)}_{\lambda_k, \lambda_k} = w_i \left[ P_k (1 - P_k) T_k^2 + \frac{P_k \, \mathrm{VarV}_k}{\lambda_k^3} \right]$$

The first term is a Fisher-information-type squared-gradient term; the second is positive and captures curvature from $T_k$’s dependence on $\lambda_k$ through the within-nest utility variance.

Diagonal $k = k_i$ (chosen nest, only when $k_i \geq 0$): $$H^{(i)}_{\lambda_{k_i}, \lambda_{k_i}} = w_i \left[ P_{k_i}(1 - P_{k_i}) T_{k_i}^2 - \frac{(1 - P_{k_i}) \mathrm{VarV}_{k_i}}{\lambda_{k_i}^3} + \frac{\mathrm{VarV}_{k_i}}{\lambda_{k_i}^4} + \frac{2(\bar{V}_{k_i} - V_{c_a})}{\lambda_{k_i}^3} \right]$$

The chosen-nest diagonal differs from the non-chosen formula because the $(\bar{V}_{k_i} - V_{c_a})/\lambda_{k_i}^2$ term in $g_{\lambda_{k_i}}$ contributes $-\mathrm{VarV}_{k_i}/\lambda_{k_i}^4 - 2(\bar{V}_{k_i} - V_{c_a})/\lambda_{k_i}^3$ when differentiated. Since the chosen alternative tends to have above-average utility within its nest ($V_{c_a} \geq \bar{V}_{k_i}$ in expectation), the last term is typically negative, making the chosen-nest diagonal somewhat smaller than the non-chosen counterpart.

When the outside option is chosen ($k_i = -1$), all lambda diagonals use the non-chosen formula.

Code reference: src/nestlogit.cpp

---

4.8 Single-Pass Implementation and Auxiliary Accumulations

All six blocks are accumulated in a single parallel pass over individuals. The per-nest auxiliary quantities below are accumulated in an inner loop over alternatives at $O(m_i)$ cost before any block is computed.

Per-nest scalars (zero-initialized per individual, indexed by nest $k$):

Quantity	Formula	Blocks that use it
$\bar{V}_k$	$\sum_{j \in B_k} P(j\mid k) V_j$	$T_k$, lambda-lambda, lambda-V cross
$\widetilde{V}^2_k$	$\sum_{j \in B_k} P(j\mid k) V_j^2$	$\mathrm{VarV}_k$ for lambda-lambda diagonal
$T_k$	$\log I_k - \bar{V}_k/\lambda_k$	lambda-lambda (all entries), lambda-V cross
$\mathrm{VarV}_k$	$\widetilde{V}^2_k - \bar{V}_k^2$	lambda-lambda diagonal

Per-nest K-vectors:

Quantity	Formula	Blocks that use it
$\mathbf{m}_k$	$\sum_{j \in B_k} P_{ij} \mathbf{x}_j$	beta-beta Term B
$\mathbf{c}_k$	$\sum_{j \in B_k} P(j\mid k) \mathbf{x}_j$	beta-beta Terms B and C

Per-nest $K \times K$ matrices:

Quantity	Formula	Blocks that use it
$\mathbf{S}^{(2)}_k$	$\sum_{j \in B_k} P_{ij} \mathbf{x}_j \mathbf{x}_j^T$	beta-beta Term B
$\mathbf{Q}_k$	$\sum_{j \in B_k} P(j\mid k) \mathbf{x}_j \mathbf{x}_j^T$	beta-beta Term C (chosen nest $k_i$ only)

Global per-individual: $\bar{\mathbf{x}}_i = \sum_j P_{ij} \mathbf{x}_j = \sum_k \mathbf{m}_k$ (beta-beta Term A).

OpenMP structure. The outer loop over individuals is parallelized with dynamic scheduling. Each thread maintains a thread-local $P \times P$ accumulator local_H, merged into global_H via #pragma omp critical after the individual loop. The nest-to-theta-index map is copied thread-privately to avoid false sharing.

Symmetrization and sanitization. After the parallel loop: $H \leftarrow (H + H^T)/2$ enforces exact symmetry; any remaining non-finite entries are zeroed.

Singleton-nest handling. Singleton nests ($\lambda_k = 1$, theta index $= -1$) contribute nothing to the lambda-lambda, lambda-beta, or lambda-delta blocks. Their beta-beta contribution via Term B simplifies to $\mathbf{S}^{(2)}_k$ (since $1 - 1/\lambda_k = 0$).

Per-individual weights. Every contribution is multiplied by $w_i$ (WESML weight; equals 1 for unweighted data).

Code reference: src/nestlogit.cpp

---

5. Implementation Details

5.1 Parameter Vector Structure

The full parameter vector $\theta$ is organized as:

Block	Indices	Length	Description
$\beta$	$[0, K)$	$K$	Coefficients for design matrix $X$
$\lambda$	$[K, K + K_\lambda)$	$K_\lambda$	Non-singleton nest dissimilarity parameters
$\delta$	$[K + K_\lambda, K + K_\lambda + J_{\text{asc}})$	$0$, $J-1$, or $J$	Free alternative-specific constants

where $K_\lambda$ is the number of non-singleton nests and $J_{\text{asc}}=0$ when use_asc = FALSE, $J-1$ without an outside option, or $J$ with an outside option.

Code reference: src/nestlogit.cpp

5.2 Singleton Nest Handling

At startup, the implementation identifies which nests are singletons (contain exactly one alternative). These nests: - Have their $\lambda$ fixed to 1 (no contribution to the upper-level logit denominator beyond their single utility term). - Are excluded from the estimated parameter vector. - Are tracked via nest_k_to_theta_idx, an integer vector of length $K$ where entry $k$ holds the corresponding index in $\theta$ for non-singleton nests and $-1$ for singletons.

This ensures that the NL model degenerates gracefully to MNL on a per-nest basis when nests are singletons.

Code reference: src/nestlogit.cpp

5.3 Data Organization

• Design matrix $X$: Stacked matrix of dimension $(\sum_i M_i) \times K$, where $M_i$ is the number of inside alternatives for individual $i$.
• alt_idx: 1-based indices mapping rows of $X$ to global alternative IDs.
• nest_idx: 1-based nest assignment of length $J$ (number of unique inside alternatives), not length $\sum M_i$. Accessed via nest_idx.elem(alt_idx0_i) to handle varying choice sets.
• choice_idx: 1-based local row position of the chosen inside alternative within each situation’s block; 0 for the outside option.
• Prefix sums $S$: Used for efficient slicing into the stacked data: $S_i = \sum_{j < i} M_j$.

Code reference: src/nestlogit.cpp

5.4 Base Utility Pre-computation

Prior to the parallel loop, all systematic utilities $V_{ij} = X_{ij}\beta + \delta_j$ are computed in a single pass:

base_util = X * beta                 // single BLAS dgemm
if use_asc: base_util += delta[alt_idx0]

This avoids redundant matrix-vector products inside the per-individual loop.

Code reference: src/nestlogit.cpp

5.5 OpenMP Parallelization

The implementation parallelizes over individuals using OpenMP: - Each thread maintains local accumulators for log-likelihood and gradient. - nest_k_to_theta_idx is copied thread-privately to avoid false sharing. - Thread results are combined using #pragma omp critical sections. - Dynamic scheduling balances load across individuals with varying choice set sizes.

Code reference: src/nestlogit.cpp

5.6 Negated Objectives

The likelihood/gradient kernel returns: - $-\ell(\theta)$ (negated log-likelihood) - $-\nabla\ell(\theta)$ (negated gradient)

for compatibility with minimization routines (e.g., nloptr) that minimize a loss function rather than maximize a likelihood. The analytical and numeric Hessian kernels return $-\nabla^2\ell$; the BHHH kernel instead returns the positive-semidefinite OPG information approximation $\sum_iw_is_is_i^T$.

5.7 Finite-Value Guards

The likelihood kernel replaces a non-finite chosen log probability by $-10^{10}$. If the final negated objective is non-finite it returns $10^{10}$ and a zero gradient; isolated non-finite gradient entries are zeroed. The analytical Hessian skips a situation whose chosen log probability is non-finite, symmetrizes the aggregate, and replaces any remaining non-finite entries by zero. These rules are optimizer failure guards, not part of the regular finite-parameter nested-logit criterion; encountering them calls for rescaling or revising the specification.

Code reference: src/nestlogit.cpp

---

6. Post-Estimation Quantities

6.1 Choice Probabilities (predict)

The joint probability is the decomposition from §2.1:

$$P_{ij} = P(j \mid k(j)) \cdot P_{k(j)}$$

with $P(j \mid k)$ and $P_k$ given in §§2.3–2.4. nl_predict returns individual-level probabilities $P_{ij}$ and systematic utilities $V_{ij}$ for every observation row. nl_predict_shares aggregates over individuals to obtain weighted market shares:

$$\hat{s}_j = \frac{\sum_{i=1}^{N} w_i \, P_{ij}}{\sum_{i=1}^{N} w_i}$$

Code reference: src/nestlogit.cpp

---

6.2 Elasticities

6.2.1 Semi-elasticity with respect to utility

Let $r = k(j)$ be the nest of the responding alternative $j$ and $s = k(a)$ be the nest of the perturbed alternative $a$. From §3.1, the derivative of $\log P_{ij}$ with respect to the systematic utility of any alternative $a$ is:

$$\frac{\partial \log P_{ij}}{\partial V_{ia}} = \frac{\mathbf{1}_{a=j}}{\lambda_r} + \left(1 - \frac{1}{\lambda_r}\right)\mathbf{1}_{s=r} \, P(a \mid r) - P_{ia}$$

This is the general expression; the gradient formula in §3.1 applies it only to the chosen alternative, whereas here we evaluate it for all pairs $(j, a)$ to construct the full elasticity matrix.

6.2.2 Elasticity of $P_{ij}$ with respect to covariate $k$ of alternative $a$

Since $\partial V_{ia} / \partial x_{ia,k} = \beta_k$:

$$E_{ij,ia}^{k} = \frac{\partial P_{ij}}{\partial x_{ia,k}} \cdot \frac{x_{ia,k}}{P_{ij}} = \beta_k \, x_{ia,k} \cdot \frac{\partial \log P_{ij}}{\partial V_{ia}}$$

Substituting the three cases:

Case	Condition	Semi-elasticity $\partial \log P_{ij}/\partial V_{ia}$
Own	$a = j$	$\dfrac{1}{\lambda_r} + \left(1 - \dfrac{1}{\lambda_r}\right) P(j \mid r) - P_{ij}$
Cross, same nest	$s = r,\; a \neq j$	$\left(1 - \dfrac{1}{\lambda_r}\right) P(a \mid r) - P_{ia}$
Cross, different nest	$s \neq r$	$-P_{ia}$

Reduction to MNL. When $\lambda_k = 1$ for all nests the within-nest term $(1 - 1/\lambda_r)$ vanishes and the conditional probability drops out, giving $\partial \log P_{ij}/\partial V_{ia} = \mathbf{1}_{a=j} - P_{ia}$, which is exactly the MNL semi-elasticity (§3.1 of the MNL document).

Outside option. When an outside option is present it occupies its own singleton nest with $\lambda_0 = 1$ and $x_{i0,k} = 0$. Its row (responding = outside option) has semi-elasticity $-P_{ia}$ for every inside alternative $a$ (different nest, §3.1); its column (perturbed = outside option) is identically zero because $\beta_k \cdot 0 = 0$.

The implementation accumulates $w_i \cdot E_{ij,ia}^k$ across individuals and divides by $\sum_i w_i$ to produce a weighted average of individual probability elasticities. This is not generally the elasticity of an aggregate share under a common market-level perturbation when attributes vary across situations. Rows index responding alternatives and columns perturbed alternatives.

Code reference: src/nestlogit.cpp

---

6.3 Diversion Ratios

6.3.1 Definition

The utility-based diversion ratio from losing alternative $j$ to receiving alternative $k$ is:

$$DR(j \to k) = -\frac{\partial Q_k / \partial V_{ij}}{\partial Q_j / \partial V_{ij}}$$

where $Q_m = \sum_i w_i P_{im}$ is aggregate demand and the perturbation is through the systematic utility of the losing alternative $j$.

6.3.2 Partial derivatives of demand

Using $\partial P_{im}/\partial V_{ij} = P_{im} \cdot (\partial \log P_{im}/\partial V_{ij})$ and the semi-elasticity from §6.2.1, the three cases with $r = k(j)$:

$$\frac{\partial P_{im}}{\partial V_{ij}} = \begin{cases} P_{ij}\!\left[\dfrac{1}{\lambda_r} + \left(1 - \dfrac{1}{\lambda_r}\right) P(j \mid r) - P_{ij}\right] & m = j \\[6pt] P_{im}\!\left[\left(1 - \dfrac{1}{\lambda_r}\right) P(j \mid r) - P_{ij}\right] & k(m) = r,\; m \neq j \\[6pt] -P_{im} \, P_{ij} & k(m) \neq r \end{cases}$$

6.3.3 Utility Perturbation and Covariate Interpretation

The implemented matrix is defined with respect to a common additive perturbation of $V_{ij}$. If a fixed-coefficient attribute of alternative $j$ receives the same additive change $dx$ in every situation, $dV_{ij}=\beta_r dx$, so the common scalar cancels:

$$DR(j \to k) = -\frac{\beta_r dx\sum_i w_i \,(\partial P_{ik}/\partial V_{ij})} {\beta_r dx\sum_i w_i \,(\partial P_{ij}/\partial V_{ij})} = -\frac{\sum_i w_i \,(\partial P_{ik}/\partial V_{ij})} {\sum_i w_i \,(\partial P_{ij}/\partial V_{ij})}$$

It is in this additive-shift sense that NL diversion is independent of the chosen fixed covariate. For a proportional perturbation $dx_i=x_{ij,r}d\log x$, the values $x_{ij,r}$ remain inside both population sums and generally do not cancel. The kernel does not compute that attribute-specific proportional diversion. In MXL even a common additive perturbation of a random-coefficient attribute carries draw-specific $\beta_{ir}^s$, so the coefficient cannot be factored out.

6.3.4 Population aggregation and column-sum property

Define:

$$\text{num}(k, j) = \sum_{i=1}^{N} w_i \left(-\frac{\partial P_{ik}}{\partial V_{ij}}\right), \qquad \text{den}(j) = \sum_{i=1}^{N} w_i \frac{\partial P_{ij}}{\partial V_{ij}}$$

so that $DR(j \to k) = \text{num}(k,j) \,/\, \text{den}(j)$ for $k \neq j$.

Columns sum to 1. The fundamental identity $\sum_{m} \partial P_{im}/\partial V_{ij} = 0$ (probabilities sum to 1 for every individual and every perturbation direction, including the outside option when present) gives:

$$\sum_{k \neq j} \frac{\partial P_{ik}}{\partial V_{ij}} = -\frac{\partial P_{ij}}{\partial V_{ij}} \implies \sum_{k \neq j} \text{num}(k, j) = \text{den}(j) \implies \sum_{k \neq j} DR(j \to k) = 1$$

When an outside option is present it acts as both a source (it can lose demand) and a destination (it can receive diverted demand), so it appears in both the row index $k$ and the column index $j$ of the $J_{\text{total}} \times J_{\text{total}}$ diversion matrix, and the column sums still equal 1 over all alternatives including the outside option. Because its systematic utility is fixed at zero in estimation, the outside-as-source column is a derivative with respect to a hypothetical perturbation of that normalized utility, not an available outside-option covariate.

6.3.5 Reduction to MNL at $\lambda = 1$

When all $\lambda_k = 1$ the cross-same-nest case collapses to the cross-different-nest case ($\partial P_{im}/\partial V_{ij} = -P_{im} P_{ij}$ for $m \neq j$) and the own case gives $\partial P_{ij}/\partial V_{ij} = P_{ij}(1 - P_{ij})$. Substituting:

$$DR(j \to k)\big|_{\lambda=1} = \frac{\sum_i w_i \, P_{ij} \, P_{ik}}{\sum_i w_i \, P_{ij} (1 - P_{ij})}$$

which is the MNL diversion ratio (§6.3 of the MNL document).

Code reference: src/nestlogit.cpp

---

6.4 BLP Share Inversion

6.4.1 Problem statement

Given observed market shares $s_j^*$, find the ASC vector $\delta$ such that the NL model reproduces those shares:

$$\hat{s}_j(\delta) = s_j^* \quad \forall j$$

6.4.2 Damped contraction mapping

Under the usual feasibility, connected-substitutes, and normalization conditions, the BLP fixed-point iteration is

$$\delta^{(t+1)} = \delta^{(t)} + \kappa \bigl(\log s^* - \log \hat{s}(\delta^{(t)})\bigr)$$

with $\kappa=1$ giving the undamped update. Existence, uniqueness, and global convergence are properties of the share mapping and target—not an unconditional guarantee for every numerical input. Predicted shares $\hat{s}(\delta^{(t)})$ are computed via the NL formula:

$$\hat{s}_j(\delta^{(t)}) = \frac{\sum_{i=1}^{N} w_i \, P_{ij}(\delta^{(t)})}{\sum_{i=1}^{N} w_i}$$

where $P_{ij}(\delta^{(t)}) = P(j \mid k(j)) \cdot P_{k(j)}$ is evaluated at the current $\delta^{(t)}$ with $\beta$ and $\lambda$ held fixed.

6.4.3 Damping for nested models

For strongly nested data (small $\lambda_k$) the plain update can overshoot numerically. Choosing $0<\kappa<1$ may stabilize iteration. A possible empirical starting heuristic is

$$\kappa \approx \min_k \lambda_k$$

but this is not a theorem or an automatic rule in the package. The kernel accepts the numeric damping value as supplied and does not constrain it to $(0,1]$; the analyst must inspect convergence and resulting share residuals.

6.4.4 Baseline re-normalization

Without an outside option, delta is a full length-$J$ vector and is centered initially and on return by subtracting its first inside entry, so the returned baseline is zero. With an outside option, supplied and returned delta contain the $J$ inside ASCs only. The kernel works internally with a prepended outside entry, pins that entry at zero after each update, and returns the inside subvector without subtracting the first inside ASC.

All target shares must be strictly positive and have length $J$ (no outside) or $J+1$ (outside first). The kernel does not check that they sum to one. Convergence uses the largest absolute change in the full working vector; reaching max_iter prints a warning and returns the last iterate without a separate convergence flag.

Note on the lambda argument. nl_blp_contraction expects a full length-$n_{\text{nests}}$ vector where singleton nests have their entry set to 1. This differs from the estimated parameter block $\lambda$ in $\theta$, which contains only non-singleton entries. The caller is responsible for expanding the estimated lambdas into the full vector before invoking the contraction.

Code reference: src/nestlogit.cpp

6.5 Willingness to Pay

Utility is linear in covariates, so willingness to pay is the same marginal rate of substitution as in the MNL model: for attribute $k$ and price coefficient $\alpha = \beta_p$,

$$\mathrm{WTP}_k = -\frac{\beta_k}{\alpha}, \qquad \widehat{\mathrm{SE}} = \sqrt{\nabla g^T \, \hat{V}_{(k,p)} \, \nabla g}, \quad \nabla g = \begin{pmatrix} -1/\alpha \\ \beta_k/\alpha^2 \end{pmatrix}$$

with $\hat{V}_{(k,p)}$ the $2 \times 2$ block of the coefficient covariance. The nest parameters $\lambda$ rescale within-nest substitution but do not enter the ratio of marginal utilities, so they affect WTP only through the estimated covariance of $(\beta_k, \alpha)$. ASCs may also be used as the numerator ($-\delta_j/\alpha$). See §8 of the MNL document for the full derivation and caveats.

Code reference: R/wtp.R

6.6 Consumer Surplus and the Nested Logsum

For the nested logit, expected maximum utility is the two-level logsum below plus the scenario-invariant extreme-value location constant (omitted by logsum(), as in MNL):

$$\mathrm{logsum}_i = \log D_i = \log\left( \sum_{l} \exp\big(\lambda_l \log I_l\big) \; [+\, 1] \right), \qquad \log I_l = \log \sum_{j \in B_l \cap C_i} \exp\!\left(\frac{V_{ij}}{\lambda_l}\right)$$

where the $+1 = \exp(0)$ term appears when an outside option is present and singleton nests have $\lambda_l = 1$ (§5.2). Both levels are computed with the max-subtraction trick (§2.7).

Two useful identities:

• Reduction to MNL. With all $\lambda_l = 1$, $\lambda_l \log I_l = \log \sum_{j \in B_l} \exp(V_{ij})$ and the nested logsum collapses to the MNL logsum $\log \sum_j \exp(V_{ij})$.
• Outside-option link. Since $P_{i0} = 1 / D_i$ (§2.5), the logsum equals $-\log P_{i0}$ — a direct consistency check against the probability kernel.

Expected consumer surplus follows Train (2009, Ch. 3) exactly as in the MNL case:

$$\mathbb{E}[CS_i] = \frac{\mathrm{logsum}_i}{-\alpha}$$

with the same caveats: no income effects, and CS levels inherit the ASC normalization, so only differences across scenarios (counterfactual newdata vs. baseline) are meaningful. The implementation reports point estimates for NL (the delta-method SE of the mean CS is currently provided for MNL only; Krinsky–Robb resampling of the coefficients is the practical interval method here). The routine does not enforce $\alpha<0$; interpreting $-\alpha$ as positive marginal utility of income requires the analyst to verify the estimated sign.

Code reference: R/surplus.R

6.7 Goodness of Fit

The measures are model-agnostic and identical to §10 of the MNL document, computed from the NL log-likelihood and predicted probabilities:

• McFadden pseudo R²: $R^2 = 1 - \ell(\hat\theta)/\ell_0$ and $R^2_{\text{adj}} = 1 - (\ell(\hat\theta) - K)/\ell_0$, with $K$ counting all estimated parameters including the $\lambda$ block.
• Equal-shares null (default, exact for unbalanced sets and weights): $\ell_0 = -\sum_i w_i \log(M_i + \mathbf{1}_{\text{outside}})$.
• Market-shares null (constants-only closed form): $\ell_0 = c\sum_j N_j \log(s_j)$, where $c$ is the common weight; requires identical choice-set composition across individuals and uniform weights. The constants-only model is itself an MNL—at $\beta = 0$ and free ASCs the nest structure is not identified from constant shares—so the same closed form applies.
• Hit rate: weighted share of choice situations whose observed choice has the highest predicted NL probability; with an outside option the outside good competes with $P_{i0} = 1 - \sum_j P_{ij}$.

Code reference: R/gof.R

---

7. Choice-Based Sampling and WESML Weighting

7.1 Endogenous Stratified (Choice-Based) Sampling

Under random (exogenous) sampling each choice situation is drawn independently of its outcome. Under choice-based (endogenous stratified) sampling the strata are defined by the chosen alternative $j_i$ and situations are drawn with stratum-specific frequencies — for example, deliberately oversampling individuals who chose a rare alternative. Let $Q(j)$ be the population share of alternative $j$ and $H(j)$ the sample share of choosers of $j$.

Unlike the MNL special case (a full set of ASCs leaves the slopes consistent), the nested logit has no such free pass: the dissimilarity parameters $\lambda_l$ enter the within- and between-nest probabilities nonlinearly, so choice-based sampling generally biases the slopes $\beta$, the constants $\delta$, and the $\lambda$ block. Weighting is therefore required for consistent estimation of the population parameters.

7.2 The WESML Weighted Log-Likelihood

Manski and Lerman (1977) restore consistency by weighting each situation’s contribution by the ratio of population to sample share of its chosen alternative,

$$w_i = \frac{Q(j_i)}{H(j_i)}, \qquad \ell^{W}(\theta) = \sum_{i=1}^{N} w_i \log P_{i j_i}(\theta).$$

The maximizer $\hat\theta$ is invariant to multiplying every $w_i$ by a common positive constant, so the weights may be normalized to mean 1 without changing the estimates.

7.3 Robust (Sandwich) Variance: the $w^2$ Meat

WESML is a weighted M-estimator solving $\sum_i w_i\, s_i(\hat\theta) = 0$, where $s_i = \partial \log P_{i j_i}/\partial\theta$ is the full per-situation score over the $\beta$, $\lambda$ and $\delta$ blocks ($$3) and $H_i$ its Hessian ($$4). Its robust (Huber–White) asymptotic variance is the sandwich

$$V = A^{-1} B A^{-1}, \qquad A = \sum_{i} w_i\,(-H_i), \qquad B = \sum_{i} w_i^{2}\, s_i s_i^{\top}.$$

The weight enters the bread $A$ linearly, but the meat $B$ carries the weight squared: the contribution of situation $i$ to the estimating equation is $w_i s_i$, so its variance involves $(w_i s_i)(w_i s_i)^{\top} = w_i^{2}\, s_i s_i^{\top}$. The score $s_i$ is assembled exactly as in the gradient kernel ($$3.3–3.4): the $\beta$ and $\delta$ blocks from grad_vec, and the $\lambda_l$ blocks from the per-nest expression $(\,\mathbf{1}\{l = l_i\} - P_{il}\,)\big[\log I_{il} - \lambda_l^{-1}\!\sum_{j\in B_l} P_{ij|l} V_{ij}\big] + \mathbf{1}\{l=l_i\}\,\lambda_l^{-2}\big(\sum_{j\in B_l} P_{ij|l}V_{ij} - V_{ij_i}\big)$, with singleton nests ($\lambda_l \equiv 1$) contributing no score.

As in the MNL/MXL cases, both naive variances are wrong under weighting: the inverse-Hessian $A^{-1}$ assumes the information-matrix equality $A = B$ (which fails once the $w_i$ are non-degenerate), and ordinary BHHH/OPG $\big(\sum_i w_i\, s_i s_i^{\top}\big)^{-1}$ uses the weight to the first power, not the second. The sandwich is invariant to the weight scale: $A \to cA$, $B \to c^2 B$, hence $V \to A^{-1} B A^{-1}$.

7.4 Scope: Robust vs. Design-Based Variance

The estimator implemented here is the robust weighted-M-estimator (Huber–White) variance, with uncentered meat $B = \sum_i w_i^{2}\, s_i s_i^{\top}$ — the asymptotic variance under variable-probability sampling. The smaller stratum-centered (design-based) variance for a fixed-quota design is out of scope; see Manski & McFadden (1981) and Cosslett (1981).

7.5 Implementation and Cluster-Robust Extension

The common post-estimation path obtains the $N\times p$ matrix $S$ of weight-free scores from nl_scores_parallel (singleton lambda rows are absent because those parameters are not in $\theta$), then forms

$$A=\texttt{nl\_loglik\_hessian\_parallel}(\widehat\theta;w), \qquad B=(\operatorname{diag}(w)S)^T(\operatorname{diag}(w)S).$$

For a fit whose selected method is the numeric NL Hessian, the post-hoc bread uses the numeric oracle instead. The eager fit-time robust path is algebraically identical and may obtain $B$ by calling the BHHH accumulator with weights = w^2. For cluster labels $g(i)$, se_method = "cluster" uses

$$B_{\mathrm{cluster}}=\sum_g \left(\sum_{i:g(i)=g}w_is_i\right) \left(\sum_{i:g(i)=g}w_is_i\right)^T,$$

with no finite-cluster correction. These matrices are combined as $V=A^{-1}BA^{-1}$. User-facing entry points: wesml_weights() computes the Manski–Lerman weights from the population shares $Q$, while sample_by_choice() draws a choice-based sample and attaches those weights; run_nestlogit(..., se_method = "sandwich") estimates with the robust variance; wesml_vcov() returns it post hoc. Passing run_nestlogit(..., weights_col = ) collapses a row-level weight column to one weight per choice situation (validated constant within id), and a provenance guard prevents a WESML-labeled sample from being silently fit unweighted.

Code reference: nl_scores_parallel() / nl_bhhh_parallel() in src/nestlogit.cpp; compute_sandwich_vcov(), .score_meat(), and .assemble_score_vcov() in R/classes.R; wesml_weights() / sample_by_choice() in R/sampling.R.

---

References

• McFadden, D. (1978). Modelling the choice of residential location. In A. Karlqvist et al. (Eds.), Spatial Interaction Theory and Planning Models (pp. 75-96). North-Holland.
• Ben-Akiva, M., & Lerman, S. R. (1985). Discrete Choice Analysis: Theory and Application to Travel Demand. MIT Press.
• Koppelman, F. S., & Wen, C.-H. (1998). Alternative nested logit models: structure, properties and estimation. Transportation Research Part B, 32(5), 289-298.
• Hensher, D. A., & Greene, W. H. (2002). Specification and estimation of the nested logit model: alternative normalisations. Transportation Research Part B, 36(1), 1-17.
• Heiss, F. (2002). Structural choice analysis with nested logit models. The Stata Journal, 2(3), 227-252.
• Berry, S., Levinsohn, J., & Pakes, A. (1995). Automobile prices in market equilibrium. Econometrica, 63(4), 841–890.
• Train, K. E. (2009). Discrete Choice Methods with Simulation (2nd ed.). Cambridge University Press. Chapter 4.
• Manski, C. F., & Lerman, S. R. (1977). The estimation of choice probabilities from choice based samples. Econometrica, 45(8), 1977–1988.
• Manski, C. F., & McFadden, D. (1981). Alternative estimators and sample designs for discrete choice analysis. In C. F. Manski & D. McFadden (Eds.), Structural Analysis of Discrete Data with Econometric Applications (pp. 2–50). MIT Press.
• Cosslett, S. R. (1981). Maximum likelihood estimator for choice-based samples. Econometrica, 49(5), 1289–1316.