{'fraction_non_alphanumeric': 0.05063044642422472, 'fraction_numerical': 0.013387858429482499, 'mean_word_length': 4.32087812459526, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Temporal-Difference estimation of dynamic discrete choice models

We study the use of Temporal-Difference learning for estimating the structural parameters in dynamic discrete choice models. Our algorithms are based on the conditional choice probability approach but use functional approximations to estimate various terms in the pseudo-likelihood function. We suggest two approaches: The first - linear semi-gradient - provides approximations to the recursive terms using basis functions. The second - Approximate Value Iteration - builds a sequence of approximations to the recursive terms by solving non-parametric estimation problems. Our approaches are fast and naturally allow for continuous and/or high-dimensional state spaces. Furthermore, they do not require specification of transition densities. In dynamic games, they avoid integrating over other players' actions, further heightening the computational advantage. Our proposals can be paired with popular existing methods such as pseudo-maximum-likelihood, and we propose locally robust corrections for the latter to achieve parametric rates of convergence. Monte Carlo simulations confirm the properties of our algorithms in practice.


Introduction
Recent years have seen a number of important developments in the field of Reinforcement Learning (RL) for computation of value functions. The goal of this paper is to study the use of a popular RL technique, Temporal-Difference (TD) learning, for estimation and inference in Dynamic Discrete Choice (DDC) models.
DDC models are frequently used to describe the inter-temporal choices of forward-looking individuals in a variety of contexts. In these models, agents maximize their expected future payoff through repeated choice amongst a set of discrete alternatives. Based on a revealed preference argument, structural estimation proceeds by using microdata on choices and outcomes to recover the underlying model parameters. 1 A key challenge in this literature is the complexity of estimation. Uncovering the structural parameters originally required an explicit solution to a dynamic programming problem in addition to the optimization of an estimation criterion. A key advance has been Hotz and Miller's (1993) Conditional Choice Probability (CCP) algorithm which avoids the repeated solution of the inter-temporal optimization problem by taking advantage of a mapping between value function differences and conditional choice probabilities.
Unfortunately, the standard CCP algorithm is computationally infeasible when the underlying states are continuous and/or the state space is high-dimensional. Such state spaces are common in applications. One frequently employed approach to tackle continuous state spaces is through state space discretization, e.g., Kalouptsidi (2014) and Almagro and Domınguez-Iino (2019) use aggregation and clustering methods to do this. However, it is not always clear how to perform such a discretization in practice, and moreover, it introduces bias into the parameter estimates. An alternative is to employ functional approximations for the value functions. For instance, Barwick and Pathak (2015) and Kalouptsidi (2018) use estimated transition densities and numerical/analytical integration to approximate the value functions using linear sieves and LASSO, respectively. However, the theoretical properties of these methods when using machine learning methods (such as LASSO) are as yet unknown, and they still require estimation of transition densities, which is not straightforward, along with numerical integration, which can be computationally expensive. 2 The aim of this paper is to develop tractable algorithms for CCP estimation when the state variables are continuous and/or the state space is large. Such algorithms should possess three properties: First, they should be fast to compute even under high-dimensional state spaces. Second, they should avoid state space discretization, and instead rely on functional approximation of value functions. Third, they should avoid estimation of transition densities which are difficult to parameterize and estimate under continuous states.
In this paper, we suggest two methods, based on the idea of TD learning, that satisfy all the above properties. The methods involve two different techniques for estimating various recursive terms (which are akin to value functions) that arise in CCP estimation.
The first approach, the linear semi-gradient method, provides functional approximations to the recursive terms using basis functions. This approach simply involves inverting a low-dimensional matrix, where the dimension is the number of basis functions being used. Thus the computational cost is trivial in most settings. Furthermore, it only requires the observed sequences of current and future state-action pairs as input and estimation of transition densities is not needed. The second approach, Approximate Value Iteration (AVI), builds a sequence of approximations to the value terms by solving a non-parametric estimation problem in each step. Almost any machine learning (ML) method devised for prediction can be used for approximation under this method, including (but not limited to) LASSO, random forests and neural networks. To the best of our knowledge, the AVI estimator is the first estimator for DDC models that can be applied with any ML method. Hence, it naturally allows for very high-dimensional state spaces.
Again, no estimation of transition densities is required. We derive the non-parametric rates of convergence for estimation of the value terms under both these methods. With the estimates of these functions at hand, estimation of the structural parameters can proceed with standard methods such as pseudo-maximum-likelihood estimation (PMLE, Aguirregabiria and Mira (2002)) or minimum distance estimation.
The focus of this paper is on the estimation of structural parameters. To this end, our procedures allow one to avoid modeling state transitions. Performing counterfactual analysis may still require estimating the transition density, but we argue that our techniques remain useful, even for this purpose, for two reasons: First, counterfactuals often involve new transition densities which are different from the ones that enter the estimation of the structural parameters, see Kalouptsidi (2018) for an example. Our methods therefore avoid estimation of the original transition densities that may not be needed for counterfactual analysis. Second, with continuous states, decoupling the estimation of structural parameters and transition densities may be beneficial for robustness and efficiency reasons. For instance, it is common to employ AR (e.g., Aguirregabiria and Mira, 2007;Kalouptsidi, 2014) or VAR (Barwick and Pathak, 2015) specifications for transition densities. However, even within these specifications there are a number of choices to be made (e.g., dimension of VARs, distribution of error terms etc) and the estimates of the structural parameters may not be robust to these choices. More importantly, even when non-parametric estimates of transition densities are available, simply plugging them into the second-stage PMLE criterion would seriously degrade the rate of convergence of structural parameters. One would need to adjust the PMLE to account for the non-parametric first stage, but it is not known what form this adjustment takes.
By contrast, our proposals employ non-parametric estimates of value functions, and as described below, we derive the necessary adjustments to account for this. To carry out counterfactual analysis, we therefore suggest plugging in our estimates of the structural parameters, which do not rely on non-parametric estimates of transition densities and are robust to mis-specification, together with non-parametric estimates of the transition densities.
The above discussion highlights that in continuous state spaces, estimation of structural parameters is inherently a problem of semi-parametric estimation. In fact, even under discrete states, Aguirregabiria and Mira (2002) show that estimation of transition densities affects the variance of the structural parameter estimates. If the state variables are continuous, existing two-step CCP methods, such as the PMLE will no longer converge at parametric rates. We therefore derive a locally robust estimator by adding a correction term to the PMLE criterion function that accounts for the non-parametric estimation of value function terms using either of our TD methods. We emphasize that this construction is novel and does not directly follow from existing results, e.g., in Chernozhukov et al. (2022). Its computation is particularly straightforward under the linear semi-gradient approach. The resulting estimator converges at parametric rates under continuous states and unrestricted transition densities.
Our TD estimators are thus consistent, computationally very cheap, and they converge at parametric rates. Importantly, they provide a feasible estimation method when the states are continuous and/or the state space is large. This is particularly important for the estimation of dynamic discrete games. Even with discrete states, existing methods for the estimation of dynamic games (Bajari et al., 2007;Aguirregabiria and Mira, 2007;Pesendorfer and Schmidt-Dengler, 2008) require integrating out other players' actions. With many players, or under continuous states, this can get quite cumbersome. By contrast, our procedure works directly with the joint empirical distribution of the states and their sample successors. Thus the 'integrating out' is done implicitly within the sample expectations.
Finally, we also incorporate permanent unobserved heterogeneity into our methods by combining the TD estimation with an Expectation-Maximization (EM) algorithm (Dempster et al., 1977).
A range of Monte Carlo studies confirm the workings of our algorithms. First, we present simulations based on the dynamic firm entry problem described in Aguirregabiria and Magesan (2018). The model has seven structural parameters of interest and five continuous state variables. Existing methods usually struggle at this dimensionality; certainly, discretization of the state space would not work too well. We show simulations using both the linear semi-gradient method and the AVI method where we estimate the value function terms using random forests, and derive results both with and without locally robust correction. Our findings suggest that while the linear semi-gradient has very little bias even without the locally robust correction, the AVI method has slightly higher small sample bias which is substantially reduced when using our locally robust estimator. Overall, our estimators perform very well in this setting and they outperform CCP estimators that employ discretization, leading to a 4 to 11-fold reduction in average mean squared error across the seven structural parameters of interest.
Second, we test our algorithms for dynamic discrete games based on a firm entry game similar to that outlined in Aguirregabiria and Mira (2007). We use the linear semigradient method here and, as before, our estimates are closely centered around the true parameters. Since this approach requires the selection of a set of basis functions for the functional approximations, we present results for different sets of basis functions (a second, third and fourth order polynomial) in this model. Our findings suggest that the choice of basis functions has only a small effect on the performance of the estimator.
Moreover, we show that a simple cross-validation procedure may be used to select the preferred set of functions.
1.1. Related literature. The paper contributes to the literature on estimation of DDC models. Rust (1987) is the seminal work in this literature. Motivated by computational considerations, Hotz and Miller (1993) propose the CCP algorithm. The CCP idea has subsequently been refined by Hotz et al. (1994) who suggest a simulation-based method, and Aguirregabiria and Mira (2002) who develop a pseudo-likelihood estimator. Arcidiacono and Miller (2011) introduce and exploit the property of finite dependence to speed up CCP estimation. Despite these advances, the estimation of DDC models remains constrained by its computational complexity, particularly in the large class of models where finite dependence does not hold. Estimation of dynamic discrete games is particularly affected by these issues as the strategic interaction of agents means that the state space increases exponentially with the number of players. It is also uncommon for finite dependence to hold in dynamic games.
The standard CCP algorithm is a two-step method, and is known to suffer from severe bias in finite samples. Aguirregabiria and Mira (2002;2007) address this issue by presenting a recursive CCP algorithm, the nested pseudo-likelihood (NPL) algorithm.
Under discrete states, the first-step estimation does not affect the rate of convergence, but shows up in form of higher-order bias for the structural parameters. The NPL algorithm can ameliorate this, see Kasahara and Shimotsu (2008). With continuous states, however, estimation of the transition density introduces bias that is the dominant term in determining the rate of convergence. This motivates the construction of our locally robust estimator which gets rid of this bias and restores parametric rates of convergence. Ackerberg et al. (2014) and Chernozhukov et al. (2022) consider semi-parametric estimation using ML methods when either finite dependence or a "terminal action" property holds (Hotz and Miller, 1993). Chernozhukov et al. (2022) also derive locally robust corrections for this setting. In both cases, the PMLE criterion can be written as a function of choice probabilities only (transition densities are not required); the authors employ non-parametric estimates for choice probabilities and correct for this estimation in the second stage. Computation and estimation is thus relatively simpler under finite dependence. By contrast, our methods are applicable to the more general and difficult setting where finite dependence may not apply. Nevertheless, the computational speed of our linear semi-gradient procedure is comparable to methods that exploit finite dependence.
For dynamic games, Semenova (2018) allows for high-dimensional state spaces, but the parameters are only partially identified.
In making use of TD learning, our methods relate to the literature on RL, particularly batch RL. Batch RL describes learning about how to map states into actions so as to maximize an expected payoff, using a fixed set of data (also called a batch); see Lange et al. (2012) for a survey. 3 It is closely related to the idea of 'experience replay', an important ingredient of RL algorithms that achieve human level play in Atari games (Mnih et al., 2015). A key step in RL, including batch RL, is the estimation of value functions. TD learning methods, first formulated by Sutton (1988), are the most commonly used set of algorithms for this purpose. We study non-parametric estimation of value functions using two TD methods: semi-gradients and AVI. Our analysis builds on the techniques developed by Tsitsiklis and Van Roy (1997) for linear semi-gradients, and Munos and Szepesvári (2008) for AVI. While Tsitsiklis and Van Roy (1997) focus on online learning (i.e., where collection of data and estimation of value functions is conducted simultaneously), we translate their methods to batch learning. With regards to Munos and Szepesvári (2008), we differ in employing assumptions that are more common to econometrics and our characterization of the rates is also different (compare Theorem 2 in their paper with our Theorem 3).
TD methods are distinct from other value function approximation methods developed in economics, e.g., parametric policy iteration (Hall et al., 2000), simulation and interpolation (Keane and Wolpin 1994), and sieve value function iteration (Arcidiacono et al., 2013). The last of these is similar in spirit to AVI with linear functional approximations. However, our semi-gradient method provides a linear approximation in a single step without any need for iterations, and we analyze AVI under generic machine learning methods. Our approximation results, and the technical arguments leading to them, are thus different from Arcidiacono et al. (2013); in fact, their setting is different too as the authors focus on estimating the 'optimal' value function, while the recursive terms in our setting are more akin to a value function under a fixed policy.
The remainder of this paper is organized as follows. Section 2 outlines the setup of the DDC model and fixes notation. Section 3 describes our TD estimation method for the functional approximations of the value functions. Section 4 proves its theoretical properties and describes the second-step estimation of the structural parameters under discrete and continuous state variables. Section 5 discusses the estimation of dynamic discrete games. Section 6 provides Monte Carlo simulations for our algorithm. Section 7 concludes. The Online Appendix discusses extensions to permanent unobserved heterogeneity.

Setup
We start with a single agent DDC model. In particular, we consider an infinite horizon model in discrete time with i = 1, . . . , n agents. We assume that the individuals are homogeneous, relegating extensions for unobserved heterogeneity to Online Appendix B.6. In each period, an agent chooses among A mutually exclusive actions, each of which is denoted by a. The payoff from the action depends on the current state x. Choosing action a when the state is x gives the agent an instantaneous utility of z(a, x) θ + e, where z(a, x) is some known vector valued function of a, x and e is an idiosyncratic error term. We denote the realization of the state of an individual i at time t by x it , and her corresponding action and error terms by a it and e it . We assume that e it is an iid draw from some known distribution g e (·). Let (a , x ) denote the one-period ahead random variables immediately following the actions and states (a, x), where x ∼ K(·|a, x), with K(·|a, x) denoting the transition density given a, x (more precisely, it is the Markov kernel). We do not make any parametric assumptions about K(·|a, x). The utility from future periods is discounted by β.
The econometrician observes a panel consisting of state-action pairs for all individuals, a i1 ), . . . , (x iT , a iT )}, for T periods (note, however, that the agent maximizes an infinite horizon objective, not a fixed T one). Typically T n in applications, so we work within an asymptotic regime where n → ∞ but T is fixed. Using this data, the econometrician aims to recover the structural parameters θ * .
In this paper, we study the CCP approach for estimating θ * (Hotz and Miller, 1993).
CCP methods are based on the conditional choice probabilities of choosing action a given state x. We denote these by P t (a|x) for a given period t but henceforth drop the subscript t with the idea that it can be made a part of the state variable x, if needed (we should also add that some of our theoretical results are based on assuming stationarity, i.e P t (a|x) is independent of t). Denote e(a, x) as expected value of the idiosyncratic error term e given that action a was chosen. Hotz and Miller (1993) show that if the distribution of e follows a Generalized Extreme Value (GEV) distribution, it is possible to express e(a, x) as a function of the choice probabilities P (a|x), i.e e(a, x) = G(P (a|x)). We assume that e follows a Type I Extreme Value distribution, which is perhaps the most common choice in applications. In this case e(a, x) = γ − ln P (a|x), where γ is the Euler constant.
The standard procedure in the CCP approach is as follows: Under the given distributional assumptions, the parameters are obtained as the maximizers of the pseudolikelihood function where for any (a, x), h(.) and g(.) solve the following recursive expressions: Here, E[.] denotes the expectation over the distribution of (a , x ) conditional on (a, x).
Note that E [h(a , x )|a, x] is a function of K(·|a, x), P (.|x). Both h(a, x) and g(a, x) have a 'value-function' form, which turns out to be useful for our approach.
Clearly, h(.) and g(.) are functions of K(·|·) and P (.|.). Since the latter are unknown, the current literature proceeds by first estimating these as (K,P ). Typically,K is obtained by MLE based on a parametric form of K(x |a, x; θ f ), whileP is estimated nonparametrically using either a blocking scheme or kernel regression. Then, given (K,P ), the values of h(.) and g(.) can be estimated by solving the recursive equations (2.1).
In the next section, we propose an alternative algorithm for maximizing Q(θ) that directly estimates h(·) and g(·) in a single step without requiring any knowledge about or estimation of K(·|·).
Notation. We assume that the distribution of (a it , x it ) is time stationary. This greatly simplifies our notation. It is not necessary for our results on the approximation properties of our TD methods, see Appendix A, but we do require it for deriving a locally robust estimator. Let P denote the stationary population (i.e, in the limit as n → ∞) distribution of (a, x, a , x ), and E[·] the corresponding expectation over P. Define E n [·] as the expectation over the empirical distribution, P n , of (a, x, a , x ). In particular, i.e we always drop the last time period in the summation index even if f (·) does not depend on a , x .
Let H denote the space of all square integrable functions over the domain A × X of (a, x). Define the pseudo-norm · 2 over H as f 2 := E[|f (a, x)| 2 ] 1/2 for all f ∈ H. We use |·| to denote the usual Euclidean norm on a Euclidean space.

Temporal-difference estimation
This section presents our TD method for estimating h(.) and g (.). Note that h(·) is a vector, of the same dimension as θ. Our methods provide functional approximations separately for each component h (j) of h. To simplify notation, we drop the superscript j indexing the elements of h(.) and proceed as if the latter, and therefore θ * , is a scalar.
However, it should be taken as implicit that all our results hold for general h(.), as long as each of its elements are treated separately.
For any candidate function, f (a, x), for h(a, x), denote the TD error by and the dynamic programming operator by Clearly, h(a, x) is the unique fixed point of Γ z [·]. TD estimation involves approximating h(a, x) using a functional class F, where each element h(·; ω) of F is indexed by a finitedimensional vector ω. The aim in TD estimation is to ostensibly minimize the meansquared TD error However, this minimization problem is neither computationally feasible nor is it proven to converge when h / ∈ F. Instead, two approaches are commonly used.
The first approach, the semi-gradient method, involves updating ω using the semi- for some small value of γ. As the name suggests, the above is not a complete gradient as the derivative does not take into account how ω affects the 'target', i.e., the future value h(a , x ; ω). Nevertheless, for linear functional classes F, it is possible to explicitly characterize the limit point of the updates, ω * , and compute it directly. Section 3.1 describes this in greater detail. In the RL literature, it is popular to employ semigradients with neural networks as the functional class F, but it appears difficult to extend our theoretical analysis to this setting.
The second approach, Approximate Value Iteration (AVI; Munos and Szepesvári, 2008), employs the idea of 'target networks'. In this approach, the parameters in the future value of h, i.e h (a , x ; ω), are fixed at the current ω, and the functional parameters iteratively updated as Clearly, the semi-gradient method and AVI are closely related: if one were to solve the minimization problem in (3.2) using gradient descent, the updates (within each iteration) would look similar to (3.1) except for fixing the value of ω in h(a , x ; ω) at the past values.
After the gradient updates converge, i.e at the end of the iteration, h(a , x ; ω) is revised with the new ω. The semi-gradient approach can thus be considered a one-step variant of AVI. Section 3.2 describes AVI in greater detail. We characterize the theoretical properties of AVI under general functional classes F including neural networks, random forests, LASSO etc.
The approximation to g follows similarly after replacing 3.1. Semi-gradients. The semi-gradient approach with linear F is particularly well suited for computation. Let φ(a, x) consist of a set of basis functions over the domain (a, x). Then the linear approximation class is F ≡ {φ(a, x) ω : ω ∈ R k }. Denote the projection operator onto F by P φ : For linear basis functions, it can be shown, e.g Tsitsiklis and Van Roy (1997), that the sequence of functional approximations h(a, x; ω j ) := φ(a, x) ω j converges to h * := φ(a, x) ω * , defined as the fixed point of the projected dynamic programming operator . Based on this characterization, we show in Lemma 1 in Usingω, we obtain an estimate of h(·) asĥ(a, x) = φ(a, x) ω.
We now turn to the estimation of g(·). As with h(·), we approximate g(·) using basis functions r(a, x), which may generally be different from φ(a, x). Let P r denote the projection operator onto the space {r(a, x) ξ : ξ ∈ R k }. The limit of the semi-gradient iterations is g * (a, x) := r(a, x) ξ * , defined as the fixed point of P r Γ e [·]. Assuming e(a, x) is known, we obtain the following characterization of ξ * in analogy with (3.3): (3.5) In the above, e(a, x) = γ − ln P (a|x) is a function of choice probabilities, which are unknown. Denote η(a, x) := P (a|x). Suppose that we have access to a non-parametric estimatorη of η. This can be obtained, e.g., through series or kernel regression. We can then plug in this estimate to obtain e(a, x;η) := γ − lnη(a, x). This in turn enables us to estimate ξ * usingξ, computed aŝ (3.6) Using the above, we obtain an estimate of g(·) asĝ(a, x) = r(a, x) ξ .
Interestingly, estimation of ξ * is unaffected to a first order by the estimation ofη, even though the latter converges to the true η at non-parametric rates (see Section 4 for a formal statement). This is because an orthogonality property holds for the estimation of ξ, in that where ∂ η · denotes the Fréchet derivative with respect to η. To show (3.7), expand Note that computation ofω andξ only involves solving linear equations of dimension dim(φ) and dim(r), respectively. This is computationally very cheap. Usingĥ(a, x) and g(a, x), we can in turn estimate θ * in many different ways. For instance, we can plug them into the PMLE estimator . (3.9) However, such plug-in estimates are sub-optimal. In Section 4.2, we suggest a locally robust version of (3.9).
Suppose that the underlying states and actions are discrete, and that our algorithm uses basis functions comprised of the set of all discrete elements of x, a. (3.10) The process can be started with an arbitrary initialization, e.g.,ĥ 1 (a, x) = z(a, x). The maximum number of iterations, J, is only limited by computational feasibility.
The minimization problem in (3.10) is equivalent to a prediction problem using the functional class F, where the outcomes are z(a, x) + βĥ j (a , x ) evaluated at the various sample draws of (a, x, a , x ). Hence, the estimation target forĥ j+1 is the conditional ). In other words, eachĥ j+1 is a nonparametric approximation to Γ z [ĥ j ], and in this manner AVI builds a series of approximate value function iterations (as its name suggests).
The interpretation of (3.10) as a prediction problem also enables us to employ any machine learning method devised for prediction, including (but not limited to) LASSO, random forests and neural networks. Our theoretical results show that it is possible to estimate h at suitably fast rates under very weak assumptions on the non-parametric estimation rates of machine learning methods.
The estimation procedure for g(·) is similar: we construct a sequence of approximations As in Section 3.1, it will be shown that the estimation error of η is first-order ignorable for the estimation of g. Usingĥ(a, x) andĝ(a, x), we can, as before, estimate θ * in many different ways, including the PMLE estimator 3.9.
Compared to the semi-gradient approach, AVI is computationally more expensive as it requires solving J prediction problems (in Section 4.1 we show that in the worst case J ≈ ln n, but this can be substantially reduced through good initializations). However, semi-gradient methods require differentiable classes of functions (e.g., random forests are not allowed) and it appears difficult to characterize their theoretical properties in the case of nonlinear basis functions.

Tuning parameters.
Both the semi-gradient and AVI methods will require choosing tuning parameters. For AVI this is straightforward: as each iteration is a nonparametric estimation problem, the tuning parameters can be chosen in the usual manner, e.g., using cross-validation. In the case of linear semi-gradient methods, the tuning parameters are the dimensions k φ = dim(φ) and k r = dim(r) of the basis functions. In analogy with AVI, we propose selecting both through a procedure akin to cross-validation.
The value of ω is estimated using a training sample and its performance evaluated on a hold-out or test sample, where the performance is measured in terms of the empirical The values of k φ , k r are chosen to minimize the mean squared TD error.

Nonlinear utility functions.
We take the utility function to be linear in θ for simplicity -and because it is the most common choice in practice -but this is easily relaxed. Denote the nonlinear utility by z(a, x; θ). Then, θ * is the maximizer of the pseudo-likelihood criterion where, for each θ, h(.; θ) is defined analogously to (2.1) with z(a, x; θ) replacing z(a, x).
We can use our TD procedures to obtain a functional approximationĥ(a, x; θ) for h(a, x; θ) at each θ (estimation of g(·) remains unchanged as it does not depend on θ). As argued earlier, computation ofĥ(·; θ) can be very fast. Appealingly, the matrix ] employed in the linear semi-gradient estimate (3.4) does not feature z(a, x; θ), and therefore only has to be inverted once. We can then plug in the values ofĥ(.; θ) andĝ(·) to estimate θ * aŝ .
A locally robust counterpart toQ(θ) can be derived in the same manner as in Section 4.2. However, computation ofθ is more involved as it is no longer convex optimization.
3.5. Unobserved heterogeneity. In Online Appendix B.6 we incorporate permanent unobserved heterogeneity into our models by pairing our TD methods with the sequential Expectation-Maximization (EM) algorithm (Arcidiacono and Jones, 2003). The resulting algorithm can handle discrete heterogeneity in both individual utilities and transition densities. Online Appendix B.7.2 provides Monte-Carlo evidence suggesting that the algorithm works well in practice.

Theoretical Properties of TD estimators
4.1. Estimation of non-parametric terms. We characterize rates of convergence for estimation of h(·) and g(·) under both semi-gradients and AVI.
4.1.1. Linear semi-gradients. We impose the following assumptions for estimation of h(·).

Assumption 1. (i) The basis vector φ(a, x) is linearly independent
Assumption 1(i) rules out multi-collinearity in the basis functions. This is easily satisfied. Assumption 1(ii) ensures that the basis functions are bounded. This is again a mild requirement and is easily satisfied if either the domain of (a, x) is compact, or the basis functions are chosen appropriately (e.g., a Fourier basis). Assumption 1(iii) is a standard condition on the rate of approximation of h(a, x) using a basis approximation.
The value of α is related to the smoothness of h(·). Newey (1997) shows that for splines and power series, α = r/d, where r is the number of continuous derivatives of h(a, ·) and d is the dimension of x. Similar results can also be derived for other approximating functions such as Fourier series, wavelets and Bernstein polynomials. The smoothness for each a ∈ A. Assumption 1(iv) requires z(a, x) to be bounded. Finally, Assumption 1(v) specifies the rate at which the dimension of the basis functions are allowed to grow.
The rate requirements are mild, and are the same as those employed for standard series estimation. For the theoretical properties, the exact rate of k φ is not relevant up to a first order since we propose estimators of θ * that are locally robust to estimation of h(·).
We then have the following theorem on the estimation of h(a, x):
We prove Theorem 1 in Appendix A.1 by adapting the results of Tsitsiklis and Van Roy (1997). The first part of Theorem 1 assures that both population and empirical TD fixed points exist. The second and third parts of Theorem 1 imply that the approximation bias and MSE of linear semi-gradients are analogous to those of standard series estimation apart from a (1 − β) −1 factor.
For the estimation ofξ we make use of cross-fitting as a technical device to obtain easyto-verify assumptions for the estimation ofη. This entails the following: we randomly partition the data into two folds. We estimateξ separately for each fold usingη estimated from the opposite fold. The final estimate of ξ * is the weighted average ofξ from both the folds. We think of cross-fitting in this context as a convenient assumption for the proofs, and do not believe it is necessary in practice.
We impose the following assumptions for the estimation of g(a, x). Let k r denote the dimension of r(a, x).

Assumption 2. (i) The basis vector r(a, x) is linearly independent, and the eigenvalues of E[r(a, x)r(a, x)
] are uniformly bounded away from zero for all k r .
(v) k r → ∞ and k 2 r /n → 0 as n → ∞. (vi)ξ is estimated from a cross-fitting procedure described above. The conditional choice Assumption 2 is a direct analogue of Assumption 1, except for the last part which provides regularity conditions when η(·) is estimated. These conditions are typical for locally robust estimates and only require the non-parametric function η(a, x) to be estimable at faster than n −1/4 rates. This is easily verified for most non-parametric estimation methods such as kernel or series regression. Under these assumptions, we have the following analogue of Theorem 1, which we prove in Appendix A.2 .
The L 2 error for the difference between g(a, x) and r(a, x) ξ is bounded as Iterating the above gives (4.1) Equation (4.1) can be considered a special case of error propagation (Munos and Szepesvári, 2008).
Recall thatĥ 1 is an arbitrary initialization. It is thus straightforward to provide conditions under which h −ĥ 1 2 is bounded by some constant M 1 . As for the second term in (4.1), recall from the discussion in Section 3.2 that the minimization problem (3.10) corresponds to non-parametric estimation of Γ z [ĥ j−1 ] using the functional class F. Most machine learning methods come with guarantees on the non-parametric estimation rate We now describe our assumptions for AVI. Let X denote the d-dimensional space of x, and define W γ,∞ (X ) as the Hölder ball with smoothness parameter γ: Assumption 3(i) is a standard requirement in non-parametric estimation. The assumption of γ-Hölder continuity is taken from Farrell et al. (2021). Assumption 3(ii) is a mild condition on the initializationĥ 1 . Assumption 3(iii), which is novel to this paper, is a crucial smoothness condition requiring the operator Γ z [·](·, a) to map all bounded f onto W γ,∞ (X ). In Online Appendix B.2, we show that this is satisfied if z(a, ·) and K(x |a, ·) are γ-Hölder continuous.
Assumption 3(iv) is a high-level condition on the machine learning (ML) method F.
The requirement of bounded f implies that the ML method cannot diverge in the l ∞ sense, see Farrell et al. (2021) for a discussion of this in the context of multi-layer perceptrons (MLPs). The second part of Assumption 3(iv) implies that the ML method is able to non-parametrically approximate all functions in W γ,∞ (X ) at the rate of at least n −c .
Most ML methods are proven to satisfy this. Consider, for instance, the class F of MLPs of width W and depth L; MLPs and, more generally, Neural Networks are widely used in RL. The results of Farrell et al. (2021) imply that for W n d 2(γ+d) ln 2 n and L ln n, Thus, Assumption 3(iv) is satisfied for MLPs. See Biau (2012) for related results on random forests.
Assumptions 3(iii) and 3(iv) imply that one can estimate Γ Combined with (4.1), this proves: Theorem 3. Suppose Assumptions 3(i) to 3(iv) hold. Then, for all n large enough, The first term in the expression for E h −ĥ J 2 from Theorem 3 is the statistical rate of estimation of h. The second term is the numerical error, which is seen to decline exponentially with the number of iterations J. Setting J ln n will ensure the numerical error is smaller than the statistical rate of convergence. The number of iterations can be further reduced using a good initialization that makes M 1 small. For instance,ĥ 1 can be the linear semi-gradient estimator. This is fast to compute and ensures M 1 = o p (1).
Incidentally, Theorem 3 justifies the use of neural networks for batch RL; to the best of our knowledge this appears to be new even in the RL literature.
Turning to estimation ofĝ, we again assume cross-fitting is employed as in Theorem 2, i.e.,η is computed from one half of the data, andĝ is computed using AVI on the other half, takingη as given.

Assumption 4. (i) Assumptions 3(i) -3(iv) hold after replacing
(ii)ĝ is estimated from a cross-fitting procedure. The conditional choice probability Assumption 4(ii) is similar to Assumption 2(vi) with the additional requirement that η andη be α-Hölder continuous, the latter with probability approaching one. Note that η would be α-Hölder continuous if the non-parametric estimator consistently estimates not only η but also its first α derivatives.
Sinceη(a, ·) ∈ W γ,∞ (X ), the first term, g −ĝ J 2 , can be bounded by a similar rate as in Theorem 3 (recall that we takeη as given under cross-fitting). As for the second term, g −g 2 , observe that under Assumption 4, where the ' ' holds with probability approaching one, and uses see the discussion following (3.7). Since η −η 2 2 = o p (n −1/2 ), the above proves: Theorem 4. Suppose Assumptions 4(i) to 4(ii) hold. Then, with probability approaching one,

Estimation of structural parameters. Estimation of h(a, x) and g(a, x) is in-
herently non-parametric. This is because h(a, x) and g(a, x) are functions of two nonparametric terms: the choice probabilities η(a, x), and the transition densities K(x |a, x).
The TD estimators implicitly take both into account. Under discrete states, the firststep estimation error does not affect the rates of convergence of structural parameters (see Online Appendix B.1). With continuous states, however, the first-step non-parametric estimation error does affect the estimation of θ * to a first order when using the PMLE criterion (3.9). This is because the estimates for K(x |a, x) and θ * are not orthogonal under a PMLE, which extends to the lack of orthogonality between the estimates for h(a, x), g(a, x) and θ * . Consequently, the PMLE estimator with plug-in values ofĥ and g will converge at slower than parametric rates.
We can recover √ n-consistent estimation by adjusting the PMLE criterion to account for the first-stage estimation of h and g. Denote (ã,x) := (a, x, a , x ) and m(a, x; θ, h, g) := The PMLE estimator with plug-in estimates solves E n [m(a, x; θ,ĥ,ĝ)] = 0, but this is not robust to estimation of h, g. Let denote the continuation value given (a, x). Also, define λ(a, x; θ) as the fixed point of the 'backward' dynamic programming operator where (a − , x − ) denotes the past actions and states preceding (a, x), and In Online Appendix B.3, we show that the locally robust moment corresponding to m(a, x; θ, h, g) is given by The construction of the locally robust moment (4.4) is new. But it is infeasible since λ(·), h(·) and g(·) are unknown. However, we can replace these quantities with consistent estimates. We have already described how to estimate h(·), g(·). Letθ denote the plug-in estimator of θ * using (3.9); note thatθ consistently estimates θ * but is not efficient. An estimator,λ(·), of λ(·) can then be obtained by applying either of our TD estimation methods on (4.2), withθ,ĥ,ĝ plugged in in place of θ, h, g. For instance, using AVI, we could obtain iterative approximations {λ (j) , j = 1, . . . , J} for λ(·) usinĝ Plugging inλ(·),ĥ,ĝ into (4.4), we obtain the feasible locally robust moment Based on the above, we can obtain a locally robust estimator,θ, as the solution to E n [ζ n (ã,x; θ,ĥ,ĝ)] = 0. We recommend obtaining this estimate using cross-fitting, see Section 4.2.1 for more details. Compared to the plug-in estimate (3.9), our locally robust estimator requires computation of λ(·), but even this can be avoided when linear semi-gradients are used for estimating h, g (see Online Appendix B.4). Solving E n [ζ n (ã,x; θ,ĥ,ĝ)] = 0 is also computationally easy; the correction term is a constant, and ∂ θ ζ n (ã,x; θ,ĥ,ĝ) = ∂ θ m(a, x; θ,ĥ,ĝ) is negative definite (as the PMLE criterion is concave), so solving this is no harder than solving the original moment condition without a correction term.

4.2.1.
√ n-consistent estimation. We focus on the general construction of the locally robust estimator,θ, using (4.5). As mentioned in the previous sub-section, we advocate cross-fitting to obtain this estimator. This entails computingθ,λ,ĥ,ĝ using one half of the sample, say N 2 , and plugging them into the locally robust moment to computê is the empirical expectation using only observations from the other half of the sample N 1 . Following the analysis of Chernozhukov et al. (2022), it can be shown that this estimator has the same limiting distribution as the one based on (4.4). In particular, it achieves parametric rates of convergence. We state the regularity conditions below (for the remainder of this section we allow θ * to be vector valued): Assumption 5(i) implies θ * is identified. Assumption 5(ii) is a mild regularity condition that is similar to Assumption 5 in Chernozhukov et al. (2022). We refer to Chernozhukov et al. (2018) for more details on three-way splitting.
We are now ready to state the main result of this section.
Theorem 5. Suppose that either Assumptions 1, 2 & 5 or 3-5 hold. Then the estimator, θ of θ * , based on (4.5) is √ n-consistent, and satisfies The proof of the above theorem follows by verifying the regularity conditions of Chernozhukov et al. (2022, Theorem 9). Since these are more or less straightforward to verify given our previous results, we omit the details. For inference onθ, the covariance matrix Chernozhukov et al. (2022) provide conditions under whichV is consistent for V ; these are straightforward to translate to our context but we omit them for brevity. Alternatively, one could employ the bootstrap.

Estimation of dynamic discrete games
So far we have considered applications of our algorithm to single-agent models, where we have argued that there are substantial computational and statistical gains from using our procedure. These gains are magnified when extended to the estimation of dynamic discrete games. As in Section 3, we take θ * to be scalar to simplify the notation; all our results continue to hold for vector valued θ, as long as each dimension is treated separately. Evolution of the states in the next period is determined by the transition density K(x |a, x) where a : = (a 1 , . . . , a N ) denotes the actions of all the players. We denote by x tm the state at market m in time period t, by a tm the vector of actions by all players at time t in market m, and by a itm the action of player i at time t in market m. We also let P i (a i |x t ) denote the choice probability of player i taking action a i when the state is x t , and define e(a i , x) := γ − ln P i (a i |x).
As in the single agent case, the parameters θ * can be obtained as solutions to the pseudo-likelihood function: where h i (.) and g i (.) are now player-specific, and given by In contrast to (2.1) in the single agent case, the expectation now averages over the actions of the other players as well.
Previous literature estimates θ * using a two-step procedure: In the first step, the conditional choice probabilities P i (a i |x t ) are calculated non-parametrically. These, along with estimates of K(.) are then used to recursively solve for h i (.) and g i (.) using equation and for any function f (·), we define x tm , a t+1m , x t+1m ).
( 5.4) Similarly, the AVI iterations for h i (·), g i (·) are given bŷ ( 5.5) If the players are symmetric (z i (a i , a −i , x) does not depend on player i) we can obtain computationally faster and more precise estimates by pooling across players.
Importantly, neither of the estimation strategies (5.3) nor (5.5) require partialling out other players' actions, leading to a tremendous reduction of computation. Intuitively, the procedures take expectations over other players' actions 'internally' using the empirical distribution. The non-parametric estimatesĥ i ,ĝ i can be plugged into the PMLE criterion (5.1) to obtain an estimate for θ aŝ .
It is straightforward to construct a locally robust estimator for θ in analogy with that for single-agent models. We describe this in Online Appendix B.5.
By the same reasoning as in Online Appendix B.1, it possible to show that with discrete states, h i (.) and g i (.) are numerically identical to the estimates obtained by plugging in cell estimatesP j (·|x) andK(.) in (5.2). This implies the psuedo-likelihood with plugin estimates for h(.) and g(.) is not efficient even with discrete states, as discussed by Aguirregabiria and Mira (2007). However the values of h(.) and g(.) can be plugged into other, more efficient objectives, such as our locally robust estimator or the minimum distance estimator of Pesendorfer and Schmidt-Dengler (2008). With continuous states, one would need to employ locally robust corrections even for the minimum distance estimator to recover parametric rates of convergence for θ. The locally robust correction term can be constructed in a similar way as that for the PMLE criterion.

Simulations
We run Monte Carlo simulations to test our estimation methods. We start by presenting results for a DDC model using both the linear semi-gradient approach, and the AVI approach with random forests as the prediction method. For each of these methods, we compare findings from the locally robust version of our estimators to those obtained without correction. Our simulations for this part are based on the firm entry problem described in Aguirregabiria and Magesan (2018).
In a second set of Monte Carlo simulations, we test our estimation method for dynamic discrete games. For this part, we present results from the linear semi-gradient approach using different sets of basis functions to approximate the value function terms, and employ the cross-validation method described in Section 3.3 to select the preferred set of basis functions. Our simulations are based on the dynamic firm entry game used in Aguirregabiria and Mira (2007).
Finally, we run additional simulations based on the famous Rust (1987) bus engine replacement problem (see Online Appendix B.7.2). Using this model, we also provide simulation results for a case with permanent unobserved heterogeneity.
6.1. Firm entry problem. Consider the following dynamic firm entry problem described in Aguirregabiria and Magesan (2018). A firm decides whether to enter (a t = 1) or not enter (a t = 0) in a market for t = 1, ..., T time periods. The payoff when entering is given by Π where V P t , F C t and EC t denote the firm's variable profit, fixed cost and entry cost, and ε t is a transitory shock that follows a logistic distribution. Variable profit is given by where ω t denotes the firm's productivity shock, and z 1t , z 2t are exogenous state variables affecting the price-cost margin in the market. The fixed cost is given by F C t = θ F C 0 + θ F C 1 z 3t , and the entry cost is given by EC t = (θ EC 0 + θ EC 1 z 4t )(1 − a t−1 ), where z 3t , z 4t are further exogenous state variables, and a t−1 denotes the entry decision in period t − 1 which is an endogenous state variable. The payoff of not entering is normalized to zero. The param- are the structural parameters of interest.
To carry out the simulations, we choose values for the structural parameters θ * (θ V P for the autoregressive processes of z jt and ω t (γ j 0 = 0.0, γ j 1 = 0.6, γ ω 0 = 0.2, γ ω 1 = 0.6), and discretize the exogenous state variables to obtain a transition matrix with a 6-point support following Tauchen (1986). The resulting dimension of the state space is 2 × 6 5 = 15, 552. The discretization of the support is for simulations only; our estimation algorithms treat these variables as continuous and do not require any prior knowledge of how they evolve (the knowledge of AR(1) dynamics is also not used). We iterate on the value function to obtain the vector of choice probabilities for each combination of the states, and use these to derive the steady-state distribution of the state variables. Using the steady-state distributions, we generate data for 3, 000 firms, with T = 2 time periods. Table 1 shows results for 1000 simulations using the linear semi-gradient method. Panel B of Table 1  We present results generated with and without the locally robust correction. For the results without correction, we obtain estimates for θ * using (3.9). To generate the locally robust estimates, we use moment equation (B.11) for the linear semi-gradient method, and moment equation (4.5) for the AVI method where we employ a random forest to derive an estimate for the λ(a, x,θ) term contained in the locally robust moment. As before, the AVI method for estimation of λ(·) iterated 70 times. We also use the sample 4 For the ω's relating to parameters θ V P 0 , θ V P 1 , θ V P 2 , θ F C 0 , θ F C 1 , and for ξ, the terms include a constant, the exogenous state variables and their interactions up to a second order, the player's binary choice a t and the interactions of a t with all terms in the exogenous states. Given the set-up of the model, we treat the interactions z 1t exp(ω t ) and z 2t exp(ω t ) instead of z 1t and z 2t themselves as state variables. In addition to the terms included above, the ω's relating to parameters θ EC 0 and θ EC 1 also contain the terms (1 − a t−1 ) and (1 − a t−1 )z 4t , respectively. The total number of terms included is 42 (43 for θ EC 0 and θ EC 1 ).

Simulation results -firm entry problem. Panel A of
splitting method described in Section 4.2.1 for the locally robust estimators, and we obtain the finalθ as weighted average of the θ * estimates from the two samples. 5 Panel A of Table 1 shows that our linear semi-gradient results are closely centered around the true values. While the locally robust estimator should in theory be preferable, we find that it produces results which are similar and if anything have slightly higher mean squared error than the non-robust version of our algorithm. In fact, there is very little bias, and the distribution of the estimates under the non-robust version is already very close to normal, see Online Appendix B.7 for the plots of the finite sample distributions.
On the other hand, locally robust methods are associated with higher variance due to sample splitting. So there appears to be no gain from the locally robust method here.
Panel B of Table 1 shows that the AVI method produces estimates with slightly higher small sample bias than the linear semi-gradient. In this case, the locally robust method is more useful as it decreases the bias significantly, especially for estimation of the entry cost parameters.
6.1.2. Comparison with existing methods. We compare our findings from Section 6.1 to those that would be obtained with a standard CCP estimator where the state variables are discretized and the transition and choice probabilities are estimated using cell values.
We start our simulations by generating data as outlined in Section 6.1. We then discretize the state space by creating dummy variables for each state variable z 1t , z 2t , z 3t , z 4t and exp(ω t ). To make the estimation feasible, the state space needs to be restricted further. A common approach is to use K-means clustering here, but this is not appropriate in the given simulation setting where the state variables are independent by construction. We therefore restrict the state space grid by combining variables z 1t and z 2t into a binary variable taking value one whenever both individual dummies take value one. The resulting state space consists of four binary variables, implying 16 cells in the state space grid. We try alternative feasible ways of discretizing the state space, but find that these do not lead to improvements over the chosen method. We run 1000 simulations, and the results are shown in Table 2.
It can be seen that, compared to the results from Table 1, the CCP estimator leads to substantially larger bias in some of the estimated parameters. Column (4) shows that the  (2) and (4) report the empirical mean and standard deviation for the estimated parameters. Columns (2)-(3) are based on the estimation method without correction function, columns (4)-(5) report results using the locally robust estimator. The mean squared errors are reported in columns (3) and (5), respectively.
corresponding mean squared errors are large and generally exceed those obtained using our estimators. This is particularly true for parameters θ V P 1 and θ V P 2 , highlighting the challenges inherent in the discretization of continuous state spaces. Overall, the average mean squared error increases from 0.01 − 0.04 across all parameters in Table 1 to 0.16 in Table 2. These results show that our estimation methods lead to important improvements over existing methods in models with as few as five continuous state variables.
6.2. Firm entry game. Consider the following firm market entry game, which is similar to that described in Aguirregabiria and Mira (2007 where ln(S tm ) is a measure of consumer market size of market m in period t, and ε itm is a transitory shock that follows a logistic distribution. We assume that ln(S tm ) is continuous and follows an AR(1) process, where the parameters are the same across markets: ln(S tm ) = α + λ ln(S (t−1)m ) + u tm . The error term u tm is assumed to follow a normal N (0, 1) distribution. The profit of not entering is normalized to zero, and the discount factor β is 0.95. The parameters θ * ≡ {θ RS , θ RN , θ F C , θ EC } are the structural parameters of interest. The state variables in this setting are given by the current market demand variable S tm , as well as the vector of all firms' previous entry decisions a (t−1)m = a i(t−1)m : i = 1, ..., 5 .
To carry out the simulations, we choose values for the structural parameters θ * (θ RS = 1, θ RN = 1, θ F C = 1.7, θ EC = 1), and for the autoregressive process for log market size (α = 1.5, λ = 0.5). We discretize ln(S tm ) and obtain a transition matrix for the discretized variable with a 10-point support following the method by Tauchen (1986). As in the Monte Carlo experiments for the firm entry problem in Section 6.1, the discretization is for simulations of the data only and we treat the state variables as continuous in our estimations. We then solve for the Markov-Perfect-Equilibrium of the game. This is done by finding the firms' conditional value functions ν j (S tm , a (t−1)m ) for each of the 2 5 × 10 = 320 possible combinations of the state variables through repeated iteration, and using these to derive the equilibrium choice probabilities p(S tm , a (t−1)m ). Based on the equilibrium probabilities, we compute the equilibrium distribution of state variables.
Using the equilibrium distributions, we generate data for 1000 and for 3000 markets, with T = 2 time periods. For the ω s relating to parameters θ RS , θ RN , θ F C and for ξ, the terms include a constant, terms up to the second/third/fourth order in the state variables ln(S tm ) and ln(1 + j =i a j(t−1)m ), the player's binary choice a itm and the interactions of a itm with all terms in the state variables. The total number of terms is 12/20/30. In addition to these terms, the ω s relating to parameter θ EC also contain the term (1 − a i(t−1)m )a itm .
ω and ξ using equation 5.3. 7 Finally, we obtain estimates for the θ * parameters as the solutions to the pseudo-likelihood function 5.1. 8 The results are shown in Table 3. Panels A, B and C present simulations for the same dataset using different basis functions to parameterize the value function terms h (a, x) and g (a, x). Column (2) shows that even with 1000 markets our algorithm produces parameter estimates that are closely centered around the true values. The results are generally similar across Panels A to C, although the bias and mean squared error tends to be lowest for the second order polynomial, and highest for the fourth order polynomial.
This is especially the case for the parameter on the number of market entrants, θ RN .
To assess these differences formally, we use the cross-validation procedure described in Section 3.3. The procedure is applied to ten random samples of market size 1000, and we find that the TD error criterion consistently selects the second order polynomial as the optimal set of basis functions. This suggests that the proposed cross-validation method provides useful guidance when choosing a set of basis functions in practice.
In a similar version of the firm entry game, Aguirregabiria and Mira (2007) use the NPL algorithm and derive results comparable to ours. Our simulations rely on slightly larger sample sizes, but note that these differences are likely due to the fact that our algorithm implicitly estimates the transition densities for market size and therefore does not rely on the true transition density to be used in the estimation. For a direct comparison of our results with those obtained using the NPL algorithm, one would need to obtain a non-parametric estimate of the transition density which is not trivial in practice.
As expected, columns (4) and (5) show that increasing the market size generally reduces the small sample bias in the estimated parameters, and leads to a fall in the empirical standard deviations. In addition to being smaller, the mean squared error across Panels A to C is also more similar across the three sets of basis functions. As before, we employ the cross-validation method described above to compare these specifications more formally.
In line with the estimation results, we find that all three polynomials now produce very similar sets of mean squared TD error, even though the second order polynomial continues to be the one that is selected by the criterion. 9 While we view this as further evidence that 7 Given the symmetric set-up of the game, we pool the data across players in this application. 8 We run our simulations on a MacBook Pro with an M1 chip and 16 GB of RAM. The approximate computation times for one estimation round are: 1.5 seconds (second order polynomial, 1000 markets); 3.7 seconds (second order polynomial, 3000 markets); 2.5 seconds (third order polynomial, 1000 markets); 6.7 seconds (third order polynomial, 3000 markets); 4 seconds (fourth order polynomial, 1000 markets); 12 seconds (fourth order polynomial, 3000 markets). 9 As before, we compute the TD criterion for ten random samples.  h(a, x) and g(a, x). Column (1) shows the true parameter values in the model. Columns (2) and (4) report the empirical mean and standard deviations for the estimated parameters, based on a sample of 1000 and 3000 markets, respectively. The mean squared errors are reported in columns (3) and (5). All results are based on the estimation method without correction function.
the proposed cross-validation method can provide useful guidance to choose a suitable set of basis functions, more importantly, the small differences in the results across panels A to C also suggest that our methods prove fairly robust to this choice in practice.

Conclusions
We propose two new estimators for DDC models which overcome previous computational and statistical limitations by combining traditional CCP estimation approaches with the idea of TD learning from the RL literature. The first approach, linear semigradient, makes use of simple matrix inversion techniques, is computationally very cheap and therefore fast. The second approach, Approximate Value Iteration, can be easily combined with any machine learning method devised for prediction. Unlike previous estimation methods, our methods are able to easily handle continuous and/or high-dimensional state spaces in settings where a finite dependence property does not hold. This is of particular importance for the estimation of dynamic discrete games. We also propose a locally robust estimator to account for the non-parametric estimation in the first stage.
We prove the statistical properties of our estimator and show that it is consistent and converges at parametric rates. A range of Monte Carlo simulations using a dynamic firm entry problem, a dynamic firm entry game and a version of the famous Rust (1987) engine replacement problem show that the proposed algorithms work well in practice.  Control, vol. 42, no. 5, pp. 674-690, 1997.

Appendix A. Proofs of main results
In what follows we drop the functional argument (a, x) when the context is clear and x ) for different functions f .

Lemma 1.
There exists a unique fixed point to the operator P φ Γ z . If Assumption 1(i) holds, this fixed point is given by φ ω * , where ω * is such that Proof. First, note that Γ z , and therefore P φ Γ z , are both contraction maps with the contraction factor β. This implies P φ Γ z has a unique fixed point. Clearly, this fixed point must lie in the space L φ . Let us denote this as φ ω * . Now, for any function f ∈ L φ , Since φ ω * is the fixed point, we must have But φ is linearly independent and E[φφ ] −1 is non-singular, by Assumption 1(i). Hence it must be the case that E [φ (z + βφ ω * − φ ω * )] = 0.
For the next lemma, we use the following definition: a square, possibly asymmetric, if A is negative definite, it is also invertible. This holds even if the matrix is asymmetric, see e.g Johnson (1970).

Lemma 2. Under Assumption 1(i), the matrix
, and is therefore invertible.
Proof. The proof is adapted from Tsitsiklis and van Roy (1997). Recall the definition of φ ω * as the fixed point to P φ Γ z [·] from Lemma 1. We now show that Observe that since the second expression in the first equation is 0. Now, where the last equality holds since Since P φ Γ z [·] is a contraction mapping with contraction factor β, it follows In view of the above, This completes the proof of the lemma.

Lemma 3. Suppose that Assumption 1(i) holds. Then,
Proof. Recall that h(·, ·) is the unique fixed point of Γ z , and φ ω * is the unique fixed point of P φ Γ z . The operator Γ z is a contraction mapping with contraction factor β.
Furthermore, the projection operator P φ is linear, and P φ [f ] 2 ≤ f 2 for any function Rearranging the above expression proves the desired claim.
For the proofs of Theorems 1-2, we work within a more general setting than in the main text, by letting the distribution of (a it , x it ) be time-varying. Let P t denote the population distribution of (a, x) at time t. Also, let P denote the probability distribution of the process {(a 1 , x 1 ), . . . , (a T , x T )}. Note that P ≡ P 1 × · · · × P T . We will denote E[·] as the expectation over P . Furthermore, we use the o p (·) and O p (·) notations to denote convergence in probability, and bounded in probability, respectively, under the probability distribution P .
We also need to extend the definitions of P and E[·]: Let P denote the relative frequency of occurrence of (a, x, a , x ) in the data as n → ∞. Let E[·] denote the corresponding expectation over P. Note that P is different from P since the latter is the distribution of (a, x, a, x ) after dropping the time index. However, the two are related since for any . These updated definitions of P and E[·] are applicable wherever these notations are used in the main text.
Note that due to the Markov process assumption, the conditional distribution P (a t+1 , x t+1 |a t , x t ) is always independent of t (indeed, one could always consider t as also a part of x). Hence, Also note that time stationarity of (a it , x it ), if it holds, implies P t ≡ P and for all t.
The second claim follows directly from Lemma 3 and Assumption 1(iii).
To prove the last claim, we first show that with probability approaching one, . We then have Aω * = b andÂω =b.
We can combine the two equations to get The above implies with probability approaching one, for any constant c ∈ (0, 1). Given (A.3) and (A.4), with probability approaching one.

Now observe that
where the final inequality follows from the second claim of this theorem and (A.2). The last claim then follows from the above along with the fact that, by Assumption 1(iv),

Hence (A.5) follows once we shoŵ
We now prove (A.6). Denoting N 2 the set of observations in the second fold: Define the above as R 1nt + R 2nt . First consider the term R 1nt . Define Clearly, E[δ it |N 2 ] = 0. We then have Now for any (a, x, a , x ), we can note from the definition of ψ(·) that with probability approaching one, where the second inequality follows from Assumption 2(iii), and the third follows from 2(v). 10 In view of (A.7) and (A.8), there exists C < ∞ such that where the last equality follows by Assumption 2(v). This proves Next consider the term R 2nt . We note that E[ψ it (η)] is twice Fréchet differentiable.
In the main text we have shown that ∂ η E[ψ it (η)] = 0 (c.f equation (3.7)). Furthermore, following some straightforward algebra it is possible to show |∂ 2 η E[ψ it (η)]| ≤ C 1 √ k r , for some C 1 < ∞, as long as η is bounded away from 0 (as assured by Assumption 2(v)). Hence and the values ofĥ(a, x) andh(a, x) coincide exactly. Thus, the two algorithms give identical results (a similar property also holds for g(a, x)). Since our estimatesĥ(a, x) coincide with those from the standard CCP estimators, the resulting estimateθ is also exactly the same.
When the states are discrete, Aguirregabiria and Mira (2002) show that the estimation of η is orthogonal to the estimation of θ * . It is important to note, however, that the estimation of the Markov kernel K(x |a, x) is not orthogonal to the estimation of θ * . Now, with discrete states, any estimate,K(x |a, x), of K(x |a, x) converges at parametric rates, so √ n-consistent estimation of θ is still possible. However, as we show in Section 4.2, this creates issues with continuous states.

B.2. Smoothness properties of dynamic-programming operators and fixed points.
The following result provides a sufficient condition for Assumption 3(iii).
Lemma 4. Suppose that max |p|≤γ sup a,x |∂ p x z(a, x)| < C and max |p|≤γ sup a,x |∂ p Proof. For the first claim, note that Hence, the claim is satisfied for M 0 ≥ C/β.
We now turn to the second claim. For any |p| ≤ γ, we have We now provide sufficient conditions for h(·) to be γ-Hölder continuous. The result provides a sufficient condition for Assumption 3(i). It can also be used as a justification for Assumption 1(iii); see the discussion following Assumption 1. We now show h(a, ·) ∈ W γ,∞ for all a. For any f ∈ S 0 , and |p| ≤ γ, where M := max{M 0 , C + (1 − β)M 0 C}. It is straightforward to write down analogous versions of Lemmas 4 and 5 for Γ g [·] and g(a, ·). We omit the details.
B.3. Verification of the locally robust moment function. The locally robust correction needs to satisfy two properties. First, it must be mean zero. This is easily verified.
Second, it must satisfy the zero Gâteaux derivative requirement, i.e., after adding the correction term, the moment condition should satisfy for all γ ∈ H, where H is the set of all square integrable functions over the domain A×X .
We now verify this below: Let ψ h (·) denote the solution of the functional equation and similarly, ψ g (·), the solution of the functional equation Following some tedious but straightforward algebra, we can show ψ g (·) = ψ(·) and ψ h (·) = θψ(·), where ψ(·) is defined in (4.3) in the main text. Now, the doubly robust moment condition in (4.4) can be expanded as ζ (ã,x; θ, h, g) = m(a, x; θ, h, g) where we have defined λ h (a, x; θ) = θλ(a, x; θ) and λ g (a, x, ; θ) = λ(a, x; θ) for the second equality. Using the definitions of ψ(·), ψ h (·), ψ g (·) and λ(·), we observe that λ h (a, x; θ) is the fixed point of the 'backward' dynamic programming operator (B.5) and λ g (a, x; θ) is the fixed point of We now verify (B.2). To this end, observe that for any square integrable γ, where the second equality follows from the definition of ψ h (·) in (B.3). Since λ h (·) is the fixed point of Γ † h,θ [·], we can expand the third term in (B.7) as a, x; θ, h, g) We make use of (B.8) and (B.9) to construct a locally robust moment for θ * . Following Newey (1994) and Chernozhukov et al. (2022), this is given by We can now construct a locally robust estimator for θ * based on (B.10). Following Chernozhukov et al. (2022), we employ a cross-fitting procedure by randomly splitting the data into two samples N 1 and N 2 . We computeω andξ using one of the samples, say N 2 . We also compute,θ, a preliminary consistent estimator of θ * by applying the plug-in estimator (3.9) on observations in N 2 . Denote by E (1) n [·] the empirical expectation using only the observations in the first sample, N 1 . We then obtainθ as the solution to the moment equation The use of cross-fitting is critical. If we had used the entire sample to estimate all of θ * , ω * and ξ * , we would have E n [ϕ h (ã,x,ω)] = 0 and E n [ϕ g (ã,x,ξ)] = 0, which implies E n ζ n (ã,x, θ,ω,ξ) = E n m (a, x, θ,ω,ξ) . As noted by Chernozhukov et al. (2022), cross-fitting gets rid of the 'own observation bias' that is the source of the degeneracy here. The solution,θ, of (B.11) is the locally robust pseduo-MLE estimator of θ * .
Even though the estimator in (B.11) is predicated on h(x, a) and g(x, a) being truly finite dimensional, the work of Chernozhukov et al. (2022) suggests that ζ n (ã,x; θ, v) should have the same form as ζ n (ã,x; θ, h, g) from (4.5). Suppose that the same vector of basis functions is used for both h and g, i.e., φ(a, x) ≡ r(a, x). Then, substituting the expressions for ϕ h (ã,x, ω), ϕ g (ã,x, ξ), we obtain (after some algebra) The locally robust moment for θ is where Here, λ i (a i , x; θ) is the fixed point to For computation, we employ cross-fitting as in the single-agent setting and randomly split the markets into two samples N 1 and N 2 . We computeĥ i ,ĝ i ,λ i for all players using one of the samples, say N 2 . Letθ denote the plug-in estimator of θ * obtained using N 2 .
Alternatively, if the players were symmetric, i.e, z i (·) did not depend on i, we can pool all the players and compute a commonĥ,ĝ,λ.
The locally robust estimator (B.12) has the same form as (B.11), except for there being separate correction terms for the estimatesĥ i ,ĝ i of each player i. Its theoretical properties are thus equivalent to, and can be derived in the same manner, as those for single-agent models.
B.6. Incorporating permanent unobserved heterogeneity. We incorporate permanent unobserved heterogeneity into our models by pairing the techniques from Section 3 with the sequential Expectation-Maximization (EM) algorithm (Arcidiacono and Jones, 2003). The use of the sequential EM algorithm in CCP estimation under unobserved heterogeneity was first advocated by Arcidiacono and Miller (2011), and we employ a similar approach.
Suppose that in addition to the observed state x, and the choice-specific shock e, individuals also base their choice decisions on a random state variable s which is known to the individual, but unobserved to the econometrician. As is common in the literature, we assume a finite set of unobserved states indexed by {1, 2, . . . , k, ...K}. The number of states is also assumed to be known a priori. Let π k denote the population probability P (s = k). The value of s for an individual is assumed to be permanent and not change with time. However, we do not place any restrictions on the transition density K(x |a, x, s), which is allowed to change with s.
To simplify the exposition, we only describe the basic version of the algorithm without local robustness corrections as in Section 4.2. It is straightforward to incorporate the correction term into the algorithm, but it comes at the expense of higher computational times.
Suppose that the per-period utility is given by z(a, x, s)θ. For each k, define h k (a, x) and g k (a, x) as the solutions to To simplify notation, let h itk := h k (a it , x it ) and g itk := g k (a it , x it ). If these quantities were known, one can estimate (θ, π) by maximizing the integrated pseudo-likelihood (B.13) As before, we have chosen to make h(·) uni-dimensional to simplify the notation. We emphasize that the methods suggested here could be employed even if there was no heterogeneity in individual utilities, but the transition density were heterogenous across individuals. This is because even just the latter would result in heterogenous values for h(a, x), as it is a function of transition densities.
To make (B.13) usable, h k (a, x), g k (a, x) would have to be estimated. Similar to the motivation of TD methods in Section 2, the heuristic we employ is to approximate h k (at each k) by minimizing the mean-squared TD error in also taking the expectation over the distribution of the unobserved state s, and P (s = k|a, x) := P r(s = k|a 1 , x 1 , . . . , a T , x T ).
If P (s = k|a, x) were known, we could use them as weights in the semi-gradient and AVI approaches. In particular, for linear semi-gradient methods, the estimates for h k , g k and we have used the notation p ik = P (s = k|a i , x i ), z itk := z(a it , x it , k), φ it := φ(a it , x it ), r it := r(a it , x it ), andė it+1k is the current estimate (described below) of e it+1k := γ − ln P (a it |x it , s i = k). Similarly, for AVI, we could obtain iterative approximations The above involves solving a sequence of weighted non-parametric estimation problems, which most machine learning methods can easily handle. The estimatesĝ (j) k (·) of g k (·) can be obtained in a similar fashion.
Estimation of θ * , h k , g k using (B.13) and (B.14) or (B.15) requires knowledge of the unknown quantities π k and p ik along withė it+1k . Furthermore, even if π k were known, maximizing the integrated likelihood function (B.13) is computationally very expensive.
The sequential EM algorithm of Arcidiacono and Jones (2003) solves both issues. 11 To describe the procedure, let Denote byπ k andp ik the estimates for π k and p ik . The algorithm consists of two steps: the M-step and the E-step. We first describe the M-step. Here, we update the estimates for h k , g k and θ * based on the current estimates for π k , p ik and e it+1k . To this end, note that we can estimate h k , g k using either (B.14) or (B.15). From these we can in-turn updateθ asθ Next, given (θ,ĥ k ,ĝ k ), we updateπ k ,p ik andė it+1k for all i, k. This is the E-step of the EM algorithm. This step consists of three parts. In the first part, we use the current θ,ĥ k ,ĝ k andπ k to updatep ik for each i, k using Bayes' rule: .

(B.17)
11 Maximizing (B.13) is not equivalent to Full Information Maximum Likelihood (FIML). As in Arcidiacono and Jones (2003), the identification and asymptotic properties of θ, π are in fact determined by constructing moment conditions that correspond to the first order conditions from maximizing Q(θ, π), augmented with additional moments identifying h k , g k . Together, these moment conditions, which motivate the sequential EM algorithm, can in turn be related to the identification properties of FIML. In the second part, we updateπ k , for each k, aŝ Finally, we also updateė it+1k for all i, t, k aṡ e it+1k ←− γ − ln l it+1k (θ,ĥ k ,ĝ k ).

(B.19)
The E and M steps are iterated until convergence.
It is also possible to extend our methods to allow for Markovian unobserved heterogeneity, by employing a variant of the classical Baum-Welch algorithm. The computational and statistical details of such a procedure are however more involved and will be described elsewhere.  Consider the following version of the Rust (1987) bus engine replacement problem which is adapted from Arcidiacono and Miller (2011). Each period t = 1, ..., T ; T < ∞, Harold Zurcher decides whether to replace the engine of a bus (a t = 0), or keep it (a t = 1).
Denote his action by j ∈ {0, 1}. Each bus is characterized by a permanent type s ∈ {1, 2}, and the mileage accumulated since the last engine replacement x t ∈ {1, 2, ...}. Harold Zurcher observes both s and x t . The econometrician observes mileage x t , and we make different assumptions on the observability of bus type s. Mileage increases by one unit if the engine is kept in period t and is set to zero if the engine is replaced. The current-period payoff for keeping the engine is given by θ 0 + θ 1 x t + θ 2 s + e 1t , where θ * ≡ {θ 0 , θ 1 , θ 2 } are the structural parameters of interest, and e jt is a choice-specific transitory shock that follows a Type 1 Extreme Value distribution.
The current-period payoff of replacing the engine is normalized to e 0t , and the discount factor is set to 0.9 in this application.
To carry out the simulations, we choose values for the structural parameters (θ 0 = 2, θ 1 = −0.15, θ 2 = 1), and recursively derive the value functions for each possible combination of x, s and t. We then use these to compute the conditional replacement probabilities for the same set of combinations of variables.
We provide results using the linear semi-gradient method to approximate h(a, x) and g(a, x). Our first set of results treats the bus type s as known to the econometrician. We then provide findings for a setting with permanent unobserved heterogeneity. All results are based on 1000 simulations with 1000 buses and T = 30 time periods each. Each round of the simulations begins by generating data for 1000 buses and 2000 time periods.
The mileage of each bus is set to zero in t = 0. We then simulate the choices a t using the conditional replacement probabilities. Finally, we restrict the generated data to 30 time periods between t = 1000 and t = 1030. This is to ensure that our data is close to being drawn from a stationary model. In all simulations, we parameterize h(a, x) and  (2) and (4) report the empirical mean and standard deviations for the estimated parameters. Columns (2)-(3) are based on the estimation method without correction function, columns (4)-(5) report results for the locally robust estimator. The mean squared error are reported in columns (3) and (5), respectively.
g(a, x) using a third order polynomial in s, x t and a t . 12 The choice probabilities η are estimated using a logit model that is a function of the state variables s and x t , where the same polynomial is used as before. 13 Table B.1 shows the results treating bus type s as known, with and without locally robust correction. It can be seen that our estimator produces parameter estimates that are closely centered around the true values. These results are comparable to those found by Arcidiacono and Miller (2011) in a similar version of the bus engine replacement problem. However, in contrast to their CCP method, our estimator does not exploit a finite dependence property. When comparing the results from our locally robust estimator in column (4) to the results from the suboptimal estimator in column (2), it can be seen that the absolute bias is smaller for all three parameter estimates. However the variance of the locally robust estimator is higher due to the sample splitting employed in the locally robust procedure (we used two-fold cross-fitting). Overall, while in theory the locally robust estimator is preferable to the non-robust version, we find that in practice there is very little difference between the two versions of the algorithm.
12 In total, there are k φ = k r = 16 terms. These include a constant, the binary variables s and a t , all x t terms up to a third order, and pairwise and triple interactions between the x t terms and the binary variables. 13 We run our simulations on a MacBook Pro with an M1 chip and 16 GB of RAM. The approximate computation times for one estimation round is 4 seconds without locally robust correction, and 14 seconds with locally robust correction.  (2) reports the empirical mean and standard deviations for the estimated parameters.
The mean squared errors are reported in column (3). The results are based on the estimation method without correction function.
In a final set of simulations, we introduce permanent unobserved heterogeneity into our setting by assuming that the permanent bus type s ∈ {1, 2} is unknown to the researcher.
To generate results for these simulations, we follow the steps outlined in Section B.6 where we pair our techniques with a sequential EM algorithm (Arcidiacono and Jones, 2003), without using a locally robust correction. 14 The results are shown in Table B.2. Once again, our algorithm produces parameter estimates that are closely centered around the true values. Compared to the results without permanent unobserved heterogeneity, the standard deviation of our estimates is slightly higher due to the uncertainty around the bus type s.
14 We run our simulations on a MacBook Pro with an M1 chip and 16 GB of RAM. The approximate computation times is 305 seconds for one estimation round.