Sequential Estimation of Dynamic Discrete Choice Models with Unobserved Heterogeneity

Sequential Estimation of Dynamic Discrete Choice Models with Unobserved Heterogeneity

\pagenumbering{arabic} \setcounter{page}{1}

\onehalfspacing

Introduction

Dynamic discrete choice (DDC) models are widely used to study forward-looking behavior. In many of these applications, unobserved heterogeneity is empirically central keane1997career,francesconi2002joint,gowrisankaran2012dynamics,wang2015impact. At the same time, incorporating unobserved heterogeneity is computationally costly because it requires repeatedly solving fixed-point equations for all unobserved types. This burden is substantial in single-agent models with large state spaces, and it is magnified further in dynamic games aguirregabiria2010dynamic,arcidiacono2011practical. A natural question is whether one can reduce the cost of these inner fixed-point solves without sacrificing statistical efficiency.

Existing approaches reduce this burden in different ways. Full-solution methods solve the Bellman equation at each parameter evaluation rust1987optimal, while Conditional Choice Probability (CCP)-based estimators replace the dynamic programming problem with fixed-point problems in CCP spaces hotz1993conditional,aguirregabiria2002swapping,aguirregabiria2007sequential,arcidiacono2011conditional. Mathematical Programming with Equilibrium Constraints approach su2012constrained,dube2012improving treats the fixed-point problem as a constraint in a joint optimization problem. Other approaches exploit special structure of the model, such as finite dependence arcidiacono2011conditional,aguirregabiria2023solving, and index invertibility bunting2025faster. bugni2021iterated study truncation of outer iterations. adusumilli2025temporal develop Temporal Difference Estimation, which uses a sieve to approximate the value function. dearing2025efficient develop Efficient Pseudo-Likelihood (EPL), which is efficient and improves the convergence properties of the sequential estimator in dynamic games. These methods substantially expand the range of feasible applications, but incorporating unobserved heterogeneity remains computationally demanding. In particular, it is unclear how truncating the inner fixed-point solver affects the statistical and computational properties of the resulting estimator.

This paper develops a unified framework for sequential estimation of DDC models with unobserved heterogeneity. We propose the EM-NPL($q$) estimator, which combines the sequential pseudo-likelihood framework of aguirregabiria2002swapping,aguirregabiria2007sequential with the finite-mixture Expectation-Maximization (EM) approach of arcidiacono2011conditional. At each outer iteration, the estimator applies only $q$ steps of a chosen inner algorithm to approximate the fixed-point object required in the M-step, and then updates the structural parameters and mixture weights. The framework is modular: it accommodates different fixed-point equations including Bellman equation rust1987optimal, policy valuation hotz1993conditional,aguirregabiria2007sequential, Euler equation aguirregabiria2023solving, and EPL dearing2025efficient. It can be paired with different numerical algorithms such as successive approximation, Generalized Minimal Residual (GMRES), or Newton-type methods.

For the workhorse class of linear-in-parameters DDC models, we show EM-NPL($q$) is numerically identical to the EM-NPL estimator that solves the inner fixed-point problem to convergence, for any $q \geq 1$. Therefore, the number of inner iterations and the choice of inner solver only affect computation, but not statistical properties. This result is practically important because it decouples numerical implementation from asymptotic properties in the models most often used in empirical work.

We also establish consistency and asymptotic normality of the estimator, characterize local convergence of the EM-NPL($q$) algorithm, and derive an iteration-count bound that links the number of outer iterations to the approximation error of the truncated inner solver. In particular, we show that if the EM-NPL algorithm converges locally, then EM-NPL($q$) also converges locally for sufficiently large $q$. Together, these results clarify when inner truncation changes only runtime and when it may also affect convergence behavior.

We provide practical implementation guidance, which is organized around three stages. First, benchmark candidate inner solvers on the initial fixed-point problem and discard dominated algorithms. Second, select the fixed-point equation and inner algorithm based on three model features: the functional form of the utility function (linear-in-parameters models and linear systems are suited to GMRES), whether finite dependence holds (which makes the Euler equation method aguirregabiria2023solving available), and the discount factor (which governs how quickly successive approximation converges). Third, choose $q$ to balance inner accuracy against the number of outer iterations: for superlinear/quadratic solvers such as GMRES/Newton, $q = 4$ to $6$ typically suffices, whereas successive approximation may require $q = 8$ or more when $\beta$ is large. In linear-in-parameters models, the choice of $q$ is computational, not statistical, and can be guided by runtime considerations alone.

We apply EM-NPL($q$) to two Monte Carlo simulations: a single-agent entry--exit model with unobserved heterogeneity and a dynamic game. In the single-agent model, truncated implementations reduce runtime by about 20--40% relative to exact inner convergence and are often several times faster than standard alternatives. In dynamic games, policy-valuation GMRES is the fastest method. The EPL estimator is about 10 times slower than the policy-valuation approach, while it has the smaller mean squared error as it is efficient. However, our method speeds up EPL estimator by about 40%. The simulations therefore provide concrete guidance on how to choose the fixed-point equation, inner solver, and truncation level.

We then apply EM-NPL($q$) to household demand for cola using Worldpanel by Numerator Take Home data. The model has a state space of 546{,}875 points, and we estimate up to ten latent types separately for three income groups, making the application a demanding test of computational feasibility. The computational time is approximately linear in the number of types, and it takes about 20 minutes to estimate a model with ten types. The results show that ignoring unobserved heterogeneity can understate long-run own-price elasticities by up to about 60%, short-run own-price elasticities by up to about 85%, and compensating variation from a soda tax by about 90%. The application shows that accounting for unobserved heterogeneity is not only substantively important but also computationally feasible.

Notation. All limits below are taken as $N \rightarrow \infty$, unless stated otherwise. Let $:=$ denote “equals by definition.” Let $\mathbbm{1}\{A\}$ denote an indicator function that takes the value of one when $A$ is true and zero otherwise.

Outline. (ref) describes the finite-mixture dynamic discrete choice model. (ref) describes the EM-NPL($q$) estimator. (ref) presents the theoretical properties. (ref) provides implementation guidance. (ref) provides Monte Carlo simulations. (ref) presents the cola-demand application. (ref) concludes. All proofs are in Appendix (ref).

Model

We consider a stationary dynamic discrete choice (DDC) model with unobserved heterogeneity, following kasahara2009nonparametric and arcidiacono2011conditional. Heterogeneity is modeled as a finite mixture with $M$ unobserved types.

A panel of $N$ markets is observed, each with $J$ firms indexed by $j = 1, \ldots, J$ and observed for $T$ periods. The unobserved heterogeneity operates at the market level: market $i$'s type $m_i \in \{1,\ldots,M\}$ is drawn once at the start and remains fixed throughout the sample. Let $\pi_m = \Pr(m_i = m)$ denote the mixing weight, with $\pi_m > 0$ for all $m$ and $\sum_{m=1}^{M} \pi_m = 1$. In what follows, we suppress the market index $i$.

Time is discrete and infinite, indexed by $t = 0, 1, 2, \ldots$. Denote the action space by $a \in \mathcal{A} = \{0, 1, \ldots, |\mathcal{A}|-1\}$ where action 0 denotes the outside option. In each period $t$, given a vector of state variables $x_{t} \in \mathcal{X}$ observable to all firms, and firm $j$'s private information $\varepsilon_{jt}:= (\varepsilon_{jt}(a))_{a \in \mathcal{A}} \in \mathbb{R}^{|\mathcal{A}|}$, firm $j$ chooses an action $a_{jt} \in \mathcal{A}$ simultaneously with other firms. The transition of the state variables is controlled by a first-order Markov process with transition probability $f(x_{t+1}, \varepsilon_{t+1}|x_{t}, a_{t}, \varepsilon_{t}; \theta_{f})$, common across types, where $a_{t} = (a_{1t}, \ldots, a_{Jt})$ is the vector of actions taken by all firms in period $t$, $\varepsilon_{t} = (\varepsilon_{1t}, \ldots, \varepsilon_{Jt})$ and $\theta_{f} \in \Theta_{f}$ is a vector of finite-dimensional parameters.

Each firm has type-specific flow utility $U_{j}(x_{t}, a_{t}, \varepsilon_{jt}; \theta^{m})$ where $\theta^{m} \in \Theta \subseteq \mathbb{R}^{d}$ is a vector of structural parameters, and chooses its action to maximize the expected discounted utility,

where $\beta \in (0, 1)$ is the discount factor. When $M = 1$, this reduces to the standard single-type DDC model. We impose the following standard assumptions.

The transition density is common across types, and $\theta_{f}$ can be estimated separately from $\theta^m$. Thus, we focus on the estimation of $\theta^{m}$ and suppress the dependence on $\theta_{f}$.

We consider Markov Perfect Equilibrium where each firm's strategy depends only on the current state variables. We only consider stationary equilibria and suppress the time subscript $t$. For each type $m$, let $\sigma^{m} := (\sigma_{1}^{m}, \ldots, \sigma_{J}^{m})$ be a strategy profile where $\sigma_{j}^{m}: \mathcal{X} \times \mathbb{R}^{|\mathcal{A}|} \to \mathcal{A}$ maps the current state variables and private information to an action. The Conditional Choice Probability (CCP) for firm $j$ in market of type $m$ is

For a market of type $m$ and firm $j$, let $\sigma_{-j}^{m}$ be the strategy profile of other firms, and $P_{-j}^{\sigma^{m}}$ be the corresponding CCP profile. The expected utility and transition density for firm $j$ are

Under the extreme value assumption on private values, the integrated value function for firm $j$ given $\sigma^{m}_{-j}$ is the solution to the Bellman equation:

where $\kappa$ is the Euler constant. The conditional value function for firm $j$ is

Then, CCPs take the logit form:

Framework

Following the sequential estimation approach of aguirregabiria2002swapping,aguirregabiria2007sequential and the Expectation--Maximization (EM) framework of arcidiacono2011conditional, we define the EM-NPL($q$) estimator as a fixed point of a system of maximum pseudo-likelihood and fixed-point equations. The conventional approach solves the fixed-point equations exactly; our estimator replaces the exact solver with $q$ iterations of an inner algorithm. We first introduce the fixed-point equations that arise in estimation, then define the EM-NPL($q$) estimator and the algorithm used to compute it.

Fixed-Point Equations for Estimation

In this section, we suppress the type index $m$ for notational simplicity. Let $P := (P_{1}, \ldots, P_{J})$ stack the CCPs across firms. At each maximization step, given $(\tilde{\theta}, P)$ from the previous iteration, we need to solve for the nuisance parameter $Y$ (e.g., the value function) that satisfies the fixed-point equation:

Throughout the paper, we assume (ref) has a unique fixed point for any given $(\theta, \tilde{\theta}, P)$. The mapping $G$ can only depend on the current candidate $\theta$ as in the NPL estimator of aguirregabiria2007sequential, or it may also depend on the previous $\tilde{\theta}$ as in the Efficient Pseudo-Likelihood (EPL) estimator of dearing2025efficient. Different choices of $G$ and $Y$ give rise to different estimators; we present examples below.

The EM-NPL($q$) Estimator

Following aguirregabiria2007sequential and arcidiacono2011conditional, we define the EM-NPL($q$) estimator as a fixed point of a system of equations that maximizes the pseudo-likelihood.

For each firm $j$, the CCP mapping $\Lambda_j: \Theta \times \mathcal{Y} \times \mathcal{P} \to \Delta(\mathcal{A})^{|\mathcal{X}|}$ maps structural parameters, the nuisance parameter, and the CCP profile to CCPs for firm $j$ using the logit formula. We write $\Lambda = (\Lambda_1, \ldots, \Lambda_J)$ for the stacked mapping across all firms. Let $\Gamma^{q}(\theta, \tilde{\theta}, P; Y)$ denote a generic inner algorithm that applies $q$ steps to approximate the solution of (ref), starting from the initial guess $Y^m$ while holding $(\theta, \tilde{\theta}, P)$ fixed. Let $\boldsymbol{\theta} = (\theta^1, \ldots, \theta^M)$, $\mathbf{Y} = (Y^1, \ldots, Y^M)$, $\mathbf{P} = (P^1, \ldots, P^M)$, and $\boldsymbol{\pi} = (\pi_1, \ldots, \pi_M)$. With a slight abuse of notation, we write $\Gamma^q(\boldsymbol{\theta}, \tilde{\boldsymbol{\theta}}, \mathbf{P}; \mathbf{Y})$ and $\Lambda(\boldsymbol{\theta}, \mathbf{Y}, \mathbf{P})$ for the type-by-type stacked mappings, i.e., $\Gamma^q(\boldsymbol{\theta}, \tilde{\boldsymbol{\theta}}, \mathbf{P}; \mathbf{Y}) := \bigl(\Gamma^q(\theta^1, \tilde{\theta}^1, P^1; Y^1), \ldots, \Gamma^q(\theta^M, \tilde{\theta}^M, P^M; Y^M)\bigr)$ and analogously for $\Lambda$. Let $\mathcal{Q}_N^q(\boldsymbol{\theta}, \boldsymbol{\pi};\, \tilde{\boldsymbol{\theta}}, \mathbf{Y}, \mathbf{P})$ denote the pseudo-likelihood function with $(\tilde{\boldsymbol{\theta}}, \mathbf{Y}, \mathbf{P})$ from the previous iteration and $\Gamma^q$ in place of the exact solution to the fixed-point equation:

The EM-NPL($q$) estimator is defined as the element of the fixed-point set $\mathcal{E}_N$ that maximizes the pseudo-likelihood.

The fixed-point conditions in (ref) require: (i) $(\boldsymbol{\theta}, \boldsymbol{\pi})$ maximize the pseudo-likelihood given $(\boldsymbol{\theta}, \mathbf{Y}, \mathbf{P})$ where $\boldsymbol{\theta}$ is the previous iterate of the structural parameters; (ii) $\mathbf{Y}$ is a fixed point of the inner algorithm $\Gamma^q$; and (iii) $\mathbf{P}$ is also a fixed point of the mapping $\Lambda$.

To compute the EM-NPL($q$) estimator, we define the EM-NPL($q$) algorithm following aguirregabiria2007sequential and arcidiacono2011conditional.\footnote{Algorithm (ref) requires initial values $(\boldsymbol{\pi}^{(0)}, \boldsymbol{\theta}^{(0)}, \mathbf{Y}^{(0)}, \mathbf{P}^{(0)})$. See (ref).}

Upon convergence, the output belongs to the fixed-point set $\mathcal{E}_N$ in (ref). (ref) establishes local convergence of the algorithm to the EM-NPL($q$) estimator given a consistent initial estimator.

Inner Algorithms

Any iterative method for solving the fixed-point equation (ref) can serve as the inner algorithm $\Gamma$. We present three examples; for a comprehensive textbook, see judd1998numerical.

Theoretical Properties

This section develops the theoretical properties of the EM-NPL($q$) estimator. We first establish the truncation invariance property, present consistency and asymptotic normality, and finally characterize local convergence of the EM-NPL($q$) algorithm.

Truncation Invariance

Recall that $q$ is the number of inner iterations per outer step. A natural concern is whether truncating the inner solver at $q < \infty$ distorts the estimator. We show that for the workhorse class of DDC models with linear-in-parameters utility, it does not. That is, the EM-NPL($q$) estimator is numerically identical to the EM-NPL estimator for any $q \geq 1$. The key observation is that when the utility function is linear in parameters, the policy valuation ((ref)) and EPL ((ref)) fixed-point equations can be written in a form that does not depend on current $\theta$. We call this $\theta$-separability.

Under $\theta$-separability, $Y$ can be solved once for given $(\tilde{\theta}, P)$ from the previous iteration, rather than re-solved for each candidate $\theta$ during the M-step. It implies that an EM-NPL($q$) fixed-point is also an EM-NPL($\infty$) fixed point. Therefore, the EM-NPL($q$) estimator is numerically identical to the EM-NPL estimator. We further impose the following condition on the inner algorithm and the fixed-point mapping:

(ref)(i) is mild and ensures that the nuisance parameter is pinned down once the structural parameters and CCPs are known. For Examples (ref), (ref), and (ref), the operator $G$ is a contraction, so the fixed point is unique. (ref)(ii)(iii) is mild and can be verified for common algorithms. For example, when $G$ is a contraction, SA satisfies the condition. For linear systems, GMRES also satisfies the condition (see saad2003iterative). The EM-NPL estimator solves the fixed point exactly, i.e., $q=\infty$. The following theorem shows that the EM-NPL($q$) estimator is numerically identical to the EM-NPL estimator for any $q \geq 1$.

(ref) has two important consequences for inference. First, there is no efficiency loss from inner truncation. Standard NPL standard errors are valid without any adjustment for the inner truncation. Second, the choice of $q$ and the inner algorithm $\Gamma$ become purely computational decisions: they affect the speed of convergence of the EM-NPL($q$) algorithm, but not the statistical properties of the estimator. This decoupling between computation and inference is what makes the EM-NPL($q$) framework practical: the researcher can set $q$ to minimize the overall runtime without worrying about the statistical consequences.

Large-Sample Properties

We establish consistency and asymptotic normality of the EM-NPL($q$) estimator by adapting aguirregabiria2007sequential Proposition 2 and arcidiacono2011conditional Theorem 2. Define the population pseudo-likelihood:

The population EM-NPL($q$) fixed point satisfies the following system of equations:

We impose the following identification conditions.

(ref) is the identification assumption, which can be established by kasahara2009nonparametric. Throughout the paper, we assume that the number of types $M$ is known. In practice, one may apply the sequential testing procedure of kasahara2014non. (ref) is similar to the isolated fixed point condition in aguirregabiria2007sequential Proposition 2.

Next, we impose regularity conditions to establish consistency and asymptotic normality. These conditions are similar to those in aguirregabiria2007sequential.

(ref) is similar to the assumptions in aguirregabiria2007sequential Proposition 2. Let $\boldsymbol{\alpha} = (\boldsymbol{\theta}, \boldsymbol{\pi}) \in \mathbb{R}^{d_\alpha}$ where $d_\alpha = M \cdot d_\theta + M - 1$ collect all structural parameters and mixing probabilities. Write the per-observation log-likelihood contribution as

and the score $s_i^q(\boldsymbol{\alpha};\, \tilde{\boldsymbol{\theta}}, \mathbf{Y}, \mathbf{P}) := \nabla_{\boldsymbol{\alpha}} \ell_i^q(\boldsymbol{\alpha};\, \tilde{\boldsymbol{\theta}}, \mathbf{Y}, \mathbf{P})$. Define the information matrix $\Omega_{\boldsymbol{\alpha}\boldsymbol{\alpha}}^{q} = -\nabla_{\boldsymbol{\alpha}\boldsymbol{\alpha}}^2 \mathcal{Q}^q(\boldsymbol{\alpha}^*;\, \boldsymbol{\theta}^*, \mathbf{Y}^*, \mathbf{P}^*)$.

Under Assumptions (ref) and (ref), the implicit function theorem applies. Let $\mathbf{Y}(\boldsymbol{\alpha})$ and $\mathbf{P}(\boldsymbol{\alpha})$ denote the solutions to $Y^m = \Gamma^q(\theta^m, \theta^m, P^m; Y^m)$ and $P^m = \Lambda(\theta^m, Y^m, P^m)$ for all $m$ in a neighborhood of $\boldsymbol{\alpha}^* = (\boldsymbol{\theta}^*, \boldsymbol{\pi}^*)$, and let $\mathbf{Y}_{\boldsymbol{\alpha}}$ and $\mathbf{P}_{\boldsymbol{\alpha}}$ denote the total derivatives of $\mathbf{Y}(\boldsymbol{\alpha})$ and $\mathbf{P}(\boldsymbol{\alpha})$ with respect to $\boldsymbol{\alpha}$. Since $\boldsymbol{\pi}$ does not enter these constraints, $\mathbf{Y}_{\boldsymbol{\pi}} = 0$ and $\mathbf{P}_{\boldsymbol{\pi}} = 0$. Define the correction term $ S^{q} := \Omega^{q}_{\boldsymbol{\alpha}\tilde{\boldsymbol{\theta}}}\,(I_{Md_\theta},\, 0) + \Omega^{q}_{\boldsymbol{\alpha} \mathbf{Y}}\, \mathbf{Y}_{\boldsymbol{\alpha}} + \Omega^{q}_{\boldsymbol{\alpha} \mathbf{P}}\, \mathbf{P}_{\boldsymbol{\alpha}}$ where $\Omega^{q}_{\boldsymbol{\alpha}\tilde{\boldsymbol{\theta}}} := -\nabla_{\boldsymbol{\alpha}\tilde{\boldsymbol{\theta}}}^2 \mathcal{Q}^q$, $\Omega^{q}_{\boldsymbol{\alpha} \mathbf{Y}} := -\nabla_{\boldsymbol{\alpha} \mathbf{Y}}^2 \mathcal{Q}^q$, and $\Omega^{q}_{\boldsymbol{\alpha} \mathbf{P}} := -\nabla_{\boldsymbol{\alpha} \mathbf{P}}^2 \mathcal{Q}^q$, $(I_{Md_\theta},\, 0)$ is the $Md_\theta \times d_\alpha$ selector extracting the $\boldsymbol{\theta}$-block from $\boldsymbol{\alpha}$, and $\mathbf{Y}_{\boldsymbol{\alpha}}$, $\mathbf{P}_{\boldsymbol{\alpha}}$ are the stacked total-derivative matrices with $\mathbf{Y}_{\boldsymbol{\alpha}} = \mathbf{Y}_{\boldsymbol{\theta}}\, (I_{Md_\theta},\, 0)$ and $\mathbf{P}_{\boldsymbol{\alpha}} = \mathbf{P}_{\boldsymbol{\theta}}\, (I_{Md_\theta},\, 0)$.

The term $\Omega^{q}_{\boldsymbol{\alpha}\tilde{\boldsymbol{\theta}}}$ arises because the fixed-point mapping $G$ may depend on the previous-iteration parameter $\tilde{\boldsymbol{\theta}}$, as in the EPL estimator ((ref)). When $G$ depends only on the current $\boldsymbol{\theta}$, these terms vanish and the correction simplifies to $S^q = \Omega^q_{\boldsymbol{\alpha} \mathbf{Y}}\, \mathbf{Y}_{\boldsymbol{\alpha}} + \Omega^q_{\boldsymbol{\alpha} \mathbf{P}}\, \mathbf{P}_{\boldsymbol{\alpha}}$.

(ref) has two regimes for inference. Under (ref), (ref) implies that the EM-NPL($q$) estimator is numerically identical to the EM-NPL estimator for any $q \geq 1$, so $\Sigma^q$ reduces to the standard NPL asymptotic variance and standard errors can be computed using the usual NPL formulas. When (ref) does not hold, $\Sigma^q$ depends on the cross-derivative matrices $\Omega^q_{\boldsymbol{\alpha}\tilde{\boldsymbol{\theta}}}$, $\Omega^q_{\boldsymbol{\alpha}\mathbf{Y}}$, and $\Omega^q_{\boldsymbol{\alpha}\mathbf{P}}$ together, which can be cumbersome to compute in practice. In that case, nonparametric bootstrap may be applied.

To compare the asymptotic variance of EM-NPL($q$) with that of EM-NPL, we impose two additional conditions on the approximation error of the inner algorithm.

(ref) requires that the approximation error is upper bounded by a function $f(q)$, and the quality of the initial guess $Y$. For example, if SA is used for the inner algorithm, then $f(q) = \rho^q$ where $\rho$ is the Lipschitz constant of the operator $G$.

(ref) shows that the asymptotic variance of the EM-NPL($q$) estimator converges to the asymptotic variance of the EM-NPL estimator at the same rate $f(q)$ that governs the inner-algorithm approximation error. For example, if SA is used for the inner algorithm, then $f(q) = \rho^q$ where $\rho$ is the Lipschitz constant of the operator $G$. For superlinear algorithms like GMRES, $f(q)$ can decay much faster than SA. Therefore, we can achieve a good approximation to the EM-NPL asymptotic variance with a small $q$ when using superlinear algorithms.

Local Convergence

We now characterize the local convergence properties of the EM-NPL($q$) algorithm.

The rate matrix $R^q$ has a block-triangular structure that reflects the sequential nature of the EM-NPL($q$) algorithm. The first row corresponds to the E-step: the mixing weights $\boldsymbol{\pi}$ depend on the previous CCPs and themselves (through $E_{\mathbf{P}}$ and $E_{\boldsymbol{\pi}}$). The second row corresponds to the M-step: the sensitivity of $\boldsymbol{\theta}$ to the previous iterates $(\tilde{\boldsymbol{\theta}}, \mathbf{Y}, \mathbf{P}, \boldsymbol{\pi})$, governed by the implicit function $H$ blocks. The third and fourth rows compose the M-step output with the inner algorithm $\Gamma^q$ and the CCP mapping $\Lambda$.

The inner truncation level $q$ enters the rate matrix through $\Gamma_{\mathbf{Y}}^q, \Gamma_{\boldsymbol{\theta}}^q, \Gamma_{\boldsymbol{\tilde{\theta}}}^q$ and $\Gamma_{\mathbf{P}}^q$. When ($q = \infty$), $\Gamma_{\mathbf{Y}}^{\infty} = 0$ because the fixed point equation is solved exactly. For finite $q$, the magnitude $\|\Gamma_{\mathbf{Y}}^q\|, \|\Gamma_{\boldsymbol{\theta}}^q - \Gamma_{\boldsymbol{\theta}}^{\infty}\|, \|\Gamma_{\boldsymbol{\tilde{\theta}}}^q - \Gamma_{\boldsymbol{\tilde{\theta}}}^{\infty}\|, \|\Gamma_{\mathbf{P}}^q - \Gamma_{\mathbf{P}}^{\infty}\|$ is controlled by the approximation error of the inner algorithm. This is the mechanism through which inner truncation affects the convergence speed of the outer loop. Under (ref), $\Gamma_{\boldsymbol{\theta}}^q = 0$ because the fixed-point mapping does not depend on the current structural parameter. When, in addition, $G$ does not depend on $\tilde{\boldsymbol{\theta}}$ (as in the policy valuation mapping of (ref)), $\Gamma_{\boldsymbol{\tilde{\theta}}}^q = 0$ as well.

The EM-NPL($q$) algorithm converges locally whenever $\rho(R^q) < 1$. When the EM-NPL($\infty$) algorithm converges locally, i.e., $\rho(R^{\infty}) < 1$, the approximation error of the inner algorithm induces a perturbation to the rate matrix $R^q$. If the perturbation is small enough, then $\rho(R^q)$ remains less than 1 and the EM-NPL($q$) algorithm also converges locally. Assumptions (ref) and (ref), stated above, formalize this perturbation argument.

(ref) shows that if the EM-NPL algorithm converges locally, then the EM-NPL($q$) algorithm also converges locally for all sufficiently large $q$. The key observation is that the approximation error of the inner algorithm induces a perturbation to the matrix $R^q$ of the EM-NPL($q$) algorithm. Next we characterize how the spectral radius $\rho(R^q)$ converges to $\rho(R^{\infty})$ as $q \to \infty$, and how the number of iterations $\bar{K}_q$ to achieve $\varepsilon$-accuracy scales with $q$.

To translate the perturbation bound on $\rho(R^q)$ into an iteration-count comparison, we specialize to the $\theta$-separable case ((ref)), in which (ref) guarantees that the EM-NPL($q$) and EM-NPL algorithms share a common fixed point $x_* := x_*^q = x_*^\infty$. Starting from a common initial guess $x_0$, let $\{x_k^q\}$ and $\{x_k^\infty\}$ denote the sequences generated by EM-NPL($q$) and EM-NPL, respectively; both converge to $x_*$. Under the local convergence results,

For any $\varepsilon < \|x_0 - x_*\|_2$, define the linear upper bounds on the number of outer iterations needed to achieve $\varepsilon$-accuracy:

where we assume $\rho(R^q),\, \rho(R^\infty) \in (0, 1)$.

When $g(q)$ is small, a first-order Taylor expansion around $g = 0$ gives $\log(\rho(R^\infty) \pm g(q)) = \log\rho(R^\infty) \pm g(q)/\rho(R^\infty) + O(g(q)^2)$, so

(ref) therefore implies the leading-order expansion

where the last step uses $|\log\rho(R^\infty)| \ge 1 - \rho(R^\infty)$ for $\rho(R^\infty) \in (0, 1)$. Two practical implications follow. First, the iteration-count gap $\bar K_q - \bar K_\infty$ shrinks at rate $f(q)^{1/\nu}$, inherited from the inner-solver approximation error. For superlinear inner solvers (GMRES, Newton's method) this rate is exponential in $q$, whereas for successive approximation it closes at rate $\rho_G^{q/\nu}$ where $\rho_G < 1$ is the Lipschitz constant of $G$. Second, the gap is amplified by $1/(1 - \rho(R^\infty))$: when the exact EM-NPL converges slowly, a larger $q$ is required before EM-NPL($q$) matches its outer iteration count.

Practical Guidance

This section provides practical guidance. We organize the practical discussion around the three stages of a typical implementation: benchmarking inner solvers, selecting the algorithm--equation pair, and choosing $q$.

Algorithm and Fixed-Point Equation Selection

The choice of fixed-point equation $G$, inner algorithm $\Gamma$, and truncation level $q$ jointly determines both the statistical and computational performance of the EM-NPL($q$) estimator.

Stage 1: Benchmark Inner Solvers and Eliminate Dominated Algorithms

Every EM-NPL($q$) implementation begins by computing an initial CCP estimate $\hat{P}_{0}$. Given $\hat{P}_{0}$, the first outer iteration requires solving the fixed-point equation for the initial consistent estimator, $Y_{0}$. This step is common to all EM-NPL($q$) variants and provides a natural opportunity to benchmark candidate inner solvers at negligible additional cost.

The practitioner should solve this initial system using each candidate algorithm, e.g., successive approximation (SA), GMRES, etc., and record the wall-clock time. Algorithms that are slower on the practitioner's specific problem can be immediately discarded as dominated. This benchmark is informative because the numerical structure of the initial solver. For example, for the policy valuation equation, the condition number of $(I - \beta F_{P})$, its sparsity pattern, etc., which govern the relative performance of different algorithms, are often preserved across subsequent EM-NPL iterations.

Stage 2: Select the Fixed-Point Equation and Algorithm

Among the non-dominated algorithms from Stage 1, the practitioner next selects the combination of fixed-point equation $G$ and inner algorithm $\Gamma$ best suited to the model. This choice depends on three model features: the functional form of the utility function, the finite-dependence structure, and the discount factor $\beta$.

\paragraph{Functional form.} When the utility function is linear in parameters, $U(x,a) = \phi(x,a)'\theta$, the policy valuation problem ((ref)) and the EPL problem ((ref)) both reduce to a set of linear systems that do not depend on $\theta$. This has two advantages. First, the system needs to be solved only once per outer iteration rather than at every parameter evaluation. Second, (ref) guarantees that the EM-NPL($q$) estimator is equivalent to the EM-NPL estimator for any $q \geq 1$. This means that for single-agent DDCs with linear utility, EM-NPL($q$) is also efficient (aguirregabiria2002swapping). For dynamic games, if we solve the EPL problem ((ref)) and the utility is linear, then our EM-NPL($q$) estimator is also efficient. In this case, the choice of algorithm and $q$ can be guided solely by computational considerations, without concern for statistical efficiency. (ref) summarizes the recommended algorithm--equation pairs based on these model features.

When the utility function is nonlinear in parameters ((ref)), the inner problem is nonlinear and must be solved at every parameter evaluation. In this case, we recommend using Newton's method, as it typically converges much faster than SA. The choice between Bellman and Euler depends on the model structure and ease of implementation.

\paragraph{Finite dependence.} If the model exhibits finite dependence ((ref)), the Euler equation representation is available. The Euler fixed-point mapping has a contraction modulus strictly less than $\beta$ aguirregabiria2023solving, which makes even SA converge relatively quickly.

\paragraph{Discount factor.} The discount factor governs the difficulty of the inner fixed-point problem. When $\beta$ is small (e.g., $\beta = 0.8$), the choice of algorithm matters little. When $\beta$ is close to unity (e.g., $\beta = 0.9999$), SA becomes extremely slow.

Stage 3: Choose $q$

The final choice is the number of inner iterations $q$. Our recommendation is to start small and increase if necessary. For algorithms with superlinear convergence (GMRES, Newton), $q = 4$ to $6$ is typically sufficient for the outer EM-NPL($q$) loop to converge. For EPL, we recommend $q = 12$ to $16$ as the matrix $\nabla_v \Phi(\tilde{\theta}, \tilde{v})$ in (ref) is not a contraction. For SA, larger values ($q = 8$ to $10$ or more) may be needed, particularly when $\beta$ is large. In the linear-utility case, the choice of $q$ has no effect on asymptotic efficiency ((ref)), so the only consideration is whether the outer loop converges and how long it takes.

\paragraph{Convergence check.} If the EM-NPL($q$) algorithm fails to converge at the initial $q$, increase $q$ and re-run. The convergence results in (ref) guarantee that convergence is eventually achieved for sufficiently large $q$ whenever the standard EM-NPL algorithm converges.

\paragraph{Computational cost considerations.} The total computational cost is $\text{CT}(q) \approx h(q) \times \bar{K}_{q}$ where $h(q)$ is the cost of $q$ inner iterations and $\bar{K}_{q}$ is the number of outer iterations to converge. It is generally U-shaped in $q$: too small a $q$ requires many outer iterations (large $\bar{K}_{q}$), while too large a $q$ wastes computation on unnecessary inner precision (large $h(q)$).

Monte Carlo Simulations

In this section, we present two Monte Carlo simulations to evaluate the performance of the proposed EM-NPL$(q)$ estimator. The first simulation considers a single-firm entry and exit with unobserved heterogeneity based on aguirregabiria2023solving. The second applies the estimator to a dynamic game as described in aguirregabiria2007sequential.

Single-Agent Entry and Exit with Unobserved Heterogeneity

Design of Simulation

We consider a single-firm entry and exit problem with unobserved heterogeneity, extending the DGP of aguirregabiria2023solving to a three-type mixture model ($M = 3$). At each period $t$, a firm decides whether to exit ($a_{t}=0$) or enter ($a_{t}=1$) the market. An active firm ($a_t=1$) earns a flow payoff determined by variable profits ($VP_t$), fixed operating costs ($FC_t$), and entry cost ($EC_t$); an inactive firm receives a normalized payoff of zero:

where $VP_t = (\theta_0^{VP} + \theta_1^{VP} z_{1t} + \theta_2^{VP} z_{2t}) \exp(w_t)$, $FC_t = \theta_0^{FC} + \theta_1^{FC} z_{3t}$, and $EC_t = \theta_0^{EC} + \theta_1^{EC} z_{4t}$. The state vector is $x_t = (w_t, z_{1t}, z_{2t}, z_{3t}, z_{4t}, a_{t-1})$, where $z_{jt}$ follows an AR(1) process $z_{jt} = 0.6 z_{jt-1} + \eta_{jt}$ for $j=1,\dots,4$ with independent standard normal innovations. We consider two scenarios for the productivity process $w_{t}$: (i) finite dependence (FD): $w_t = 0.2 + 0.6 w_{t-1} + \eta_t$; (ii) non-finite dependence (NFD): $w_t = 0.2 + 0.3 a_{t-1} + 0.6 w_{t-1} + \eta_t$. The state space is discretized using Tauchen's method tauchen1986finite with 6 grid points per variable, giving $|\mathcal{X}| = 2 \times 6^5 = 15{,}552$.\footnote{For the NFD case, the grid for $w$ spans $[\min(\mu_{d=0}, \mu_{d=1}) - n_{\sigma}\sigma_{LR},\, \max(\mu_{d=0}, \mu_{d=1}) + n_{\sigma}\sigma_{LR}]$, where $\mu_d$ is the long-run mean under action $d$ and $\sigma_{LR} = \sigma_w/\sqrt{1-\gamma_2^2}$.}

Firms are drawn from one of three unobserved types with type-specific parameter vectors $\theta_1^* = (1.5, 1.5, -0.3, -0.3, -0.2, -0.3, -1)$, $\theta_2^* = (0.2, 0.2, -0.2, -3.5, -2.0, -0.5, -3)$, and $\theta_3^* = (0.8, 0.8, -1.0, -1.5, -0.8, -3.0, -1)$, with mixing weights $\pi^* = (0.5, 0.3, 0.2)$. Type 1 has high variable profits and low costs, Type 2 has low profits and high entry costs, and Type 3 is intermediate. We simulate $N = 5{,}000$ firms, each observed for $T = 20$ periods, with discount factor $\beta = \{0.8, 0.95, 0.9999\}$. We perform 500 Monte Carlo replications.

We compare six fixed-point--algorithm combinations: (i) PV_GMRES solves the policy valuation linear system using GMRES; (ii) PV_SA uses successive approximation on the same system; (iii) BM_SA applies successive approximation to the Bellman equation; (iv) BM_AD uses Anderson acceleration on the Bellman equation; (v) BM_NT applies Newton's method to the Bellman equation; and (vi) EE_SA uses successive approximation on the Euler equation (available in the FD case only). For each method, we consider $q \in \{4, 6, 8, 10\}$ and report full-convergence baselines for PV_GMRES, BM_NT, and EE_SA. Convergence is declared when $\max\{\|\hat{P}^{(t)} - \hat{P}^{(t-1)}\|_\infty,\, \|\hat{\theta}^{(t)} - \hat{\theta}^{(t-1)}\|_\infty,\, \|\hat{V}^{(t)} - \hat{V}^{(t-1)}\|_\infty / (1 + \|\hat{V}^{(t-1)}\|_\infty)\} < 10^{-3}$. The inner solver tolerance is $10^{-8}$.

Results

Tables (ref)--(ref) report results across three discount factors ($\beta = 0.80, 0.95, 0.9999$), each with finite dependence (FD) and non-finite dependence (NFD) specifications. The results confirm several key patterns. First, the MSE of all policy-valuation-based estimators (PV_GMRES and PV_SA) is invariant to $q$: across all discount factors and both FD and NFD cases, MSE remains unchanged regardless of whether $q = 4$ or full convergence is used, confirming (ref). The Euler equation estimator (EE_SA) yields slightly higher MSE (e.g., 0.026 vs.\ 0.024 at $\beta = 0.95$ in FD). Second, the EM-NPL($q$) algorithm converges in virtually all 500 replications for every method and $q$ value tested, demonstrating the robustness of the algorithm in the mixture setting. Third, PV_GMRES achieves substantial computational savings over all alternatives. At $\beta = 0.95$, PV_GMRES($q = 4$) takes 2.88s in the FD case, compared to 9.18s for PV_SA, 8.11s for BM_SA, and 15.16s for BM_NT at the same $q$, roughly 3$\times$ faster than SA-based methods and 5$\times$ faster than Newton. Truncating the inner iterations also helps: PV_GMRES($q = 4$) is approximately 26% faster than the full-convergence GMRES baseline (2.88s vs.\ 3.88s in FD) and 17% faster in NFD (2.63s vs.\ 3.15s). Fourth, the speed advantage of GMRES over Newton grows with the discount factor: at $\beta = 0.9999$, the full-convergence GMRES takes 4.36s (FD) versus 38.45s for Newton.

Dynamic Games

This section considers a dynamic game as in aguirregabiria2007sequential. A similar data generating process is also used by aguirregabiria2021imposing.

Design of Simulation

There are $J = 7$ firms competing across independent markets. In each period $t$, firm $j$ simultaneously chooses whether to be active ($a_{jt} = 1$) or inactive ($a_{jt} = 0$). The flow payoff for firm $j$ is given by:

where $s_{t} \in \{1,2,3,4,5\}$ is the market size\footnote{It follows a Markov chain with transition matrix$

$}, and $a_{j,t-1}$ is firm $j$'s action in the previous period. The state vector is $x_{t} = (s_{t}, a_{1,t-1}, \ldots, a_{J,t-1})$, giving a state space of cardinality $|\mathcal{X}| = 5 \times 2^{7} = 640$. $\theta_{RS}$ governs the effect of market size on revenue, $\theta_{RC}$ captures the strategic interaction, $\theta_{EC}$ is the entry cost parameter, and $\theta_{FC,j}$ are firm-specific fixed costs.

The true parameter values are set to $\theta_{\text{FC}} = (1.9, 1.8, 1.7, 1.6, 1.5, 1.4, 1.3)$, $\theta_{\text{EC}} = 1$, and the discount factor $\beta = 0.95$. We consider two values of $\theta_{\text{RC}} \in \{2.4, 4\}$, referred to as the mildly unstable and very unstable cases following the terminology of aguirregabiria2021imposing.

Results

We compare five estimator--algorithm combinations. All methods initialize CCPs via the frequency estimator. (i) PV_GMRES solves the policy valuation linear system $(I - \beta F_P)V = U_P$ for each firm $j$ using GMRES, where $F_P$ is the $640 \times 640$ transition matrix induced by all firms' CCPs and the exogenous Markov chain. (ii) PV_SA replaces GMRES with successive approximation on the same linear system. (iii) BM_SA replaces the policy valuation step with the Bellman equation fixing other firms' CCPs, and iterates SA on the Bellman operator for $q$ steps. (iv) EPL_GMRES uses the EPL mapping ((ref)) and solves the resulting linear system using GMRES. Note that the EPL solves a larger linear system of dimension $J|\mathcal{A}||\mathcal{X}| = 8{,}960$. (v) PV_EPL_GMRES first iterates NPL($\infty$) until convergence, then switches to EPL_GMRES($q$). For each method, we consider $q \in \{4, 8, 12, 16\}$ inner iterations. Full-convergence baselines are reported for PV_GMRES, PV_Exact (direct matrix solver), EPL_Exact, EPL_GMRES, and PV_EPL_GMRES.

Following the spectral algorithm of aguirregabiria2021imposing, we update the CCPs using a Barzilai--Borwein update: $P^{(k+1)} = P^{(k)} - \alpha_k \bigl(P^{(k)} - \Lambda(\hat\theta^{(k)}, \hat Y^{(k)}, P^{(k)})\bigr)$ where the step size $\alpha_k$ is computed as $\alpha_k = \frac{\|\Delta P^{(k)}\|^2}{\langle \Delta P^{(k)},\, \Delta \Phi^{(k)} \rangle}, \ \Delta P^{(k)} = P^{(k)} - P^{(k-1)}, \ \Delta \Phi^{(k)} = \Phi^{(k)} - \Phi^{(k-1)}$ with $\Phi^{(k)} = P^{(k)} - \Lambda(\hat\theta^{(k)}, \hat Y^{(k)}, P^{(k)})$ denoting the fixed-point residual. For $k=0$, we set $\alpha_0 = \min(1, 1/\|\Phi^{(0)}\|)$.

Convergence is declared when $\max\{\|P^{(k)} - P^{(k-1)}\|_\infty,\, \|\hat{\theta}^{(k)} - \hat{\theta}^{(k-1)}\|_\infty,\, \|Y^{(k)} - Y^{(k-1)}\|_\infty / (1 + \|Y^{(k-1)}\|_\infty)\} < 10^{-3}$. The inner solver tolerance is $10^{-8}$. All results are averaged over 500 Monte Carlo simulations with $N = 1{,}600$ independent markets.

(ref) reports the results. Three patterns emerge. First, PV_GMRES is the fastest method by a wide margin: at $\theta_{\text{RC}} = 2.4$, PV_GMRES($q = 4$) takes 2.49s, compared to 3.61s for PV_SA, 5.17s for BM_SA, and 32.71s for EPL_GMRES. The gap widens under stronger strategic interaction: at $\theta_{\text{RC}} = 4$, PV_GMRES($q = 4$) takes 5.09s versus 16.48s for BM_SA and 75.61s for EPL_GMRES($q = 4$). PV_GMRES achieves 100% convergence at every configuration. BM_SA, by contrast, sees convergence drop to 93.8% at $\theta_{\text{RC}} = 4$ with $q = 4$, along with inflated MSE (0.493 vs.\ 0.415).

Second, the EPL-based estimators achieve substantially lower MSE than the NPL-based estimators, reflecting its statistical efficiency. At $\theta_{\text{RC}} = 4$, EPL_GMRES attains MSE of 0.299 compared to 0.415 for PV_GMRES, a 28% reduction. PV_EPL_GMRES, which switches from NPL to EPL after the first iteration, achieves a comparable MSE of 0.300. At $\theta_{\text{RC}} = 2.4$, the efficiency gain is more modest: EPL_GMRES yields MSE of 0.583 versus 0.617 for PV_GMRES. However, the EPL estimators are an order of magnitude slower than PV_GMRES because the EPL fixed-point mapping $\nabla_v \Phi(\tilde{\theta}, \tilde{v})$ is not a contraction, requiring substantially more inner GMRES iterations than the policy valuation system used by PV_GMRES. For instance, EPL_GMRES needs $q \geq 12$ to stabilize computation time and convergence at $\theta_{\text{RC}} = 4$ (31.58s at $q = 12$ versus 75.61s at $q = 4$), whereas PV_GMRES performs well with $q = 4$.

Third, the choice between PV_GMRES and EPL_GMRES reflects a speed--efficiency tradeoff. When computational speed is the priority (for example, when the state space is large), PV_GMRES with a small $q$ (e.g., $q = 4$) is the recommended default. When statistical efficiency matters, EPL_GMRES with a moderate $q$ (e.g., $q = 12$) delivers meaningful efficiency gains at the cost of longer computation time.

\paragraph{Summary.} Across both Monte Carlo designs, PV_GMRES with a small truncation level ($q = 4$) is the fastest method, achieving full convergence and MSE comparable to exact solvers. Truncation invariance holds empirically for all policy-valuation-based estimators, consistent with (ref). EPL-based estimators offer lower MSE at substantially higher computational cost, and EM-NPL($q$) reduces that cost by about 40%. These findings validate the practical recommendations in (ref): start with PV_GMRES at small $q$ for speed, and switch to EPL_GMRES at moderate $q$ when efficiency is the priority.

Empirical Application

We illustrate the EM-NPL($q$) estimator with an application to household demand for cola in London. We show that ignoring unobserved heterogeneity substantially biases estimates of price elasticities and counterfactual policy effects.

Data

We construct household-level purchase histories from the Worldpanel by Numerator Take Home data for London, covering the period 2015--2019. The panel records every grocery purchase made by a representative sample of households, including price, quantity, date, and product characteristics. The Cola category combines bottled and canned colas. We classify each product by its holding company (Coca-Cola vs.\ Pepsi vs.\ Other) crossed with its diet type (Regular, Diet, or Zero/Max), yielding seven choice alternatives: Coca-Cola Regular, Coca-Cola Diet, Coca-Cola Zero, Pepsi Regular, Pepsi Diet, Pepsi Max, and Other. The six named groups jointly account for 87.1% of category trips.

We aggregate multiple purchases by the same household on the same date to the trip level: when a household makes multiple cola purchases on a single day, we retain only the brand--type combination with the highest expenditure on that trip. We restrict the sample to London households and to those with more than 10 and fewer than 500 category trips over the five-year window.

After these filters, the sample contains 75{,}559 trips. The overall repeat rate is 70.9%, and the median inter-purchase spell is 9 days, reflecting the high purchase frequency of soft drinks. Table (ref) reports the choice shares and repeat-purchase rates.

Model

We assume that each household expects to make a purchase per median inter-purchase spell of 9 days. At each purchase occasion $t$, household $i$ of type $m$ chooses action $a_{it}\in\mathcal{A}=\{0,\ldots,6\}$ to maximize the expected discounted utility, where the seven alternatives correspond to the Cola choice groups defined in Table (ref): Coca-Cola Regular, Coca-Cola Diet, Coca-Cola Zero, Pepsi Regular, Pepsi Diet, Pepsi Max, and Other. The per-period payoff for type $m$ is

where $p$ denotes the vector of current prices, $a_{-1}$ is the previous purchase, $p_a$ is the unit price of product $a$, $\gamma_a^m$ captures type-specific product fixed effects, $\alpha^m$ is the type-specific price coefficient, and $\eta$ measures state dependence, common to all types. The error $\varepsilon(a)$ is i.i.d.\ Type-I extreme value. The type-specific parameter vector is $\theta^m = (\gamma_1^m,\ldots,\gamma_6^m, \alpha^m, \eta) \in \mathbb{R}^8$. We set the discount factor to $\beta=0.99$.

Prices follow a first-order vector autoregression. We discretize the estimated VAR(1) with 5 grid points per product, yielding $5^7 = 78{,}125$ price states. Combined with seven previous-action states, the augmented state space has $7 \times 78{,}125 = 546{,}875$ elements.

We estimate a finite mixture with $M$ latent types using the EM-NPL($q$) algorithm ((ref)) with $q=10$ GMRES iterations to solve the policy valuation problem. We randomly initialize the CCPs using a second-degree polynomial sieve logit. We choose 150 random starting values, and choose the one that maximizes the likelihood. For EM-NPL($q$), the criterion is $\max\{\| \boldsymbol{\theta}^{(k)} - \boldsymbol{\theta}^{(k-1)} \|, \|\mathbf{P}^{(k)} - \mathbf{P}^{(k-1)}\|, \|\boldsymbol{\pi}^{(k)} - \boldsymbol{\pi}^{(k-1)}\|\} \leq 10^{-4}$. We set the maximum number of EM-NPL($q$) iterations to 100. We estimate models for $M \in \{1, 2, \ldots, 10\}$.

As our utility function is linear in parameters and we solve the policy valuation problem, our estimator is numerically identical to the EM-NPL estimator that solves the inner fixed-point problem to convergence. Due to the Zero-Jacobian property aguirregabiria2002swapping, the asymptotic variance of the estimator simplifies to the inverse of the Hessian of the pseudo-likelihood with respect to the structural parameters. That is, the asymptotic variance does not depend on the nuisance parameters. To obtain standard errors for the counterfactual analysis, we draw 100 times from the asymptotic distribution of the estimated parameters and report the standard deviation across draws.

Estimation Results and Counterfactual Analysis

We estimate the model separately for three income groups defined by: low (under \pounds30{,}000), middle (\pounds30{,}000--\pounds49{,}999), and high (\pounds50{,}000 and above). Estimating separately by income group allows the distribution of latent types, and hence the pattern of brand preferences and price sensitivity, to differ across income segments.

(ref) reports the computational time of the EM-NPL($q$) algorithm for each $M$. We report the average time across the three income groups and 150 random initializations.\footnote{We also implemented the sieve-MLE initialization discussed in (ref). However, it takes around 10 hours to compute the sieve-MLE estimates for $M=10$, which is prohibitively long for our application. Also, the sieve-MLE initialization does not yield higher likelihood values than the random initialization. Therefore, we use the random initialization.} The figure shows that computation time scales approximately linearly in $M$. For $M=10$, the average runtime is about 20 minutes, which is feasible for empirical applications.

Short-Run and Long-Run Price Elasticities

We compute short-run and long-run price elasticities. The short-run elasticity measures the contemporaneous response to an unanticipated, temporary 1% price increase, holding continuation values fixed. The long-run elasticity permanently shifts product $j$'s price by 1% and re-solves the dynamic programming problem. Figures (ref)--(ref) plot both elasticities as a function of $M$ for each income group, and Table (ref) reports the full long-run cross-price elasticity matrices at the benchmark specification.

In general, own-price elasticities become more negative across all three income groups as additional latent types are introduced, indicating stronger substitution once the model absorbs persistent taste heterogeneity. The elasticity paths generally stabilize around $M=5$, so we use $M=5$ as the benchmark specification.

For the low-income group (Panel A of Table (ref)), the long-run own-price elasticity for Coca-Cola Regular shifts from about $-0.91$ under homogeneity ($M=1$) to $-1.78$ (s.e.\ 0.21) at the benchmark $M=5$, while the short-run elasticity shifts from about $-0.10$ to $-0.39$. The largest movement is for Pepsi Regular, whose long-run own-price elasticity goes from about $-0.70$ to $-3.15$ (s.e.\ 0.29) and whose short-run elasticity goes from about $-0.10$ to $-0.72$. Coca-Cola Diet, in contrast, moves only modestly to $-0.35$ (s.e.\ 0.15). Substitution is strongest between products of the same holding company: a 1% increase in the price of Coca-Cola Regular raises the long-run share of Coca-Cola Zero by 0.38 and of Coca-Cola Diet by 0.18. The Other category remains nearly price-inelastic throughout.

The middle-income group (Panel B) displays a qualitatively similar pattern. The long-run own-price elasticity for Coca-Cola Regular is $-1.78$ (s.e.\ 0.28) and for Coca-Cola Diet is $-0.99$ (s.e.\ 0.16). Cross-price effects are somewhat larger than for the low-income group: Coca-Cola Zero exhibits a long-run own-price elasticity of $-1.87$ (s.e.\ 0.40) and Pepsi Regular reaches $-3.01$ (s.e.\ 0.44), reflecting stronger substitution among diet, zero-sugar, and Pepsi varieties in this segment.

The high-income group (Panel C) exhibits the largest Coca-Cola elasticity magnitudes. The long-run own-price elasticity for Coca-Cola Regular is $-2.14$ (s.e.\ 0.20), and Pepsi Max reaches $-1.47$ (s.e.\ 0.24). Pepsi Regular displays the largest long-run own-price elasticity of $-3.07$ (s.e.\ 0.29), reflecting its small market share among high-income households. The short-run elasticities are especially sensitive to $M$ in all three groups: under homogeneity, consumers appear locked into their habitual brand, suppressing the immediate price response, but as heterogeneity is introduced this spurious persistence is reduced and the short-run response becomes substantially larger.

Figures (ref)--(ref) decompose the long-run own-price elasticity of Coca-Cola Regular across latent types for each income group, where types are ranked from the least to the most price-elastic and the mixture probabilities $\pi_m$ are reported below each bar. The figures reveal substantial heterogeneity in price sensitivity across types. In the low-income group, the most prevalent type (Type 4, $\pi_4=0.310$) has an elasticity of about $-2.1$, while the least elastic types (Types 1 and 2, jointly accounting for $\pi_1+\pi_2 = 0.377$ of households) have elasticities that are essentially zero with confidence intervals that include zero, and the most elastic type (Type 5) reaches about $-2.7$. In the middle-income group, Type 1 ($\pi_1=0.148$) has a slightly positive point estimate with a wide confidence interval that includes zero, and the more elastic types range from about $-1.7$ to $-2.8$. The high-income group displays a more uniform downward shift: the most prevalent type (Type 5, $\pi_5=0.303$) is also the most elastic at about $-2.3$, while the least elastic type (Type 1, $\pi_1=0.127$) is at about $-0.5$. Types with elasticities near zero likely represent brand-loyal consumers whose Coca-Cola Regular purchases are relatively insensitive to price changes.

Soda Tax Policy Analysis

We consider a sugar-sweetened beverage tax of \pounds0.25 per litre applied to the two regular (full-sugar) products, Coca-Cola Regular and Pepsi Regular. This policy experiment is broadly comparable to the UK Soft Drinks Industry Levy introduced in 2018. We use compensating variation (CV) to measure the consumers' welfare loss, defined as the amount of money needed to leave households indifferent between the baseline and taxed environments.

The left column of Figure (ref) plots the average CV as a function of $M$ for each income group. In all three groups, CV rises sharply between $M=2$ and $M=3$, jumping from less than \pounds2 to roughly \pounds3--\pounds4, and then levels off. Beyond $M=3$, successive values of $M$ yield CV estimates that are not statistically distinguishable. At the benchmark $M=5$, the average CV is about \pounds5.70 for the low-income group, \pounds4.80 for the middle-income group, and \pounds5.10 for the high-income group---between six and thirteen times the homogeneous ($M=1$) baseline of below \pounds1 in each group. dubois2020well reports the average CV of the UK Soft Drinks Industry Levy to be about \pounds4.94, which is close to our estimates for $M\geq3$ across all three income groups.

The right column of Figure (ref) decomposes CV by latent type at the benchmark $M=5$, with types ranked by long-run Coca-Cola Regular own-price elasticity as in (ref). In each income group, a single type accounts for almost the entire aggregate welfare loss, while all other types have CV close to zero. In the low-income group, Type 4 ($\pi_4=0.310$) has a CV of about \pounds17.4, while the CV of every other type is below \pounds1.3. In the middle-income group, Type 2 ($\pi_2=0.280$) has a CV of about \pounds14.8, versus no more than \pounds1.9 for the other types. In the high-income group, Type 5 ($\pi_5=0.303$) has a CV of about \pounds15.9, versus \pounds0.7 or less for the remaining types. In each case, the type bearing the welfare cost is the most prevalent type in that income group, and its mass-weighted contribution alone accounts for roughly 85--95% of the aggregate CV. The remaining types bear small welfare loss. In summary, Figure (ref) shows that ignoring unobserved heterogeneity leads to a substantial underestimation of the welfare cost of the tax, and that the welfare burden is highly concentrated among a single latent type in each income group, which is also the most prevalent type in that group.

Figures (ref)--(ref) decompose the welfare cost by household size, grouping households into four categories by the number of children (0, 1, 2, and $\geq 3$). In the low-income group, CV rises from about \pounds4.9 for childless households to \pounds6.1 for one-child households, peaks at \pounds9.1 for two-child households, and is \pounds7.3 for households with three or more children. The middle-income group is relatively flat across households with zero to two children (CV between \pounds4.6 and \pounds5.7) and falls to about \pounds3.4 for households with three or more. The high-income group is essentially flat across all categories, with CV between \pounds4.9 and \pounds5.7. In the middle- and high-income groups, the 95% confidence intervals overlap substantially across categories, so apparent differences by household size should be interpreted cautiously. In the low-income group, the peak at two children is statistically distinguishable from childless households and is the clearest pattern in the figure.

\FloatBarrier

Conclusion

This paper develops EM-NPL($q$) for dynamic discrete choice models with finite-mixture unobserved heterogeneity. The method embeds a truncated inner solver inside EM; for linear-in-parameters models, EM-NPL($q$) is numerically identical to EM-NPL for any $q \geq 1$, and we also establish consistency, asymptotic normality, and local convergence. Simulations show that GMRES delivers large speed gains, while EPL can improve accuracy in dynamic games at higher cost. In a cola-demand application estimated separately by income group, ignoring unobserved heterogeneity substantially understates price elasticities and the welfare cost of a soda tax, with the burden falling disproportionately on low-income households.

{0pt} \setcitestyle{authoryear,round}