Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
94,616 characters · 20 sections · 56 citation commands
Leveraging Uniformization and Sparsity for Estimation and Computation of Continuous-Time Dynamic Discrete Choice Games
\pdfauthor{Jason R. Blevins}
\abstract{ Continuous-time empirical dynamic discrete choice games offer notable computational advantages over discrete-time models. This paper addresses remaining computational and econometric challenges to further improve both model solution and estimation. We establish convergence rates for value iteration and policy evaluation with fixed beliefs, and develop Newton-Kantorovich methods that exploit analytical Jacobians and sparse matrix structure. We apply uniformization both to derive a new representation of the value function that draws direct analogies to discrete-time models and to enable stable computation of the matrix exponential and its parameter derivatives for estimation with discrete-time snapshot data, a common but challenging data scenario. These methods provide a complete chain of analytical derivatives from the value function for a given equilibrium through the log-likelihood function, eliminating the need for numerical differentiation and improving finite-sample estimation accuracy and computational efficiency. Monte Carlo experiments demonstrate substantial gains in both statistical performance and computational efficiency, enabling researchers to estimate richer models of strategic interaction. While we focus on games, our methods extend to single-agent dynamic discrete choice and continuous-time Markov jump processes.}
\acknowledgments{ Python replication code is available at \url{https://github.com/jrblevin/ctcomp}. I am grateful to Victor Aguirregabiria, Adam Dearing, Terence Johnson, and Larry Samuelson as well as participants of the 2024 Empirical Methods in Game Theory Workshop at Stony Brook University, the 2025 IIOC Conference, and the 2025 Midwest Econometrics Group conference for helpful feedback and discussions. }
Understanding the dynamics of market structure and competition in oligopolistic industries is a foundational goal of modern empirical industrial organization. Beginning with the seminal works of bresnahan91empirical and berry-1992, the field has transitioned from static models of market entry to dynamic models that capture the forward-looking behaviors of firms, influenced by methodological advancements from aguirregabiria-mira-2007, bajari-benkard-levin-2007, pakes-ostrovsky-berry-2007, and pesendorfer08asymptotic. However, this shift introduced significant computational challenges, where solving discrete-time dynamic equilibrium models becomes exponentially more complex as the number of players increases. In discrete-time games with simultaneous moves, computing equilibria requires forming expectations over all possible combinations of rivals' actions, and the complexity grows exponentially with the number of players. These computational barriers have limited empirical applications and counterfactuals to relatively simple models with few players and states, preventing researchers from analyzing realistic market structures. The prevalence of snapshot data---observations recorded at evenly spaced intervals without exact transition times---has driven the focus on discrete-time models with simultaneous moves in empirical work. However, the underlying economic reality often involves sequential decision-making---whether observed or not---with distinct strategic implications.
The introduction of continuous-time dynamic discrete choice games by doraszelski-judd-2012 and \citet*{abbe-2016} was aimed at addressing the computational limitations of discrete-time models. These models offer a more granular representation of decision-making processes which unfold in a stochastic, sequential manner and align with the nature of strategic interactions in many real-world settings. doraszelski-judd-2012 showed that continuous-time modeling simplifies the computation of players' expectations over future states, thereby reducing the computational costs relative to discrete-time models. Because only one player moves at any instant in continuous time, the order of operations required to compute equilibria scales linearly with the number of players rather than exponentially. \citet*{abbe-2016} developed an empirical framework and two-step CCP estimator based on the seminal contributions of hotz93conditional and hotz94simulation. More recently, blevins-kim-2024 extended the nested pseudo likelihood estimator of aguirregabiria-mira-2007 to the case of continuous-time games and ctgames considered identification and estimation of move arrival rates in the model, which were previously held fixed.
Despite these advancements, estimation of continuous-time models faces computational challenges that have limited their adoption. First, the convergence properties of value iteration had not been formally established, leaving researchers without theoretical guarantees about algorithm behavior. Second, reliance on numerical derivatives for optimization introduces approximation errors that propagate into parameter estimates, degrading finite-sample performance. Third, for researchers using snapshot data, computation of the matrix exponential $\exp(\Delta Q)$ is necessary for likelihood evaluation but is computationally intensive, particularly as model size grows.\footnote{ In contrast, continuous-time data with observed transition times allow likelihood computation directly from the intensity matrix $Q$ without matrix exponentials billingsley61statistical.} The sparsity of $Q$ is a key feature our methods exploit for computational efficiency in both equilibrium computation and matrix exponential calculation. While the value function methods we develop apply regardless of data structure, our matrix exponential algorithms specifically target the snapshot data case.
This paper provides theoretical foundations and computational methods to enhance the solution and estimation of continuous-time dynamic discrete choice models. In Section (ref), we review the model and assumptions. In Section (ref), we consider computation of value functions and equilibria. We establish the contractivity and modulus of contraction of the Bellman optimality operator when beliefs are held fixed, filling a gap in the theoretical foundations for this class of models. Building on this, we develop a polyalgorithm that combines the global convergence properties of value iteration with the quadratic convergence of Newton-Kantorovich methods. This approach, inspired by the nested fixed point (NFXP) algorithm of rust87optimal but adapted for continuous-time settings, switches algorithms based on convergence monitoring to achieve both reliability and efficiency in solving the model. iskhakov16comment reinforce the importance of Newton-Kantorovich iterations for the performance of NFXP in discrete-time models. We then develop a new representation of the value function based on the method of uniformization, introduced by jensen-1953 for continuous-time Markov chain analysis, which allows us to draw direct comparisons to the discrete-time dynamic discrete choice literature and to establish rates of convergence for policy evaluation.
In Section (ref), we turn to computation of the log-likelihood function with snapshot data. We show how to apply uniformization to develop an efficient method for simultaneous computation of the matrix exponential and its derivatives in a way that is optimized for the sparse structure of transition rate matrices in the model. This is important for computationally efficient and accurate estimation of the model's structural parameters. Using analytical derivatives eliminates a source of approximation error that affects not just computational efficiency but also the finite sample properties of estimators in practice.
While we develop these methods in the context of continuous-time dynamic games to address their specific computational challenges, the techniques have broader applicability. Our application of uniformization for computing matrix exponentials and their derivatives extends to any finite-state continuous-time Markov chain with sparse intensity matrices, making the approach valuable for single-agent dynamic discrete choice models, duration models with time-varying covariates, and other finite-state continuous-time econometric models.
To demonstrate the practical utility of the methods developed, in Section (ref) we conduct a series of Monte Carlo experiments using a dynamic entry/exit model with stochastic market demand. These experiments focus on evaluating the computational performance of our analytical gradient and sparse matrix exponential methods. We find that analytical gradients reduce computation time by 46% and function evaluations by 84% compared to numerical gradients, while sparse matrix methods achieve speedups of up to 40$\times$ when computing only the columns of the matrix exponential needed for likelihood evaluation, rather than computing the full dense matrix. By identifying these computational efficiencies, we facilitate the estimation and simulation of larger, more complex dynamic discrete choice models in continuous time, allowing for deeper insights into strategic interactions and policy interventions in oligopolistic markets.
We model strategic interactions in a continuous-time dynamic discrete choice framework with $N$ forward-looking agents, indexed by $i = 1, \dots, N$, who operate over an infinite time horizon $t \in [0, \infty)$. The model accommodates both the stochastic nature of decision times and the dynamic considerations of agents making decisions with persistent effects.
\paragraph{State Space}
At any instant $t$, the state of the system is represented by a state vector $x \in \mathcal{X}$ that includes all payoff-relevant information that is common knowledge among the agents. This vector typically includes both agent-specific attributes (e.g., market entry status, product quality) and exogenous market conditions (e.g., demand levels, available technology) that affect decisions and payoffs. The state space is finite, with $K = |\mathcal{X}|$ denoting the total number of states. We will generally work with the integer state space $\mathcal{K} = \{1, \dots, K\}$ and refer to states by their index $k$ or as an indexed vector $x_k$, as this allows us to vectorize many expressions.
\paragraph{Decisions and Endogenous State Changes}
Decision opportunities for agent $i$ in state $k$ arrive according to independent Poisson processes with finite rate parameters $\lambda_{ik}$, reflecting limited attention or other frictions on decision frequency. At each decision time, the agent observes $x_k$ and chooses an action $j$ from the choice set $\mathcal{J}_{ik}$. For expositional clarity, we focus on a common choice set: $\mathcal{J}_{ik} = \mathcal{J} = \{0, 1, \ldots, J-1\}$. The outcome of agent $i$'s decision $j$ in state $k$ is a deterministic state transition to $k' = l(i,j,k)$.
\paragraph{Exogenous State Changes}
An artificial player called “nature” ($i = 0$) governs exogenous state transitions, which follow a Markov jump process characterized by an intensity matrix $Q_0$. The $(k,k')$ element $q_{0kk'}$ represents the rate at which the system transitions from state $k$ to state $k'$ due to exogenous factors, with the aggregate exit rate from state $k$ given by $\sum_{k' \neq k} q_{0kk'}$.
\paragraph{Payoffs}
Agent $i$ receives flow payoffs $u_{ik}$ while in state $k$, discounted at rate $\rho_i$, with present value $\int_0^\tau e^{-\rho_i t} u_{ik}\,dt = u_{ik} \left(1 - e^{-\rho_i \tau} \right)/\rho_i$ over interval $[0,\tau)$. In contrast, instantaneous choice-specific payoffs $c_{ijk}$ are incurred when agent $i$ chooses action $j$ in state $k$. These payoffs are additively separable:
where $\psi_{ijk}$ is the deterministic mean payoff and $\varepsilon_{ijk}$ is an unobserved choice-specific shock. We normalize $\psi_{i0k} = 0$ for the continuation action $j=0$.
\paragraph{Assumptions}
Before describing how agents determine their actions in equilibrium, we formalize the main assumptions of the model, following \citet*{abbe-2016} and ctgames:
The previously mentioned Assumptions (ref)--(ref) lay the groundwork for the model's structure and dynamics. Assumption (ref) clarifies the role of the action $j = 0$, designating it as a continuation action that does not alter the current state. It also stipulates actions $j > 0$ are meaningfully different from one another, for identification purposes. Assumption (ref) describes our assumptions on the idiosyncratic error terms in the model and mirrors similar assumptions used in discrete time models aguirregabiria-mira-2010.
\paragraph{Strategies, Beliefs, and Value Functions}
A stationary Markov strategy for agent $i$, denoted $\delta_i$, maps each state $k \in \mathcal{K}$ and vector of choice-specific shocks $\varepsilon_{ik} \in \R^J$ to an action in $\mathcal{J}$. Any strategy $\delta_i$ determines choice probabilities
for all choices $j$ and states $k$.
To determine an optimal strategy, player $i$ must form beliefs about the choice probabilities of rivals. Let $\varsigma_{imjk}$ denote player $i$'s belief that rival $m$ chooses action $j$ in state $k$, with $\varsigma_{im} = \left(\varsigma_{imjk}\right)_{j\in\mathcal{J},k\in\mathcal{K}}$ denoting the beliefs about rival $m$, and
denoting player $i$'s complete beliefs about all rivals.
Given beliefs $\varsigma_i$, the value function for player $i$ is $V_i^{\varsigma_i} = (V_{i1}^{\varsigma_i}, \ldots, V_{iK}^{\varsigma_i})^\top$ where each element represents the expected present value of future rewards beginning in state $k$ and making optimal decisions thereafter. Following ctgames, the Bellman equation for $V_{ik}^{\varsigma_i}$ is:
This representation shows the balance between the immediate flow payoff in state $k$, the expected values following exogenous and endogenous state transitions, and the anticipated outcomes of player $i$'s decisions, all adjusted for the probability of occurrence and discounted appropriately. The expectation is over the joint distribution of the vector of shocks $\varepsilon_{ik}$. The denominator includes the total rate of events that can occur in state $k$, which we refer to as the maximum exit rate from state $k$, denoted as
\paragraph{Markov Perfect Equilibrium}
We now define the equilibrium concept used in this paper. A profile of stationary Markov strategies $\delta = (\delta_1, \ldots, \delta_N)$ determines a choice probability profile $\sigma = (\sigma_1, \ldots, \sigma_N)$, where $\sigma_i = \{\sigma_{ijk}\}_{j,k}$. A strategy $\delta_i$ is a best response to beliefs $\varsigma_i$ if it assigns action $j$ that maximizes expected discounted utility in state $k$:
We can characterize equilibria of the model as fixed points of a nonlinear, simultaneous system of equations for the value functions $T(V) = V$, where the system operator $T: \R^{N \times K} \rightarrow \R^{N \times K}$ applies each player's Bellman optimality operator given the beliefs $\varsigma_i(V_{-i})$ implied by the value functions of other players $V_{-i}$:
Existence of MPE was established by \citet*{abbe-2016} in models with homogeneous rates and by ctgames in the present model with heterogeneous rates.
The model's dynamics follow a finite-state continuous-time Markov chain (CTMC), or Markov jump process, denoted $X(t)$ for $t \in [0, \infty)$ with values in the state space $\mathcal{X}$. \footnote{We refer the reader to karlin75first, tijms-2003, or chung-1967 for details.} We characterize such a process using its $K \times K$ intensity matrix $Q = (q_{kk'})$ where, for $k' \neq k$,
represents the instantaneous rate at which transitions occur from state $k$ to state $k'$. The diagonal elements, $q_{kk}$, are set to $-\sum_{k' \neq k} q_{kk'}$, ensuring row sums are zero. The rate parameter for the exponential distribution of holding times before exiting state $k$ is $-q_{kk}$, equal to the aggregate off-diagonal transition rates out of state $k$. When the process exits state $k$, it transitions to state $k' \neq k$ with probability $q_{kk'} / (-q_{kk})$.
For estimation with snapshot data sampled at intervals of length $\Delta$, we require the transition probability matrix $P(\Delta) = \exp(\Delta Q)$, whose $(k,k')$ element gives the probability of transitioning from state $k$ to state $k'$ over time $\Delta$. The matrix exponential is defined by its power series expansion:
Despite its theoretical appeal, directly calculating the matrix exponential via (ref) is problematic. Matrix powers can be numerically unstable, and catastrophic cancellation may occur from alternating signs since $Q$ contains both positive and negative elements. Furthermore, the direct approach fails to exploit sparsity, as matrix powers of sparse $Q$ become dense as $j$ increases. While various numerical methods have been developed to address these issues, including Pad\'{e} approximations and scaling-and-squaring methods moler78nineteen,moler-van-loan-2003, we employ uniformization jensen-1953, a technique tailored specifically for generator matrices of continuous-time Markov chains. This approach naturally exploits sparsity when computing $\exp(\Delta Q) v$ for vectors $v$---a central feature of the models we consider.\footnote{Other contributions to uniformization include grassmann-1977-ejor, grassmann-1977-compor, reibman-trivedi-1988, and sherlock-2022. See vandijk-2018-uniformization for a recent overview.}
Uniformization converts a continuous-time Markov chain with varying exit rates into a discrete-time Markov chain subordinate to a Poisson process with rate parameter $\eta$ satisfying
This is achieved by introducing self-transitions that allow the process to remain in its current state. The discrete-time transition probability matrix is constructed as:
The probability of transitioning from state $k$ to state $k'$ in time $\Delta$ is then given by
where $[\Sigma^j]_{kk'}$ is the probability of transitioning from $k$ to $k'$ in exactly $j$ jumps, which is weighted by the Poisson probability of $j$ arrivals within interval $\Delta$. While any $\eta$ satisfying (ref) is valid, in practice using the smallest valid rate $\eta = \max_{k} \abs{q_{kk}}$ minimizes the number of terms required in (ref) to achieve a given accuracy (determined by the Poisson parameter $\eta \Delta$).
We now illustrate the construction and uniformization of equilibrium intensity matrices using the market entry and exit model introduced earlier. Recall that with $N = 2$ firms and 2 demand states, we have $K = 8$ total model states. The aggregate intensity matrix $Q$ decomposes as $Q = Q_0 + Q_1 + Q_2$, where $Q_0$ captures exogenous demand transitions and $Q_i$ captures firm $i$'s endogenous decisions.
Demand switches between H and L at rate $\gamma$, so $Q_0$ has the block form:
where $I_4$ is the $4 \times 4$ identity matrix. This structure reflects the fact that exogenous demand changes do not change firm activity configurations.
For simplicity, assume that firms make decisions at the same rate across states so that $\lambda_{ik} = \lambda$ for all $i$ and $k$. Then firm 1's entry and exit decisions create transitions at rate $\lambda \sigma_{11k}$:
Similarly, $Q_2$ captures firm 2's decisions with a different sparsity pattern.
The $8 \times 8$ aggregate intensity matrix (with diagonal elements omitted for brevity) is
The aggregate $Q$ matrix exhibits significant sparsity because direct transitions are only possible between states that differ in exactly one component (either demand or one firm's entry status). In this example, each row has at most 4 non-zero elements, yielding 50% sparsity. More generally, for $N$ firms and $D$ demand states, each state connects to at most $N + D$ other states, while the total state space has $2^{N} \times D$ elements. As we show later in Table (ref), with $5$ players and $3$ demand states ($K = 96$ total states), the $Q$ matrix is over 90% sparse, and the sparsity ratio approaches 1 as $N$ grows, making sparse matrix algorithms essential for computational tractability in realistic applications.
\definecolor{firmOneColor}{RGB}{0, 0, 255} \definecolor{firmTwoColor}{RGB}{34, 139, 34} \definecolor{demandColor}{RGB}{180, 20, 60} \definecolor{stateColor}{RGB}{255, 255, 255} \definecolor{stateBorderColor}{RGB}{0, 0, 0} \definecolor{groupColor}{RGB}{248, 248, 248}
To apply the uniformization technique, we select a uniform rate $\eta$ satisfying (ref) for any equilibrium $\sigma$: $\eta = 2\lambda + \gamma$. Given the $Q$ matrix in (ref), the uniform transition probability matrix $\Sigma = I + Q/\eta$ inherits the same sparsity pattern and has the following structure. For firm $i$ switching from state $k$ to $k'$, where states $k$ and $k'$ differ only in firm $i$'s activity status, the transition probabilities are $\lambda \sigma_{ijk}/\eta$. For demand transitions, where states $k$ and $k'$ have the same firm configuration but demand levels that differ by one unit, the transition probabilities are $\gamma / \eta$. Finally, the self-transition probabilities ensure that each row of $\Sigma$ sums to one are $\Sigma_{kk} = 1 - (\gamma + \lambda\sigma_{11k} + \lambda\sigma_{21k})/\eta$.
Figure (ref) illustrates this transformation, with panel \subref{fig:entry:unif:standard} showing the original continuous-time chain and panel \subref{fig:entry:unif:uniform} showing the uniformized discrete-time chain. States are grouped by demand level, with horizontal transitions representing firm 1's decisions and vertical transitions representing firm 2's decisions. Transitions between demand states for the same firm configuration each occur with the same rate and are represented by single arrows. The uniformization preserves the original dynamics while enabling more efficient computation.
This section develops computational methods for solving continuous-time dynamic discrete choice games. We first establish contraction properties of the Bellman operator and introduce a polyalgorithm that combines value iteration with Newton-Kantorovich iterations for efficient equilibrium computation. We then present a novel uniform policy-evaluation representation of the value function that connects continuous-time models to discrete-time frameworks.
We establish that the Bellman optimality operator $T_{i}^{\varsigma_i}$ defined in (ref), holding beliefs $\varsigma_i$ fixed, is a contraction with respect to the supremum norm. While analogous results are well-known for discrete-time dynamic programming, formal convergence guarantees have not been established for the continuous-time models considered here.
The proof of this theorem, and all remaining results, is given in Appendix (ref).
While Theorem (ref) established contractivity of each player's Bellman operator $T_{i}^{\varsigma_i}$ for fixed beliefs, the equilibrium system operator $T(V)$ may not be globally contractive because beliefs $\varsigma_i(V_{-i})$ are updated simultaneously with value functions. This makes equilibrium computation challenging, since multiple equilibria may exist.\footnote{ Equilibrium multiplicity is well-recognized in static and dynamic games. See depaula-2013, borkovsky-2015, aguirregabiria-collardwexler-ryan-2021, and otsu-pesendorfer-2023.} Therefore, we focus on methods for efficiently computing an equilibrium, and we consider three approaches: value iteration for finding equilibria that are stable under best-response dynamics, Newton-Kantorovich iterations for finding regular equilibria with good initialization, and a polyalgorithm that combines both methods. The polyalgorithm, inspired by rust87optimal's approach for single-agent problems, leverages the global properties of value iteration with the rapid local convergence of Newton-Kantorovich iterations. We discuss each method in turn.
Starting from any initial value function $V^{(0)}$, value iteration updates all players' value functions simultaneously: $V^{(n+1)} = T(V^{(n)})$. Without global contractivity, multiple equilibria may arise and value iteration alone can be slow to converge, particularly when the contraction modulus is close to one. Nevertheless, value iteration remains useful for finding equilibria that are stable under the best-response dynamics implicit in $T(V)$. An equilibrium $V^*$ is stable in this sense if the spectral radius of $\frac{\partial T}{\partial V}(V^*)$ is less than unity, ensuring that value iteration starting near $V^*$ will converge to $V^*$. Unstable equilibria---where the spectral radius exceeds unity---will have small or empty basins of attraction and will not be found unless initialization is extremely close to them.
Following nfxp, we monitor the convergence rate
to assess progress. In practice, we take the target modulus to be $\beta = \max_i \beta_i$ where $\beta_i$ is defined in (ref). When $r^{(n)}$ approaches $\beta$, the iteration is converging linearly at the expected rate, suggesting that switching to Newton-Kantorovich iterations, discussed in the next section, may be beneficial.
While value iteration provides convergence to stable equilibria, it can be slow when the contraction modulus is close to one. On the other hand, Newton-Kantorovich iterations exhibit quadratic local convergence, but rely on good initialization. The method solves the fixed-point equation $V = T(V)$ by reformulating the problem as finding a zero of $F(V) \equiv V - T(V)$ and applying Newton's method:
doraszelski-escobar-2010 show that regular Markov perfect equilibria---where the Jacobian $I - \frac{\partial T}{\partial V}(V^*)$ is nonsingular---are generic in discrete-time dynamic stochastic games with finite state and action spaces. Nonsingularity ensures local uniqueness via the Inverse Function Theorem and guarantees well-defined basins of attraction for Newton-Kantorovich iterations ortega-rheinboldt-1970.
To systematically search for unstable equilibria when needed, one can use Newton-Kantorovich iterations directly with varied initial guesses. The analytical Jacobian we derive can also serve to efficiently implement alternative equilibrium-finding approaches, such as homotopy methods besanko10learning,borkovsky10homotopy, which rely on repeated Jacobian evaluations along a continuation path.\footnote{ Other approaches for finding multiple equilibria include recursive lexicographical search iskhakov16recursive, which requires additional structure on the state space and action sets. Two-step estimation methods can mitigate issues of multiple equilibria, bypassing explicit equilibrium computation blevins-kim-2024. }
The key features that make Newton-Kantorovich tractable here are the reduced computational burden of applying the Bellman operator (linear complexity in $N$ rather than exponential) and the sparse structure of continuous-time models. Adapting this approach to multi-agent settings requires solving a potentially large linear system of equations and computing the analytical Jacobian of the equilibrium system---operations that would be computationally prohibitive for discrete-time games of comparable size but remain tractable in continuous time due to both the efficiency of computing expectations over future states and the sparsity of the Jacobian.
Implementing the Jacobian $\frac{\partial T}{\partial V}$ efficiently is important. We compute the Jacobian analytically by differentiating the equilibrium operator, with explicit formulas derived in Appendix (ref). Importantly, both the Jacobian and the system matrix $I - \frac{\partial T}{\partial V}$ inherit the sparsity of $Q$, enabling the use of efficient sparse linear solvers such as GMRES saad-schultz-1986. Table (ref) illustrates this sparsity for the entry/exit model from Section (ref), reporting state space size and sparsity of both $Q$ and $\frac{\partial T}{\partial V}$ across several model sizes ($N$ players and $D$ demand states). In models with 96 or more states, over 90% of the elements of both matrices are zeros.
The memory savings from sparse storage are substantial. For dense storage, the Jacobian requires $O(K^2 N^2)$ memory, which becomes prohibitive for larger models. For example, in the $10 \times 6$ model ($K = 6,144$), the Jacobian contains 3.8 billion elements. With 99.95% sparsity, we only store the 1.9 million non-zero elements---a reduction of three orders of magnitude. Similarly, the Newton-Kantorovich linear solve requires $O(K^3 N^3)$ operations with dense methods but only $O(K N)$ with sparse iterative solvers like GMRES, exploiting the fact that matrix-vector products preserve sparsity.
In addition to improving convergence, the analytical Jacobian enables efficient computation of value function derivatives $\frac{\partial V}{\partial \theta}$ via implicit differentiation:
Since the same system matrix appears in both the Newton-Kantorovich step and implicit differentiation, we factor it once and reuse this factorization to solve $\dim(\theta)$ linear systems---far more efficient than numerical differentiation, which would require solving $\dim(\theta)$ additional equilibrium problems. These derivatives propagate through the choice probabilities and intensity matrix to provide the log-likelihood gradient for maximum likelihood estimation, as we demonstrate below in Section (ref).
Following rust87optimal, we now combine value iteration and Newton-Kantorovich iterations into a polyalgorithm that leverages the strengths of each method. Value iteration establishes a basin of attraction around a stable equilibrium, then Newton-Kantorovich rapidly converges to the solution. The algorithm monitors the convergence rate $r^{(n)}$ during value iteration and switches to Newton-Kantorovich when $r^{(n)}$ approaches the target modulus $\beta$.
Algorithm (ref) presents the complete polyalgorithm, which monitors the convergence rate and switches from value iteration to Newton-Kantorovich when the target rate is approached. The algorithm requires specifying the convergence tolerance $\epsilon_V$ for the value function, rate tolerance $\epsilon_r$ for switching, and a minimum number of iterations MIN_MONITORING_ITER before monitoring begins. For example, in our Monte Carlo experiments we use $\epsilon_V = 10^{-13}$. We apply value iteration for at least 10 iterations before convergence monitoring, and we employ a relatively aggressive switching tolerance $\epsilon_r = 0.1$.
Figure (ref) demonstrates the polyalgorithm's convergence compared to pure value iteration. It achieves convergence in significantly fewer iterations (16 vs 2,071) by leveraging quadratic convergence of Newton-Kantorovich once the contraction rate approaches the theoretical modulus. Figure (ref) shows the convergence rate relative to $\beta$, with the switching condition met at iteration 14.
We now develop a novel uniform representation of the value function based on the uniformization method from Section (ref). While the representation builds on the CCP inversion result of hotz93conditional as applied to continuous-time settings by \citet*{abbe-2016}, our uniformization-based approach reveals a clear connection to discrete-time models and enables us to establish convergence rates for policy evaluation.
Theorem (ref) reveals a direct connection to discrete-time dynamic discrete choice models through uniformization. The representation defines a contraction with modulus $\overline\beta_i$, the effective discount factor, which depends on the uniformization rate $\eta$ relative to the discount rate $\rho_i$. Under Assumptions (ref) and (ref), $\overline\beta_i \in (0,1)$. The effective flow utility $U_i(\sigma)$ combines the actual flow utility with the rate-weighted expected instantaneous payoff at the agent's next decision time, while $\Sigma(\sigma)$ is a stochastic matrix serving as the effective transition probability matrix.
The lemma establishes that (non-uniform) value iteration converges at least as fast as uniform policy evaluation. The uniform representation trades some tightness in the contraction bound for the computational advantages of working with a stochastic matrix $\Sigma$ and the connection to discrete-time methods.
This section has developed computational methods for solving continuous-time dynamic discrete choice games, including a polyalgorithm for computing equilibria and a uniform policy evaluation representation that connects to discrete-time models. With these established, we now turn to the computational aspects of estimating such games using only discrete-time data.
While this paper focuses on computational methods for estimation, ctgames provides a comprehensive treatment of identification in this framework. The identification analysis proceeds in two steps: first recovering the intensity matrix $Q$ from discrete-time transition probabilities $P(\Delta)$, then recovering structural primitives (payoffs, value functions, and heterogeneous arrival rates $\lambda_{ik}$) from $Q$.
The first step addresses the aliasing problem: multiple intensity matrices $Q$ may generate the same discrete-time transitions $P(\Delta) = \exp(\Delta Q)$. However, the structural model provides natural restrictions that ensure identification. For instance, in entry/exit models, firms cannot force rivals to enter or exit, implying zeros in specific positions of the $Q$ matrix. Similarly, exogenous state variables (like demand conditions) can only be changed by nature, not by players' actions. These exclusion restrictions, arising naturally from the economic model, restrict the domain of possible $Q$ matrices to ensure unique identification from $P(\Delta)$. Alternative approaches include using sufficiently small sampling intervals $\Delta$ where $P(\Delta) \approx I + \Delta Q$, or observing transitions at multiple distinct intervals.
Given $Q$, the structural primitives can be identified with $2K$ linear restrictions per player---substantially fewer than the exponentially many restrictions required in discrete-time models pesendorfer08asymptotic. This reduction occurs because continuous-time models avoid simultaneous moves: only one player acts at any instant, simplifying the identification problem. Common restrictions include normalizations (e.g., $\psi_{i0k} = 0$ for continuation), constant parameters across states (e.g., entry costs independent of demand), or known values in terminal states. Importantly, heterogeneous arrival rates $\lambda_{ik}$ are not automatically identified from discrete-time data and require these additional restrictions.
The sampling frequency $\Delta$ affects both identification and estimation. Smaller $\Delta$ aids identification by making $P(\Delta)$ closer to the underlying $Q$, but the computational methods we discuss next are designed to handle the matrix exponential calculations for any fixed $\Delta$.
Estimation of the model with continuous-time data is relatively straightforward, so here we focus on the more computationally difficult task of estimation with discrete-time data.\footnote{With continuous-time data where exact transition times are observed, estimation is computationally simpler as the likelihood depends directly on the intensity matrix $Q$ rather than requiring matrix exponential calculations billingsley61statistical. Our focus on discrete-time snapshot data addresses the computationally more difficult case that is common in empirical practice.} This is the typical case in many empirical applications where data are collected annually, quarterly, or monthly, but the underlying process evolves continuously. Suppose the researcher observes snapshots of markets at regularly-spaced intervals $\Delta > 0$. Let $M$ denote the number of markets and $T_m$ the number of snapshots over time observed in each market $m = 1, \dots, M$. The sample can be represented as $\lbrace k_{m,0}, \dots, k_{m,T_m} \rbrace_{m = 1}^{M}$, where $k_{mn}$ denotes the observed state for snapshot $n$ in market $m$. Given the sample, we first pre-calculate transition counts $d_{kk'}$ for all pairs of states $(k, k')$:
The log-likelihood function can then be expressed simply in terms of these counts as
where $d_{k'}^\top = ( d_{1k'}, \dots, d_{Kk'} )$ contains transition counts into state $k'$ and $e_k$ is the $k$-th basis vector.
The form of (ref) enables several important computational efficiencies. First, we need to compute only those columns of $P(\Delta,\theta)$ corresponding to destination states $k'$ with positive observed transition counts. Individual columns can be obtained through matrix-vector products of the form $P(\Delta,\theta) e_{k'}$, substantially reducing the computational burden when the number of observed destination states is much smaller than the size of the state space---a common scenario in large models with moderate sample sizes. Second, for each required column $k'$, we only need the inner product $d_{k'}^\top \ln P(\Delta,\theta) e_{k'}$ using the non-zero elements of $d_{k'}$.
The gradient of $\ell(\theta)$ will also be important for accurate and efficient optimization:
Under standard regularity conditions and assuming a unique equilibrium, the maximum likelihood estimator obtained by maximizing (ref) has the usual asymptotic properties: consistency, asymptotic normality, and efficiency newey94large. Our analytical derivatives provide the exact score function $\nabla \ell(\theta)$ in (ref), eliminating numerical approximation error in both optimization and standard error computation.
Recall that the transition matrix $P(\Delta, \theta)$ appearing in (ref) and (ref) is computed via the matrix exponential as $P(\Delta,\theta) = \exp(\Delta Q(\theta))$. Therefore, to maximize the log-likelihood function we shall require an efficient method to compute $\exp(\Delta Q(\theta)) v$, the action of the matrix exponential on a probability vector $v$, along with the derivatives $\frac{\partial}{\partial \theta}\exp(\Delta Q(\theta))v$. The following section presents a numerically stable uniformization-based algorithm to compute these quantities efficiently.
Given a valid uniformization rate $\eta$ satisfying (ref), define $\Sigma \equiv I + Q / \eta$, so that $\Delta Q = \eta \Delta \Sigma - \eta \Delta I$. Then following (ref), the action of $\exp(\Delta Q)$ on a vector $v$ can be written as
Importantly, all elements of $\Sigma$ are non-negative, meaning that this calculation will not suffer from cancellation of alternating positive and negative terms goldberg-1991. As a result, computations will be much more numerically stable for the uniformization of a rate matrix $Q$ than for the exponential of a generic matrix $A$. Furthermore, we can compute (ref) efficiently using the following recurrence:
with $w_j$ involving only a single sparse matrix vector product sherlock-2022.
In practice, we truncate the series at $j = \bar{J}_\varepsilon < \infty$ such that the approximation error is below some tolerance $\varepsilon$. Since uniformization produces a Poisson process, we determine $\bar{J}_{\varepsilon}$ using Poisson tail probabilities fox-glenn-1988, reibman-trivedi-1988. When $\Sigma$ is a stochastic matrix and $v$ is a probability vector, each term $\Sigma^j v$ is also a probability vector with $\norm{\Sigma^j v}_1 = 1$. This property allows us to bound the truncation error directly: we choose $\bar{J}_\varepsilon$ as the smallest integer such that
which gives $\bar{J}_{\varepsilon} = \text{InvPoissonCDF}(1-\varepsilon; \eta \Delta)$.\footnote{sherlock-2022 discusses accurate computation of Poisson tail probabilities for small $\varepsilon$ and modifications to handle overflow in $\sum_{j=0}^J \frac{(\eta \Delta)^j}{j!}$.}
Uniformization of the matrix exponential also leads to an efficient, recursive algorithm for computing its derivatives rupp-2024. Returning to (ref) and differentiating, noting that $\frac{\partial (\Sigma^0)}{\partial \theta} = \frac{\partial I}{\partial \theta} = 0$, so we can start the sum at $j=1$:
Note that $\frac{\partial \Sigma^j}{\partial \theta}$ denotes the derivative of the matrix $\Sigma^j$, as opposed to the $j$-th power of the derivative of $\Sigma$. Computing this directly is infeasible, but we can compute it recursively from $\frac{\partial \Sigma}{\partial \theta}$ by noting that $\Sigma^j = \Sigma \Sigma^{j-1}$ and applying the product rule:
Substituting (ref) into (ref):
where we can compute $w_{j}$ and $\delta_{j}$ recursively as follows:
A pseudocode implementation is presented below as Algorithm (ref).
Finally, we return to computing the derivatives $\frac{\partial Q}{\partial \theta}$ which are required as inputs for the recursive algorithm above. These involve both direct parameter effects and indirect effects through choice probabilities:
where the subscript notation indicates the direct effect holding $\sigma$ fixed. Each element of $Q$ has the form:
with diagonal elements determined by the row-sum constraint. For off-diagonal elements ($k \neq k'$), the total derivatives combining both direct and indirect effects are
where the first two terms ($\frac{\partial q_{0kk'}}{\partial \theta}$ and $\frac{\partial \lambda_{ik}}{\partial \theta} \sigma_{ijk}$) capture direct parameter dependence of transition rates (e.g., $\gamma$ or $\lambda_{ik}$), and the final term ($\lambda_{ik} \frac{\partial \sigma_{ijk}}{\partial \theta}$) captures the indirect effect through equilibrium choice probability derivatives, derived below. Diagonal element derivatives follow from the row-sum constraint: $\frac{\partial q_{kk}}{\partial \theta} = -\sum_{k' \neq k} \frac{\partial q_{kk'}}{\partial \theta}$.
Recall that the equilibrium choice probabilities are determined by the Markov perfect equilibrium condition. For type I extreme value errors, these take the familiar logit form:
The logit structure provides closed-form derivatives. Differentiating with respect to continuation values yields:
That is, the derivative with respect to action $j$'s own continuation value is $\sigma_{ijk}(1-\sigma_{ijk})$, while the derivative with respect to a different action's continuation value is $-\sigma_{ijk}\sigma_{ij'k}$. The derivatives $\frac{\partial \sigma}{\partial \psi}$ follow the same functional form, while $\frac{\partial \psi}{\partial \theta}$ depends on the specific parameter (e.g., entry costs).
This completes the computational pathway for evaluating both the log-likelihood function and its gradient with respect to the structural parameters, allowing exact gradient computation for maximum likelihood estimation while exploiting sparsity and uniformization for efficiency and numerical stability.
We evaluate the finite-sample properties of our estimation methods through Monte Carlo simulations, focusing on how analytical derivatives improve both statistical accuracy and computational efficiency. As is common in empirical practice, we work with discrete-time snapshot data where the econometrician observes states only at fixed intervals $\Delta$, making matrix exponential computation necessary for likelihood evaluation.
The model for our Monte Carlo experiments is based on the example model of entry and exit with $N$ firms operating in a single product market. Firms operate under stochastically varying market demand conditions with $D$ demand levels $x_{k0} \in \{ 0, \dots, D - 1 \}$. The payoff relevant states are the number of active firms, denoted $n_{k}$, and the current state of market demand, $x_{k0}$.\footnote{In principle we could use this to reduce the state space of the model along the lines of \citet*{abbe-2016}, but we kept the full state space to keep the source code easy to understand and better illustrate how the computational complexity varies with the number of firms.} For our experiments, we specify $N = 7$ players and $D = 5$ demand states, yielding $K = 2^{N} \times D = 640$ total states.
Three structural parameters determine the flow payoffs and instantaneous payoffs: $\theta_{\text{EC}}$ is the entry cost (incurred when inactive firms enter), $\theta_{\text{RN}}$ represents competitive effect, and $\theta_{\text{D}}$ represents the profitability of the market level of demand, $x_{k0}$. For firm $i$ in state $k$, the flow payoff $u_{ik}$ is
where $x_{ki} \in \{0,1\}$ is the activity indicator for firm $i$ when the market state is $k$. The exit scrap value is normalized to zero.
We generate discrete-time datasets with sampling interval $\Delta=1.0$ and two sample sizes: $T \in\{ 1000, 4000 \}$ observations per replication. The true parameter values are given in Table (ref). With $\lambda = 1.0$, each player makes on average one decision per observation interval.\footnote{Both numerical and analytical gradient approaches use the same starting values, which are chosen to be far from the true parameters: $(\theta_{\text{EC}}, \theta_{\text{RN}}, \theta_{\text{D}}, \lambda, \gamma) = (-1.0, -0.1, 1.0, 0.2, 1.0)$ compared to true values of $(-2.0, -0.5, 2.0, 1.0, 0.3)$. The value function convergence tolerance is set to $10^{-13}$ for all experiments.}
As noted in Section (ref), the equilibrium system operator $T(V)$ may not be globally contractive, allowing for multiple equilibria. However, for this specification there appears to be a unique equilibrium.\footnote{ We searched for multiple equilibria by solving the equilibrium system $T(V)$ from 10,000 random initial value vectors with components uniformly distributed over $[\min u_{ik} / \rho, \max u_{ik} / \rho]$ using Newton-Kantorovich iterations. We solved the system to a tolerance of $10^{-13}$, then clustered the choice probability vectors using a tolerance of $10^{-3}$. Each of the solutions converged to the same equilibrium. } This allows us to focus on evaluating our matrix exponential and analytical derivative methods for a given equilibrium.
We compare three estimation approaches over 100 Monte Carlo replications:
All approaches maximize the log-likelihood function $\ell(\theta)$ in (ref) using L-BFGS-B and sparse matrix routines from the Python SciPy library. For the infeasible start results, we initialize the optimization using the best parameter values from three choices: (1) the true parameter values, (2) the estimates from the numerical gradient approach, and (3) the estimates from the analytical gradient approach.
Tables (ref) and (ref) present results from 100 Monte Carlo replications for sample sizes of 1,000 and 4,000 observations, respectively. Our results reveal that beyond computational gains, analytical derivatives substantially improve small-sample estimation performance.
With the smaller sample size (Table (ref)), analytical gradients show dramatically better performance in terms of variance. The standard deviation of the $\theta_{\text{EC}}$ estimates is 0.239 with numerical gradients versus 0.204 with analytical gradients. More strikingly, the estimates of $\lambda$ show standard deviations of 0.437 versus 0.057, respectively---a nearly eight-fold improvement. This dramatic difference suggests that numerical gradient approximation errors contaminate the optimization process, preventing the optimizer from reaching the true maximum and introducing additional variance into the estimates. The analytical gradient results are very similar to the infeasible baseline results, suggesting that the analytical gradient method achieves near-optimal convergence.
The outliers from using numerical gradients can be seen clearly in the box plots shown in Figure (ref). With analytical gradients, we have fewer extreme outliers, providing substantial improvements in finite-sample performance. Both approaches have little bias despite estimation with only discrete-time data where the order of moves is unobserved, but the variance reduction from analytical gradients is substantial across all parameters.
With the larger sample size (Table (ref)), the differences in standard errors between numerical and analytical gradients diminish as expected, though analytical gradients still provide improved precision. The convergence of both methods as the sample size increases confirms the asymptotic properties discussed in Section (ref).
Beyond the statistical improvements, analytical gradients provide significant computational efficiency gains. Figure (ref) provides a detailed breakdown of computational performance across the Monte Carlo replications. Panel (a) shows that using analytical gradients yields accurate optimal log-likelihood values up to $10^{-11}$ precision, while numerical gradients exhibit larger dispersion. Panel (b) shows that both methods require approximately the same number of iterations, but panels (c) and (d) reveal large differences in the number of functional evaluations and computation times.
From Table (ref), we see that using analytical gradients reduces the average estimation time from 605.7 to 329.4 seconds---a 46% reduction---and reduces the average number of functional evaluations from 370.9 to 59.8---an 84% reduction. This is primarily due to the elimination of expensive finite-difference gradient approximations at each iteration, which involves solving the equilibrium of the model for each functional evaluation. The analytical derivatives are computed efficiently using Algorithm (ref) with implicit differentiation for the equilibrium value functions.
Finally, to focus on the performance of the sparse matrix exponential algorithm from Section (ref), in Table (ref) we compare the standard dense matrix exponential function expm from SciPy (based on Pad\'{e} approximation), with the uniformization-based sparse version described in Algorithm (ref). First we compare computation of the full matrix exponential $\exp(Q)$ for each model. For smaller models, there is no advantage to using sparse methods to compute the full matrix exponential, but for the largest model we do begin to see a speedup. For the models considered, the most significant benefit of sparse methods occurs when we only need to compute a relatively small number of columns, such as when evaluating the log-likelihood function over a finite sample. As we showed in Section (ref), we only need to compute columns corresponding to observed destination states, which is typically much smaller than $K$ for moderate sample sizes. Under the “Columns of $\exp(Q)$” heading in Table (ref), where we compute $\min\{200, K\}$ columns of $\exp(Q)$ (e.g., for a sample of 200 observations), we see speedups of up to $40.48\times$. Comparing the accuracy of the sparse method using the dense method as the baseline, we see that even for very large models the relative errors are all on the order of $10^{-8}$ or smaller, showing that the sparse uniformization-based method is both efficient and accurate for maximum likelihood estimation.
The choice between sparse and dense methods depends on both the state space dimension $K$ and the number of matrix exponential columns required. Dense methods compute the full $K \times K$ matrix exponential at cost $O(K^3)$, making them efficient for small problems or when many columns are needed. Sparse methods avoid this cubic scaling by computing individual columns at cost approximately $O(K)$ per column, becoming advantageous when computing only $c$ columns with $c \ll K$. In maximum likelihood estimation with sample size $T$, the number of required columns is bounded by the number of unique destination states observed, typically $O(\min\{T, K\})$. Thus, sparse methods are recommended when $T \ll K^2$, which encompasses most empirically relevant settings with moderate to large state spaces.
This paper provides theoretical foundations and practical methods for estimation of continuous-time dynamic discrete choice models, with applications ranging from single-agent problems to multi-agent games. Our contributions address both econometric and computational challenges that have limited the adoption of continuous-time methods despite their theoretical advantages.
From an econometric perspective, we provide formal convergence guarantees for value iteration with fixed beliefs in continuous-time settings and develop a complete analytical derivative chain for maximum likelihood estimation. These analytical derivatives are not merely a computational convenience---they eliminate approximation errors from numerical gradients that affect finite-sample performance, as our Monte Carlo experiments demonstrate with up to eight-fold reductions in parameter standard errors. This has significant implications for empirical researchers working with real-world datasets.
From a computational perspective, our application of uniformization exploits the sparse structure inherent in many economic models to achieve substantial efficiency gains. Specifically, we have (1) established theoretical convergence properties for value iteration with fixed beliefs, (2) developed a polyalgorithm combining value iteration with Newton-Kantorovich methods for robust equilibrium computation, (3) introduced a uniformization-based representation that provides direct analogies with discrete-time models, and (4) demonstrated efficient estimation with discrete-time snapshot data using analytical derivatives of sparse matrix exponentials. Our experiments show that sparse methods achieve speedups of up to 40$\times$ when computing only the columns needed for likelihood evaluation, while analytical gradients reduce computation time by 46%.
By reducing computational barriers and improving finite sample statistical properties, these methods enable researchers to estimate richer, more realistic models that were previously infeasible. The techniques apply to other models involving finite-state continuous-time Markov chains, including single-agent dynamic discrete choice, duration models with time-varying covariates, and regime-switching models. This opens new possibilities for empirical work in industrial organization, labor economics, finance, and related fields.