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.
104,426 characters · 24 sections · 78 citation commands
Temporal-Difference estimation of dynamic discrete choice models
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.\footnote{See AguirregabiriaMira2010 for a detailed survey of the literature on the estimation of DDC models.} 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 HotzMiller1993's HotzMiller1993 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., kalouptsidi2014time and almagro2019location 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, barwick2015costs and kalouptsidi2018detection 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.\footnote{Yet another alternative is to use forward Monte Carlo simulations (BBL2007, HMSS1994), but this again becomes very involved as the number of continuous state variables or players increases. The use of a finite number of Monte Carlo simulations also adds bias to the estimates.}
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, AguirregabiriaMira2002) 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 kalouptsidi2018detection 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., AguirregabiriaMira2007; kalouptsidi2014time) or VAR (barwick2015costs) 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, AguirregabiriaMira2002 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 CEINR2018. 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 (BBL2007; AguirregabiriaMira2007; PesendorferSchmidt-Dengler2008) 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 DLR1977.
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 AgMag2018. 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 AguirregabiriaMira2007. We use the linear semi-gradient 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.
The paper contributes to the literature on estimation of DDC models. Rust1987 is the seminal work in this literature. Motivated by computational considerations, HotzMiller1993 HotzMiller1993 propose the CCP algorithm. The CCP idea has subsequently been refined by HMSS1994 who suggest a simulation-based method, and AguirregabiriaMira2002 who develop a pseudo-likelihood estimator. ArcidiaconoMiller2011 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 AguirregabiriaMira2002,AguirregabiriaMira2007 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 KasaharaShimotsu2008. 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.
ackerberg2014asymptotic and CEINR2018 consider semi-parametric estimation using ML methods when either finite dependence or a “terminal action” property holds (HotzMiller1993). CEINR2018 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, semenova2018machine 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 lange2012batch for a survey.\footnote{See SuttonBarto2018 for a detailed treatment of RL in genreal.} It is closely related to the idea of `experience replay', an important ingredient of RL algorithms that achieve human level play in Atari games (mnih2015human). A key step in RL, including batch RL, is the estimation of value functions. TD learning methods, first formulated by sutton1988learning, 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 TsitsiklisRoy1997 for linear semi-gradients, and munos2008finite for AVI. While TsitsiklisRoy1997 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 munos2008finite, 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 (ref)).
TD methods are distinct from other value function approximation methods developed in economics, e.g., parametric policy iteration (benitez2000comparison), simulation and interpolation (keane1994solution), and sieve value function iteration (arcidiacono2013approximating). 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 arcidiacono2013approximating; 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 (ref) outlines the setup of the DDC model and fixes notation. Section (ref) describes our TD estimation method for the functional approximations of the value functions. Section (ref) proves its theoretical properties and describes the second-step estimation of the structural parameters under discrete and continuous state variables. Section (ref) discusses the estimation of dynamic discrete games. Section (ref) provides Monte Carlo simulations for our algorithm. Section (ref) concludes. The Online Appendix discusses extensions to permanent unobserved heterogeneity.
We start with a single agent DDC model. In particular, we consider an infinite horizon model in discrete time with $i=1,\dots,n$ agents. We assume that the individuals are homogeneous, relegating extensions for unobserved heterogeneity to Online Appendix (ref). 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)^{\intercal}\theta+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}(\cdot)$. Let $(a^{\prime},x^{\prime})$ denote the one-period ahead random variables immediately following the actions and states $(a,x)$, where $x^{\prime}\sim K(\cdot\vert a,x)$, with $K(\cdot\vert 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(\cdot|a,x)$. The utility from future periods is discounted by $\beta$.
Agent $i$ chooses actions ${\bf a}_{i}=(a_{i1},a_{i2},\dots)$ to sequentially maximize the discounted sum of payoffs \[ E\left[\sum_{t=1}^{\infty}\beta^{t}\left\{ z(x_{it},a_{it})^{\intercal}\theta^{*}+e_{it}\right\} \right]. \] The econometrician observes a panel consisting of state-action pairs for all individuals, $({\bf x}_{i},{\bf a}_{i})=\{(x_{i1},a_{i1}),\dots,(x_{iT},a_{iT})\}$, for $T$ periods (note, however, that the agent maximizes an infinite horizon objective, not a fixed $T$ one). Typically $T\ll n$ in applications, so we work within an asymptotic regime where $n\to\infty$ but $T$ is fixed. Using this data, the econometrician aims to recover the structural parameters $\theta^{*}$.
In this paper, we study the CCP approach for estimating $\theta^{*}$ (HotzMiller1993). CCP methods are based on the conditional choice probabilities of choosing action $a$ given state $x$. We denote these by $P_{t}(a\vert 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\vert x)$ is independent of $t$). Denote $e(a,x)$ as expected value of the idiosyncratic error term $e$ given that action $a$ was chosen. HotzMiller1993 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\vert x)$, i.e $e(a,x)=\mathcal{G}(P(a\vert 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)=\gamma-\ln P(a\vert x)$, where $\gamma$ 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 pseudo-likelihood function \[ Q(\theta)=\sum_{i=1}^{n}\sum_{t=1}^{T-1}\log\frac{\exp\left\{ h(a_{it},x_{it})^{\intercal}\theta+g(a_{it},x_{it})\right\} }{\sum_{a}\exp\left\{ h(a,x_{it})^{\intercal}\theta+g(a,x_{it})\right\} }, \] where for any $(a,x)$, $h(.)$ and $g(.)$ solve the following recursive expressions:
Here, $\mathbb{E}[.]$ denotes the expectation over the distribution of $(a^{\prime},x^{\prime})$ conditional on $(a,x)$. Note that $\mathbb{E}\left[h(a^{\prime},x^{\prime})\vert a,x\right]$ is a function of $K(\cdot|a,x),P(.\vert 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(\cdot|\cdot)$ and $P(.\vert.)$. Since the latter are unknown, the current literature proceeds by first estimating these as $(\hat{K},\hat{P})$. Typically, $\hat{K}$ is obtained by MLE based on a parametric form of $K(x^{\prime}\vert a,x;\theta_{f})$, while $\hat{P}$ is estimated non-parametrically using either a blocking scheme or kernel regression. Then, given $(\hat{K},\hat{P})$, the values of $h(.)$ and $g(.)$ can be estimated by solving the recursive equations ((ref)). In the next section, we propose an alternative algorithm for maximizing $Q(\theta)$ that directly estimates $h(\cdot)$ and $g(\cdot)$ in a single step without requiring any knowledge about or estimation of $K(\cdot\vert\cdot)$.
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 (ref), but we do require it for deriving a locally robust estimator. Let $\mathbb{P}$ denote the stationary population (i.e, in the limit as $n\to\infty$) distribution of $(a,x,a^{\prime},x^{\prime})$, and $\mathbb{E}[\cdot]$ the corresponding expectation over $\mathbb{P}$. Define $\mathbb{E}_{n}[\cdot]$ as the expectation over the empirical distribution, $\mathbb{P}_{n}$, of $(a,x,a^{\prime},x^{\prime})$. In particular, $\mathbb{E}_{n}[f(a,x,a^{\prime},x^{\prime})]:=(n(T-1))^{-1}\sum_{i=1}^{n}\sum_{t=1}^{T-1}f(a_{it},x_{it},a_{it+1},x_{it+1})$, i.e we always drop the last time period in the summation index even if $f(\cdot)$ does not depend on $a^{\prime},x^{\prime}$.
Let $\mathcal{H}$ denote the space of all square integrable functions over the domain $\mathcal{A}\times\mathcal{X}$ of $(a,x)$. Define the pseudo-norm $\left\Vert \cdot\right\Vert _{2}$ over $\mathcal{\mathcal{H}}$ as $\left\Vert f\right\Vert _{2}:=\mathbb{E}[\vert f(a,x)\vert^{2}]^{1/2}$ for all $f\in\mathcal{\mathcal{H}}$. We use $\left|\cdot\right|$ to denote the usual Euclidean norm on a Euclidean space.
This section presents our TD method for estimating $h(.)$ and $g(.)$. Note that $h(\cdot)$ is a vector, of the same dimension as $\theta$. 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 $\theta^{*}$, 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 \[ \delta_{z}(a,x,a^{\prime},x^{\prime};f):=z(a,x)+\beta f(a^{\prime},x^{\prime})-f(a,x), \] and the dynamic programming operator by \[ \Gamma_{z}[f](a,x):=z(a,x)+\beta\mathbb{E}[f(a^{\prime},x^{\prime})\vert a,x]. \] Clearly, $h(a,x)$ is the unique fixed point of $\Gamma_{z}[\cdot]$. TD estimation involves approximating $h(a,x)$ using a functional class $\mathcal{F}$, where each element $h(\cdot;\omega)$ of $\mathcal{F}$ is indexed by a finite-dimensional vector $\omega$. The aim in TD estimation is to ostensibly minimize the mean-squared TD error \[ \textrm{TDE}(\omega):=\mathbb{E}\left[\left\Vert z(a,x)+\beta h(a^{\prime},x^{\prime};\omega)-h(a,x;\omega)\right\Vert ^{2}\right]. \] However, this minimization problem is neither computationally feasible nor is it proven to converge when $h\notin\mathcal{F}$. Instead, two approaches are commonly used.
The first approach, the semi-gradient method, involves updating $\omega$ using the semi-gradients
for some small value of $\gamma$. As the name suggests, the above is not a complete gradient as the derivative does not take into account how $\omega$ affects the `target', i.e., the future value $h(a^{\prime},x^{\prime};\omega)$. Nevertheless, for linear functional classes $\mathcal{F}$, it is possible to explicitly characterize the limit point of the updates, $\omega^{*}$, and compute it directly. Section (ref) describes this in greater detail. In the RL literature, it is popular to employ semi-gradients with neural networks as the functional class $\mathcal{F}$, but it appears difficult to extend our theoretical analysis to this setting.
The second approach, Approximate Value Iteration (AVI; munos2008finite), employs the idea of `target networks'. In this approach, the parameters in the future value of $h$, i.e $h(a^{\prime},x^{\prime};\omega)$, are fixed at the current $\omega$, 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 ((ref)) using gradient descent, the updates (within each iteration) would look similar to ((ref)) except for fixing the value of $\omega$ in $h(a^{\prime},x^{\prime};\omega)$ at the past values. After the gradient updates converge, i.e at the end of the iteration, $h(a^{\prime},x^{\prime};\omega)$ is revised with the new $\omega$. The semi-gradient approach can thus be considered a one-step variant of AVI. Section (ref) describes AVI in greater detail. We characterize the theoretical properties of AVI under general functional classes $\mathcal{F}$ including neural networks, random forests, LASSO etc.
The approximation to $g$ follows similarly after replacing $\delta_{z}(\cdot;f),\Gamma_{z}[\cdot]$ by
The semi-gradient approach with linear $\mathcal{F}$ is particularly well suited for computation. Let $\phi(a,x)$ consist of a set of basis functions over the domain $(a,x)$. Then the linear approximation class is $\mathcal{F}\equiv\{\phi(a,x)^{\intercal}\omega:\omega\in\mathbb{R}^{k}\}$. Denote the projection operator onto $\mathcal{F}$ by $P_{\phi}$: \[ P_{\phi}[f](a,x):=\phi(a,x)^{\intercal}\mathbb{E}[\phi(a,x)\phi(a,x)^{\intercal}]^{-1}\mathbb{E}[\phi(a,x)f(a,x)]. \]
For linear basis functions, it can be shown, e.g TsitsiklisRoy1997, that the sequence of functional approximations $h(a,x;\omega_{j}):=\phi(a,x)^{\intercal}\omega_{j}$ converges to $h^{*}:=\phi(a,x)^{\intercal}\omega^{*}$, defined as the fixed point of the projected dynamic programming operator $P_{\phi}\Gamma_{z}[\cdot]$ (i.e., $P_{\phi}\Gamma_{z}[h^{*}]=h^{*}$). Based on this characterization, we show in Lemma (ref) in Appendix (ref) that $h^{*}(a,x)=\phi(a,x)^{\intercal}\omega^{*}$, where
Lemma (ref) in Appendix (ref) assures that $\mathbb{E}\left[\phi(a,x)\left(\phi(a,x)-\beta\phi(a^{\prime},x^{\prime})\right)^{\intercal}\right]$ is indeed non-singular as long as $\beta<1$ and $\mathbb{E}[\phi(a,x)\phi(a,x)^{\intercal}]$ is non-singular. As defined above, $\omega^{*}$ cannot be computed directly, since it is a function of the true expectation $\mathbb{E}[\cdot]$. We can however obtain an estimator, $\hat{\omega}$, after replacing $\mathbb{E}[\cdot]$ with the sample expectation $\mathbb{E}_{n}[\cdot]$:
Using $\hat{\omega}$, we obtain an estimate of $h(\cdot)$ as $\hat{h}(a,x)=\phi(a,x)^{\intercal}\hat{\omega}$.
We now turn to the estimation of $g(\cdot)$. As with $h(\cdot)$, we approximate $g(\cdot)$ using basis functions $r(a,x)$, which may generally be different from $\phi(a,x)$. Let $P_{r}$ denote the projection operator onto the space $\{r(a,x)^{\intercal}\xi:\xi\in\mathbb{R}^{k}\}$. The limit of the semi-gradient iterations is $g^{*}(a,x):=r(a,x)^{\intercal}\xi^{*},$ defined as the fixed point of $P_{r}\Gamma_{e}[\cdot]$. Assuming $e(a,x)$ is known, we obtain the following characterization of $\xi^{*}$ in analogy with ((ref)):
In the above, $e(a,x)=\gamma-\ln P(a\vert x)$ is a function of choice probabilities, which are unknown. Denote $\eta(a,x):=P(a\vert x)$. Suppose that we have access to a non-parametric estimator $\hat{\eta}$ of $\eta$. This can be obtained, e.g., through series or kernel regression. We can then plug in this estimate to obtain $e(a,x;\hat{\eta}):=\gamma-\ln\hat{\eta}(a,x)$. This in turn enables us to estimate $\xi^{*}$ using $\hat{\xi}$, computed as
Using the above, we obtain an estimate of $g(\cdot)$ as $\hat{g}(a,x)=r(a,x)^{\intercal}\hat{\xi}$.
Interestingly, estimation of $\xi^{*}$ is unaffected to a first order by the estimation of $\hat{\eta}$, even though the latter converges to the true $\eta$ at non-parametric rates (see Section (ref) for a formal statement). This is because an orthogonality property holds for the estimation of $\xi$, in that
where $\partial_{\eta}\cdot$ denotes the Fr\'echet derivative with respect to $\eta$. To show ((ref)), expand
where the first equality follows from the Markov property. Consider the functional $M(\tilde{\eta}):=\mathbb{E}\left[\left.\ln\tilde{\eta}(a^{\prime},x^{\prime})\right|x^{\prime}\right]$ at different candidate values $\tilde{\eta}(\cdot,\cdot)$. At the true conditional choice probability, $\eta$, $M(\tilde{\eta})$ becomes the conditional entropy of $P(a\vert x^{\prime}$) and attains its maximum. Hence, $\partial_{\eta}\mathbb{E}\left[\left.\ln\eta(a^{\prime},x^{\prime})\right|x^{\prime}\right]=0$ and ((ref)) follows from ((ref)). Consequently, $\hat{\xi}$ is a locally robust estimator for $\xi$.
Note that computation of $\hat{\omega}$ and $\hat{\xi}$ only involves solving linear equations of dimension $\textrm{dim}(\phi)$ and $\textrm{dim}(r)$, respectively. This is computationally very cheap. Using $\hat{h}(a,x)$ and $\hat{g}(a,x)$, we can in turn estimate $\theta^{*}$ in many different ways. For instance, we can plug them into the PMLE estimator
However, such plug-in estimates are sub-optimal. In Section (ref), we suggest a locally robust version of ((ref)).
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$. Then, we show in Online Appendix (ref) that the resulting estimate of $h(a,x)$ is exactly the same as that obtained from the standard CCP estimators, if both the choice and transition probabilities were estimated using cell values.
For a feasible estimation procedure using AVI, we can replace $\mathbb{E}[\cdot]$ by $\mathbb{E}_{n}[\cdot]$ in ((ref)). The procedure builds a sequence of approximations $\{\hat{h}_{j}:=h(a,x;\hat{\omega}_{j}),j=1,\dots,J\}$ for $h$, where
The process can be started with an arbitrary initialization, e.g., $\hat{h}_{1}(a,x)=z(a,x)$. The maximum number of iterations, $J$, is only limited by computational feasibility.
The minimization problem in ((ref)) is equivalent to a prediction problem using the functional class $\mathcal{F}$, where the outcomes are $z(a,x)+\beta\hat{h}_{j}(a^{\prime},x^{\prime})$ evaluated at the various sample draws of $(a,x,a^{\prime},x^{\prime})$. Hence, the estimation target for $\hat{h}_{j+1}$ is the conditional expectation $\mathbb{E}[z(a,x)+\beta\hat{h}_{j}(a^{\prime},x^{\prime})\vert a,x]\equiv\Gamma_{z}[\hat{h}_{j}](a,x)$. In other words, each $\hat{h}_{j+1}$ is a non-parametric approximation to $\Gamma_{z}[\hat{h}_{j}]$, and in this manner AVI builds a series of approximate value function iterations (as its name suggests).
The interpretation of ((ref)) 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(\cdot)$ is similar: we construct a sequence of approximations $\{\hat{g}_{j},j=1,\dots,J\}$ for $g$ as \[ \hat{g}_{j+1}=\operatorname*{arg\,min}_{f\in\mathcal{F}}\mathbb{E}_{n}\left[\left\Vert \beta e(a^{\prime},x^{\prime};\hat{\eta})+\beta\hat{g}_{j}(a^{\prime},x^{\prime})-f(a,x)\right\Vert ^{2}\right]. \] As in Section (ref), it will be shown that the estimation error of $\eta$ is first-order ignorable for the estimation of $g$. Using $\hat{h}(a,x)$ and $\hat{g}(a,x)$, we can, as before, estimate $\theta^{*}$ in many different ways, including the PMLE estimator (ref).
Compared to the semi-gradient approach, AVI is computationally more expensive as it requires solving $J$ prediction problems (in Section (ref) we show that in the worst case $J\approx\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.
Both the semi-gradient and AVI methods will require choosing tuning parameters. For AVI this is straightforward: as each iteration is a non-parametric 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_{\phi}=\textrm{dim}(\phi)$ and $k_{r}=\textrm{dim}(r)$ of the basis functions. In analogy with AVI, we propose selecting both through a procedure akin to cross-validation. The value of $\omega$ 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 mean-squared TD error $\mathbb{E}_{n,\textrm{test}}[\delta_{h}^{2}(a,x,a^{\prime},x^{\prime};\hat{h})]$ on the test dataset. The values of $k_{\phi},k_{r}$ are chosen to minimize the mean squared TD error.
We take the utility function to be linear in $\theta$ 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;\theta)$. Then, $\theta^{*}$ is the maximizer of the pseudo-likelihood criterion \[ Q(\theta)=\sum_{i=1}^{n}\sum_{t=1}^{T-1}\log\frac{\exp\left\{ h(a_{it},x_{it};\theta)+g(a_{it},x_{it})\right\} }{\sum_{a}\exp\left\{ h(a,x_{it};\theta)+g(a,x_{it})\right\} }, \] where, for each $\theta$, $h(.;\theta)$ is defined analogously to ((ref)) with $z(a,x;\theta)$ replacing $z(a,x)$.
We can use our TD procedures to obtain a functional approximation $\hat{h}(a,x;\theta)$ for $h(a,x;\theta)$ at each $\theta$ (estimation of $g(\cdot)$ remains unchanged as it does not depend on $\theta$). As argued earlier, computation of $\hat{h}(\cdot;\theta)$ can be very fast. Appealingly, the matrix $\mathbb{E}_{n}\left[\phi(a,x)\left(\phi(a,x)-\beta\phi(a^{\prime},x^{\prime})\right)^{\intercal}\right]$ employed in the linear semi-gradient estimate ((ref)) does not feature $z(a,x;\theta)$, and therefore only has to be inverted once. We can then plug in the values of $\hat{h}(.;\theta)$ and $\hat{g}(\cdot)$ to estimate $\theta^{*}$ as \[ \hat{\theta}=\arg\max_{\theta\in\Theta}\hat{Q}(\theta);\quad\hat{Q}(\theta):=\sum_{i=1}^{n}\sum_{t=1}^{T-1}\log\frac{\exp\left\{ \hat{h}(a_{it},x_{it};\theta)+\hat{g}(a_{it},x_{it})\right\} }{\sum_{a}\exp\left\{ \hat{h}(a,x_{it};\theta)+\hat{g}(a,x_{it})\right\} }. \] A locally robust counterpart to $\hat{Q}(\theta)$ can be derived in the same manner as in Section (ref). However, computation of $\hat{\theta}$ is more involved as it is no longer convex optimization.
In Online Appendix (ref) we incorporate permanent unobserved heterogeneity into our models by pairing our TD methods with the sequential Expectation-Maximization (EM) algorithm (arcidiacono2003finite). The resulting algorithm can handle discrete heterogeneity in both individual utilities and transition densities. Online Appendix (ref) provides Monte-Carlo evidence suggesting that the algorithm works well in practice.
We characterize rates of convergence for estimation of $h(\cdot)$ and $g(\cdot)$ under both semi-gradients and AVI.
We impose the following assumptions for estimation of $h(\cdot)$.
\begin{asm1} (i) The basis vector $\phi(a,x)$ is linearly independent (i.e. $\phi(a,x)^{\intercal}\omega\linebreak=0$ for all $(a,x)$ if and only if $\omega=0$). Additionally, the eigenvalues of $\mathbb{E}[\phi(a,x)\linebreak\phi(a,x)^{\intercal}]$ are uniformly bounded away from zero for all $k_{\phi}$.
(ii) $\vert\phi(a,x)\vert_{\infty}\le M$ for some $M<\infty$.
(iii) There exists $C<\infty$ and $\alpha>0$ such that $\left\Vert h-P_{\phi}[h]\right\Vert _{2}\le Ck_{\phi}^{-\alpha}$.
(iv) The domain of $(a,x)$ is a compact set, and $\vert z(a,x)\vert_{\infty}\le L$ for some $L<\infty$.
(v) $k_{\phi}\to\infty$ and $k_{\phi}^{2}/n\to0$ as $n\to\infty$.
\end{asm1}
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 $\alpha$ is related to the smoothness of $h(\cdot)$. newey1997convergence shows that for splines and power series, $\alpha=r/d$, where $r$ is the number of continuous derivatives of $h(a,\cdot)$ 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 properties of $h(a,\cdot)$ are discussed in Online Appendix (ref), where we provide primitive conditions on $z(a,x),K(x^{\prime}\vert a,x)$ that ensure existence of $r$ continuous derivatives of $h(a,\cdot)$ for each $a\in\mathcal{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_{\phi}$ is not relevant up to a first order since we propose estimators of $\theta^{*}$ that are locally robust to estimation of $h(\cdot)$.
We then have the following theorem on the estimation of $h(a,x)$:
We prove Theorem (ref) in Appendix (ref) by adapting the results of TsitsiklisRoy1997. The first part of Theorem (ref) assures that both population and empirical TD fixed points exist. The second and third parts of Theorem (ref) imply that the approximation bias and MSE of linear semi-gradients are analogous to those of standard series estimation apart from a $(1-\beta)^{-1}$ factor.
For the estimation of $\hat{\xi}$ we make use of cross-fitting as a technical device to obtain easy-to-verify assumptions for the estimation of $\hat{\eta}$. This entails the following: we randomly partition the data into two folds. We estimate $\hat{\xi}$ separately for each fold using $\hat{\eta}$ estimated from the opposite fold. The final estimate of $\xi^{*}$ is the weighted average of $\hat{\xi}$ 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)$.
\begin{asm2} (i) The basis vector $r(a,x)$ is linearly independent, and the eigenvalues of $\mathbb{E}[r(a,x)r(a,x)^{\intercal}]$ are uniformly bounded away from zero for all $k_{r}$.
(ii) $\vert r(a,x)\vert_{\infty}\le M$ for some $M<\infty$.
(iii) There exists $C<\infty$ and $\alpha>0$ such that $\left\Vert g(a,x)-P_{r}[g(a,x)]\right\Vert _{2}\le Ck_{r}^{-\alpha}$.
(iv) The domain of $(a,x)$ is a compact set, and $\vert e(a,x)\vert_{\infty}\le L<\infty$.
(v) $k_{r}\to\infty$ and $k_{r}^{2}/n\to0$ as $n\to\infty$.
(vi) $\hat{\xi}$ is estimated from a cross-fitting procedure described above. The conditional choice probability function satisfies $\eta(a,x)\ge\delta>0$, where $\delta$ is independent of $a,x$. Additionally, $\left\Vert \eta-\hat{\eta}\right\Vert _{\infty}=o_{p}(1)$ and $\left\Vert \eta-\hat{\eta}\right\Vert _{2}^{2}=o_{p}(n^{-1/2})$.
\end{asm2}
Assumption 2 is a direct analogue of Assumption 1, except for the last part which provides regularity conditions when $\eta(\cdot)$ is estimated. These conditions are typical for locally robust estimates and only require the non-parametric function $\eta(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 (ref), which we prove in Appendix (ref) .
We can expand the estimation error $\left\Vert \hat{h}_{J}-h\right\Vert _{2}$ in terms of the non-parametric estimation errors $\left\Vert \Gamma_{z}[\hat{h}_{j-1}]-h_{j}\right\Vert _{2}$ for $j=1,\dots,J$. In particular, since $\Gamma_{z}[h]=h$ and $\Gamma_{z}[\cdot]$ is a $\beta$-contraction, we have
Iterating the above gives
Equation ((ref)) can be considered a special case of error propagation (munos2008finite).
Recall that $\hat{h}_{1}$ is an arbitrary initialization. It is thus straightforward to provide conditions under which $\left\Vert h-\hat{h}_{1}\right\Vert _{2}$ is bounded by some constant $M_{1}$. As for the second term in ((ref)), recall from the discussion in Section (ref) that the minimization problem ((ref)) corresponds to non-parametric estimation of $\Gamma_{z}[\hat{h}_{j-1}]$ using the functional class $\mathcal{F}$. Most machine learning methods come with guarantees on the non-parametric estimation rate $\left\Vert \Gamma_{z}[\hat{h}_{j-1}]-\hat{h}_{j}\right\Vert _{2}$.
We now describe our assumptions for AVI. Let $\mathcal{X}$ denote the $d$-dimensional space of $x$, and define $\mathcal{W}^{\gamma,\infty}(\mathcal{X})$ as the H\"{o}lder ball with smoothness parameter $\gamma$: \[ \mathcal{W}^{\gamma,\infty}(\mathcal{X}):=\left\{ f:\max_{0<\vert p\vert\le\gamma}\sup_{x\in\mathcal{X}}\left|D^{p}f\right|\le M\right\} . \]
\begin{asm3} (i) The domain, $\mathcal{X}$, of $x$ is compact and there exist $M_{0},M<\infty$ such that $\left|h\right|_{\infty}\le M_{0}$ and $h(\cdot,a)\in\mathcal{W}^{\gamma,\infty}(\mathcal{X})$ for each $a$.
(ii) $\left|\hat{h}_{1}\right|_{\infty}\le M_{0}$ and $\left\Vert h-\hat{h}_{1}\right\Vert _{2}\le M_{1}$ for some $M_{1}<\infty$.
(iii) $\left|\Gamma_{z}[f]\right|_{\infty}\le M_{0}$ and $\Gamma_{z}[f](\cdot,a)\in\mathcal{W}^{\gamma,\infty}(\mathcal{X})$ for all $a\in\mathcal{A}$ and $\left\{ f:\left|f\right|_{\infty}\le M_{0}\right\} $.
(iv) The candidate class of functions $\mathcal{F}$ is such that $\left|f\right|_{\infty}\le M_{0}$ for all $f\in\mathcal{F}.$ Additionally, consider the non-parametric estimation problem\\ $\hat{f}=\operatorname*{arg\,min}_{\tilde{f}\in\mathcal{F}}n^{-1}\sum_{i=1}^{n}(y_{i}-\tilde{f}(a_{i},x_{i}))^{2}$, where $y_{i}$ is compactly supported and $\mathbb{E}[y_{i}\vert a_{i},x_{i}]=f(a_{i},x_{i})$ for some $f\in\mathcal{W}^{\gamma,\infty}(\mathcal{X})$. Then, uniformly over all $f\in\mathcal{W}^{\gamma,\infty}(\mathcal{X})$, $\mathbb{E}\left[\left\Vert f-\hat{f}\right\Vert _{2}\right]\le Cn^{-c}$ for constants $C<\infty$, $c>0$ independent of $n$, but which may depend on $M,M_{0},\gamma$.
\end{asm3}
Assumption 3(i) is a standard requirement in non-parametric estimation. The assumption of $\gamma$-H\"{o}lder continuity is taken from farrell2021deep. Assumption 3(ii) is a mild condition on the initialization $\hat{h}_{1}$. Assumption 3(iii), which is novel to this paper, is a crucial smoothness condition requiring the operator $\Gamma_{z}[\cdot](\cdot,a)$ to map all bounded $f$ onto $\mathcal{W}^{\gamma,\infty}(\mathcal{X})$. In Online Appendix (ref), we show that this is satisfied if $z(a,\cdot)$ and $K(x^{\prime}\vert a,\cdot)$ are $\gamma$-H\"{o}lder continuous.
Assumption 3(iv) is a high-level condition on the machine learning (ML) method $\mathcal{F}$. The requirement of bounded $f$ implies that the ML method cannot diverge in the $l_{\infty}$ sense, see farrell2021deep 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 $\mathcal{W}^{\gamma,\infty}(\mathcal{X})$ at the rate of at least $n^{-c}$. Most ML methods are proven to satisfy this. Consider, for instance, the class $\mathcal{F}$ of MLPs of width $W$ and depth $L$; MLPs and, more generally, Neural Networks are widely used in RL. The results of farrell2021deep imply that for $W\asymp n^{\frac{d}{2(\gamma+d)}}\ln^{2}n$ and $L\asymp\ln n$, \[ \sup_{f\in\mathcal{W}^{\gamma,\infty}(\mathcal{X})}\mathbb{E}\left[\left\Vert f-\hat{f}\right\Vert _{2}\right]\le C\left\{ n^{-\frac{\gamma}{2(\gamma+d)}}\ln^{4}n+\sqrt{\frac{\ln\ln n}{n}}\right\} . \] Thus, Assumption 3(iv) is satisfied for MLPs. See biau2012analysis for related results on random forests.
Assumptions 3(iii) and 3(iv) imply that one can estimate $\Gamma[f]$ for any $\left|f\right|_{\infty}\le M_{0}$ at the $n^{-c}$ rate, i.e., $\sup_{j}\mathbb{E}\left[\left\Vert \Gamma_{z}[\hat{h}_{j-1}]-\hat{h}_{j}\right\Vert _{2}\right]\le Cn^{-c}.$ Combined with ((ref)), this proves:
The first term in the expression for $\mathbb{E}\left[\left\Vert h-\hat{h}_{J}\right\Vert _{2}\right]$ from Theorem (ref) 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\asymp\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, $\hat{h}_{1}$ can be the linear semi-gradient estimator. This is fast to compute and ensures $M_{1}=o_{p}(1)$. Incidentally, Theorem (ref) 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 $\hat{g}$, we again assume cross-fitting is employed as in Theorem (ref), i.e., $\hat{\eta}$ is computed from one half of the data, and $\hat{g}$ is computed using AVI on the other half, taking $\hat{\eta}$ as given.
\begin{asm4} (i) Assumptions 3(i) - 3(iv) hold after replacing $(h,\hat{h}_{1},\Gamma_{z}[\cdot])$ with $(g,\hat{g}_{1},\Gamma_{e}[\cdot])$.
(ii) $\hat{g}$ is estimated from a cross-fitting procedure. The conditional choice probability function satisfies $\eta(a,\cdot)\in\mathcal{W}^{\gamma,\infty}(\mathcal{X})$ for all $a$, and $\eta(a,x)\ge\delta>0$, where $\delta$ is independent of $a,x$. Additionally, $\left\Vert \eta-\hat{\eta}\right\Vert _{2}^{2}=o_{p}(n^{-1/2})$ and with probability approaching one, $\hat{\eta}(a,\cdot)\in\mathcal{W}^{\gamma,\infty}(\mathcal{X})$ for all $a$. \end{asm4}
Assumption 4(ii) is similar to Assumption 2(vi) with the additional requirement that $\eta$ and $\hat{\eta}$ be $\alpha$-H\"{o}lder continuous, the latter with probability approaching one. Note that $\hat{\eta}$ would be $\alpha$-H\"{o}lder continuous if the non-parametric estimator consistently estimates not only $\eta$ but also its first $\alpha$ derivatives.
Let $\tilde{g}$ denote the fixed point of $\tilde{\Gamma}_{e}[f](a,x):=\beta\mathbb{E}[e(a^{\prime},x^{\prime};\hat{\eta})+f(a^{\prime},x^{\prime})\vert a,x]$. Decompose \[ \left\Vert g-\hat{g}_{J}\right\Vert _{2}\le\left\Vert \tilde{g}-\hat{g}_{J}\right\Vert _{2}+\left\Vert g-\tilde{g}\right\Vert _{2}. \] Since $\hat{\eta}(a,\cdot)\in\mathcal{W}^{\gamma,\infty}(\mathcal{X})$, the first term, $\left\Vert \tilde{g}-\hat{g}_{J}\right\Vert _{2}$, can be bounded by a similar rate as in Theorem (ref) (recall that we take $\hat{\eta}$ as given under cross-fitting). As for the second term, $\left\Vert g-\tilde{g}\right\Vert _{2}$, observe that under Assumption 4,
where the `$\lesssim$' holds with probability approaching one, and uses $\partial_{\eta}\mathbb{E}\left[\left.\ln\eta(a^{\prime},x^{\prime})\right|x^{\prime}\right]=0$, see the discussion following ((ref)). Since $\left\Vert \eta-\hat{\eta}\right\Vert _{2}^{2}=o_{p}(n^{-1/2})$, the above proves:
Estimation of $h(a,x)$ and $g(a,x)$ is inherently non-parametric. This is because $h(a,x)$ and $g(a,x)$ are functions of two non-parametric terms: the choice probabilities $\eta(a,x)$, and the transition densities $K(x^{\prime}\vert a,x)$. The TD estimators implicitly take both into account. Under discrete states, the first-step estimation error does not affect the rates of convergence of structural parameters (see Online Appendix (ref)). With continuous states, however, the first-step non-parametric estimation error does affect the estimation of $\theta^{*}$ to a first order when using the PMLE criterion ((ref)). This is because the estimates for $K(x^{\prime}\vert a,x)$ and $\theta^{*}$ are not orthogonal under a PMLE, which extends to the lack of orthogonality between the estimates for $h(a,x),g(a,x)$ and $\theta^{*}$. Consequently, the PMLE estimator with plug-in values of $\hat{h}$ and $\hat{g}$ will converge at slower than parametric rates.
We can recover $\sqrt{n}$-consistent estimation by adjusting the PMLE criterion to account for the first-stage estimation of $h$ and $g$. Denote $(\tilde{{\bf a}},\tilde{{\bf x}}):=(a,x,a^{\prime},x^{\prime})$ and $m(a,x;\theta,h,g):=\partial_{\theta}\ln\pi(a,x;\theta,h,g)$, where \[ \pi(a,x;\theta,h,g):=\frac{\exp\left\{ h(a,x)\theta+g(a,x)\right\} }{\sum_{\breve{a}}\exp\left\{ h(\breve{a},x)\theta+g(\breve{a},x))\right\} }. \] The PMLE estimator with plug-in estimates solves $\mathbb{E}_{n}[m(a,x;\theta,\hat{h},\hat{g})]=0$, but this is not robust to estimation of $h,g$. Let \[ V(a,x;\theta,h,g):=h(a,x)\theta+g(a,x) \] denote the continuation value given $(a,x)$. Also, define $\lambda(a,x;\theta)$ as the fixed point of the `backward' dynamic programming operator
where $(a^{-\prime},x^{-\prime})$ denotes the past actions and states preceding $(a,x)$, and
In Online Appendix (ref), we show that the locally robust moment corresponding to $m(a,x;\theta,h,g)$ is given by
The construction of the locally robust moment ((ref)) is new. But it is infeasible since $\lambda(\cdot),h(\cdot)$ and $g(\cdot)$ are unknown. However, we can replace these quantities with consistent estimates. We have already described how to estimate $h(\cdot),g(\cdot)$. Let $\tilde{\theta}$ denote the plug-in estimator of $\theta^{*}$ using ((ref)); note that $\tilde{\theta}$ consistently estimates $\theta^{*}$ but is not efficient. An estimator, $\hat{\lambda}(\cdot)$, of $\lambda(\cdot)$ can then be obtained by applying either of our TD estimation methods on ((ref)), with $\tilde{\theta},\hat{h},\hat{g}$ plugged in in place of $\theta,h,g$. For instance, using AVI, we could obtain iterative approximations $\{\hat{\lambda}^{(j)},j=1,\dots,J\}$ for $\lambda(\cdot)$ using \[ \hat{\lambda}^{(j+1)}=\operatorname*{arg\,min}_{f\in\mathcal{F}}\mathbb{E}_{n}\left[\left\Vert -\psi(a,x;\tilde{\theta},\hat{h},\hat{g})+\beta\hat{\lambda}^{(j)}(a^{-\prime},x^{-\prime})-f(a,x)\right\Vert ^{2}\right]. \] Plugging in $\hat{\lambda}(\cdot),\hat{h},\hat{g}$ into ((ref)), we obtain the feasible locally robust moment
Based on the above, we can obtain a locally robust estimator, $\hat{\theta}$, as the solution to $\mathbb{E}_{n}[\zeta_{n}(\tilde{{\bf a}},\tilde{{\bf x}};\theta,\hat{h},\hat{g})]=0$. We recommend obtaining this estimate using cross-fitting, see Section (ref) for more details. Compared to the plug-in estimate ((ref)), our locally robust estimator requires computation of $\lambda(\cdot)$, but even this can be avoided when linear semi-gradients are used for estimating $h,g$ (see Online Appendix (ref)). Solving $\mathbb{E}_{n}[\zeta_{n}(\tilde{{\bf a}},\tilde{{\bf x}};\theta,\hat{h},\hat{g})]=0$ is also computationally easy; the correction term is a constant, and $\partial_{\theta}\zeta_{n}(\tilde{{\bf a}},\tilde{{\bf x}};\theta,\hat{h},\hat{g})=\partial_{\theta}m(a,x;\theta,\hat{h},\hat{g})$ 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.
We focus on the general construction of the locally robust estimator, $\hat{\theta}$, using ((ref)). As mentioned in the previous sub-section, we advocate cross-fitting to obtain this estimator. This entails computing $\tilde{\theta},\hat{\lambda},\hat{h},\hat{g}$ using one half of the sample, say $\mathcal{N}_{2}$, and plugging them into the locally robust moment to compute $\hat{\theta}$ as the solution to $\mathbb{E}_{n}^{(1)}[\zeta_{n}(\tilde{{\bf a}},\tilde{{\bf x}};\theta,\hat{h},\hat{g})]=0$, where $\mathbb{E}_{n}^{(1)}[\cdot]$ is the empirical expectation using only observations from the other half of the sample $\mathcal{N}_{1}$. Following the analysis of CEINR2018, it can be shown that this estimator has the same limiting distribution as the one based on ((ref)). In particular, it achieves parametric rates of convergence. We state the regularity conditions below (for the remainder of this section we allow $\theta^{*}$ to be vector valued):
\begin{asm5} (i) $\theta^{*}\in\Theta$, a compact set, and $\mathbb{E}\left[\zeta(\tilde{{\bf a}},\tilde{{\bf x}};\theta,h,g)\right]=0\iff\theta=\theta^{*}$.
(ii) There exists a neighborhood, $\mathcal{N}$, of $\theta^{*}$ such that uniformly over $\theta\in\mathcal{N}$ and for $\parallel\!\tilde{h}-h\!\parallel$, $\parallel\!\tilde{g}-g\!\parallel$ sufficiently small, $\parallel\negmedspace\partial_{\theta}\zeta(\tilde{{\bf a}},\tilde{{\bf x}};\theta,\tilde{h},\tilde{g})-\partial_{\theta}\zeta(\tilde{{\bf a}},\tilde{{\bf x}};\theta^{*},\tilde{h},\tilde{g})\!\parallel\le d(\tilde{{\bf a}},\tilde{{\bf x}})\left\Vert \theta-\theta^{*}\right\Vert $, where $\mathbb{E}[d(\tilde{{\bf a}},\tilde{{\bf x}})]<\infty$. Furthermore, $G:=\mathbb{E}\left[\partial_{\theta}\zeta(\tilde{{\bf a}},\tilde{{\bf x}};\theta^{*},h,g)\right]$ is invertible.
(iii) $\forall\ \theta\in\mathcal{N}$, $\left\Vert \lambda(\cdot,\cdot;\theta)-\lambda(\cdot,\cdot;\theta^{*})\right\Vert _{2}\le d(a,x)\left\Vert \theta-\theta^{*}\right\Vert $, where $\mathbb{E}[d(a,x)]<\infty$.
(iv) $\left\Vert \hat{h}-h\right\Vert _{2}=o_{p}(n^{-1/4})$ and $\left\Vert \hat{g}-g\right\Vert _{2}=o_{p}(n^{-1/4})$.
(v) $\sup_{\theta}\left\Vert \hat{\lambda}(\cdot,\cdot;\theta)-\lambda(\cdot,\cdot;\theta)\right\Vert _{2}=o_{p}(n^{-1/4})$. \end{asm5}
Assumption 5(i) implies $\theta^{*}$ is identified. Assumption 5(ii) is a mild regularity condition that is similar to Assumption 5 in CEINR2018. Assumption 5(iii) is satisfied if $\psi(\cdot)$ is uniformly Lipschitz continuous in $\theta$. Assumption 5(iv) follows from Theorems (ref)-(ref) under suitable conditions on the degree of smoothness of $h,g$. For instance, it is satisfied for AVI with neural networks if $\gamma\ge d$. Assumption 5(v) requires $\lambda$ to be estimable at faster than $n^{-1/4}$ rates as well. If $h,g$ are known, it is straightforward to derive $n^{-1/4}$ rates as in Theorems (ref)-(ref). For plug-in estimation, we would need additional assumptions. For instance, we could employ three-way sample splitting (the first third of the sample is used to compute $\hat{h},\hat{g}$, which are then plugged into the second third of the sample to estimate $\lambda$), in which case Assumption 5(v) holds under the regularity condition \[ \mathbb{E}\left[\sup_{\theta}\left\Vert \psi(a,x;\theta,\hat{h},\hat{g})-\psi(a,x;\theta,h,g)\right\Vert \right]\le\left\Vert \hat{h}-h\right\Vert _{2}+\left\Vert \hat{g}-g\right\Vert _{2}. \] We refer to chernozhukov2018learning for more details on three-way splitting.
We are now ready to state the main result of this section.
The proof of the above theorem follows by verifying the regularity conditions of CEINR2018. Since these are more or less straightforward to verify given our previous results, we omit the details. For inference on $\hat{\theta}$, the covariance matrix $V$ can be estimated as $\hat{V}=\left(\hat{G}^{\intercal}\hat{\Omega}^{-1}\hat{G}\right)^{-1},$ where
CEINR2018 provide conditions under which $\hat{V}$ is consistent for $V$; these are straightforward to translate to our context but we omit them for brevity. Alternatively, one could employ the bootstrap.
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.
Our setup is based on AguirregabiriaMira2010. We assume a single Markov-Perfect-Equilibrium setup where multiple players $i=1,2,\dots,N$ play against each other in $M$ different markets. Each player chooses among $A$ mutually exclusive actions to maximize an infinite horizon objective. We observe the state of play for $T$ time periods, where both $T$ and the number of players $N$ are assumed fixed while $M\to\infty$. Utility of the players in any time period is affected by the actions of all the others, and a set of states $x$ that are observed by all players. The per-period utility is denoted by $z_{i}(a_{i},a_{-i},x)^{\intercal}\theta^{*}+e_{i}$ for each player $i$, for some finite-dimensional parameter $\theta^{*}$, where $a_{i}$ denotes player $i$'s action, $a_{-i}$ denotes the actions of all other players and $e_{i}$ is an idiosyncratic error term. As in Section (ref), we take $\theta^{*}$ to be scalar to simplify the notation; all our results continue to hold for vector valued $\theta$, as long as each dimension is treated separately. Evolution of the states in the next period is determined by the transition density $K(x^{\prime}\vert a,x)$ where ${\bf a}:=(a_{1},\dots,a_{N})$ denotes the actions of all the players. We denote by $x_{tm}$ the state at market $m$ in time period $t,$ by ${\bf 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}\vert 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):=\gamma-\ln P_{i}(a_{i}\vert x)$.
As in the single agent case, the parameters $\theta^{*}$ 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 ((ref)) in the single agent case, the expectation now averages over the actions of the other players as well.
Previous literature estimates $\theta^{*}$ using a two-step procedure: In the first step, the conditional choice probabilities $P_{i}(a_{i}\vert 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 ((ref)). This step requires integrating over the actions of all the other players. Finally, given the estimated values of $h_{i}(.)$ and $g_{i}(.)$, the parameter $\theta$ is estimated through either pseudo-likelihood (AguirregabiriaMira2007) or minimum distance estimation (PesendorferSchmidt-Dengler2008).
By contrast, our algorithm is a straightforward extension of those suggested in earlier sections for single-agent models. Let $\hat{\eta}_{i}(a_{i},x)$ denote a non-parametric estimate of the choice probabilities for player $i$ and denote $e(a_{i},x;\hat{\eta}_{i})=\gamma-\ln\hat{\eta}_{i}(a_{i},x)$. We apply our TD methods on the recursion ((ref)), separately for each player. The linear semi-gradient estimates are given by $\hat{h}_{i}(a_{i},x)=\phi(a_{i},x)^{\intercal}\hat{\omega}_{i}$ and $\hat{g}_{i}(a_{i},x)=r(a_{i},x)^{\intercal}\hat{\xi}$, where
and for any function $f(\cdot)$, we define
Similarly, the AVI iterations for $h_{i}(\cdot),g_{i}(\cdot)$ are given by
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 ((ref)) nor ((ref)) 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 $\hat{h}_{i},\hat{g}_{i}$ can be plugged into the PMLE criterion ((ref)) to obtain an estimate for $\theta$ as
It is straightforward to construct a locally robust estimator for $\theta$ in analogy with that for single-agent models. We describe this in Online Appendix (ref).
By the same reasoning as in Online Appendix (ref), it possible to show that with discrete states, $h_{i}(.)$ and $g_{i}(.)$ are numerically identical to the estimates obtained by plugging in cell estimates $\hat{P}_{j}(\cdot\vert x)$ and $\hat{K}(.)$ in ((ref)). This implies the psuedo-likelihood with plug-in estimates for $h(.)$ and $g(.)$ is not efficient even with discrete states, as discussed by AguirregabiriaMira2007. 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 PesendorferSchmidt-Dengler2008. With continuous states, one would need to employ locally robust corrections even for the minimum distance estimator to recover parametric rates of convergence for $\theta$. The locally robust correction term can be constructed in a similar way as that for the PMLE criterion.
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 AgMag2018.
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 (ref) to select the preferred set of basis functions. Our simulations are based on the dynamic firm entry game used in AguirregabiriaMira2007.
Finally, we run additional simulations based on the famous Rust1987 bus engine replacement problem (see Online Appendix (ref)). Using this model, we also provide simulation results for a case with permanent unobserved heterogeneity.
Consider the following dynamic firm entry problem described in AgMag2018. 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 $\Pi_{t}=VP_{t}-FC_{t}-EC_{t}+\varepsilon_{t},$ where $VP_{t}$, $FC_{t}$ and $EC_{t}$ denote the firm's variable profit, fixed cost and entry cost, and $\varepsilon_{t}$ is a transitory shock that follows a logistic distribution. Variable profit is given by $VP_{t}=(\theta_{0}^{VP}+\theta_{1}^{VP}z_{1t}+\theta_{2}^{VP}z_{2t})\exp(\omega_{t})$, where $\omega_{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 $FC_{t}=\theta_{0}^{FC}+\theta_{1}^{FC}z_{3t}$, and the entry cost is given by $EC_{t}=(\theta_{0}^{EC}+\theta_{1}^{EC}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 parameters $\theta^{*}\equiv\left\{ \theta_{0}^{VP},\theta_{1}^{VP},\theta_{2}^{VP},\theta_{0}^{FC},\theta_{1}^{FC},\theta_{0}^{EC},\theta_{1}^{EC}\right\} $ are the structural parameters of interest. The exogenous state variables $z_{jt}$ and $\omega_{t}$ are continuous and follow AR(1) processes, where $z_{jt}=\gamma_{0}^{j}+\gamma_{1}^{j}z_{jt-1}+e_{jt}$, and $\omega_{t}=\gamma_{0}^{\omega}+\gamma_{1}^{\omega}\omega_{t-1}+e_{\omega t}$. The error terms $e_{jt},e_{\omega t}$ follow normal $N(0,1)$ distributions. The discount factor $\beta$ is $0.95$.
To carry out the simulations, we choose values for the structural parameters $\theta^{*}$ ($\theta_{0}^{VP}=0.5$, $\theta_{1}^{VP}=1.0$, $\theta_{2}^{VP}=-1.0$, $\theta_{0}^{FC}=1.5$, $\theta_{1}^{FC}=1.0$, $\theta_{0}^{EC}=1.0$, $\theta_{1}^{EC}=1.0$) and for the autoregressive processes of $z_{jt}$ and $\omega_{t}$ ($\gamma_{0}^{j}=0.0$, $\gamma_{1}^{j}=0.6$, $\gamma_{0}^{\omega}=0.2$, $\gamma_{1}^{\omega}=0.6$), and discretize the exogenous state variables to obtain a transition matrix with a 6-point support following Tauchen1986. The resulting dimension of the state space is $2\times6^{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.
Panel A of Table (ref) shows results for 1000 simulations using the linear semi-gradient method. Panel B of Table (ref) presents results of 250 simulations using the AVI method. Each round of the simulations begins by generating new data, where the first-period state variables are drawn from the steady-state distribution. For the results in Panel A, we parameterize $h(a,x)$ and $g(a,x)$ using a second order polynomial in the state variables.\footnote{For the $\omega$'s relating to parameters $\theta_{0}^{VP},\theta_{1}^{VP},\theta_{2}^{VP},\theta_{0}^{FC},\theta_{1}^{FC}$, and for $\xi$, 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(\omega_{t})$ and $z_{2t}\exp(\omega_{t})$ instead of $z_{1t}$ and $z_{2t}$ themselves as state variables. In addition to the terms included above, the $\omega$'s relating to parameters $\theta_{0}^{EC}$ and $\theta_{1}^{EC}$ 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 $\theta_{0}^{EC}$ and $\theta_{1}^{EC}$).} For the results in Panel B, we approximate $h(a,x)$ and $g(a,x)$ using a random forest, where we iterate the AVI procedure 70 times for each round of the simulations. For both the linear semi-gradient and the AVI methods, we estimate the choice probabilities $\eta$ that enter $e(a',x';\eta)$ using a logit model where the explanatory variables are the same as those used as basis functions in Panel A.
We present results generated with and without the locally robust correction. For the results without correction, we obtain estimates for $\theta^{*}$ using ((ref)). To generate the locally robust estimates, we use moment equation ((ref)) for the linear semi-gradient method, and moment equation ((ref)) for the AVI method where we employ a random forest to derive an estimate for the $\lambda(a,x,\tilde{\theta})$ term contained in the locally robust moment. As before, the AVI method for estimation of $\lambda(\cdot)$ iterated 70 times. We also use the sample splitting method described in Section (ref) for the locally robust estimators, and we obtain the final $\hat{\theta}$ as weighted average of the $\theta^{*}$ estimates from the two samples.\footnote{We run our simulations on a MacBook Pro with an M1 chip and 16 GB of RAM. The computation time for one estimation round is about 4 seconds for the linear semi-gradient method without locally robust correction, and 14 seconds with locally robust correction. For the AVI method, the computation time is about 90 seconds without locally robust correction, and 315 seconds with locally robust correction.}
Panel A of Table (ref) 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 (ref) 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 (ref) 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.
We compare our findings from Section (ref) 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 (ref). We then discretize the state space by creating dummy variables for each state variable $z_{1t},z_{2t},z_{3t},z_{4t}$ and $\exp(\omega_{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 (ref).
It can be seen that, compared to the results from Table (ref), the CCP estimator leads to substantially larger bias in some of the estimated parameters. Column (4) shows that the corresponding mean squared errors are large and generally exceed those obtained using our estimators. This is particularly true for parameters $\theta_{1}^{VP}$and $\theta_{2}^{VP}$, 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 (ref) to $0.16$ in Table (ref). These results show that our estimation methods lead to important improvements over existing methods in models with as few as five continuous state variables.
Consider the following firm market entry game, which is similar to that described in AguirregabiriaMira2007. There are $i=1,...,5$ firms (players), and we observe their decision to enter ($a_{itm}=1)$ or not enter ($a_{itm}=0)$ in $m=1,...,M$ different markets for $t=1,...,T$ time periods. Denote a firm's action by $j\in\left\{ 1,0\right\} $. The payoff of each firm $i$ is affected by the decision of all the other firms whether to enter, as well as firm $i$'s previous-period entry decision. Current-period profits when entering are given by \[ \Pi_{itm}=\theta_{RS}\ln(S_{tm})-\theta_{RN}\ln(1+\sum_{j\neq i}a_{jtm})-\theta_{FC}-\theta_{EC}(1-a_{i(t-1)m})+\varepsilon_{itm}, \] where $\ln(S_{tm})$ is a measure of consumer market size of market $m$ in period $t$, and $\varepsilon_{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})=\alpha+\lambda\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 $\beta$ is $0.95$. The parameters $\theta^{*}\equiv\left\{ \theta_{RS},\theta_{RN},\theta_{FC},\theta_{EC}\right\} $ 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}=\left\{ a_{i(t-1)m}:i=1,...,5\right\} $.
To carry out the simulations, we choose values for the structural parameters $\theta^{*}$ ($\theta_{RS}=1,\theta_{RN}=1,\theta_{FC}=1.7,\theta_{EC}=1$), and for the autoregressive process for log market size ($\alpha=1.5$, $\lambda=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 Tauchen1986. As in the Monte Carlo experiments for the firm entry problem in Section (ref), 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 $\nu_{j}(S_{tm},a_{(t-1)m})$ for each of the $2^{5}\times10=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.
We present the results of $1000$ simulations based on the linear semi-gradient method, without employing the locally robust correction. Each round of the simulations begins by generating new data, where the first-period state variables are drawn from the steady-state distribution. In order to assess the sensitivity of our algorithm to different specifications for the basis functions, we parameterize $h(a,x)$ and $g(a,x)$ using different sets of polynomials in the state variables. In particular, we show results where $h(a,x)$ and $g(a,x)$ are approximated using a second, third or fourth order polynomial.\footnote{For the $\omega'$s relating to parameters $\theta_{RS},\theta_{RN},\theta_{FC}$ and for $\xi$, the terms include a constant, terms up to the second/third/fourth order in the state variables $\ln(S_{tm})$ and $\ln(1+\sum_{j\neq 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 $\omega'$s relating to parameter $\theta_{EC}$ also contain the term $(1-a_{i(t-1)m})a_{itm}$. } For all simulations, the choice probabilities $\eta$ that enter $e(a',x';\eta)$ are estimated using individual logit models for each firm, where we use a third order polynomial in the state variables as explanatory variables. We then estimate the parameters $\omega$ and $\xi$ using equation (ref).\footnote{Given the symmetric set-up of the game, we pool the data across players in this application.} Finally, we obtain estimates for the $\theta^{*}$ parameters as the solutions to the pseudo-likelihood function (ref).\footnote{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).}
The results are shown in Table (ref). 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, $\theta_{RN}$. To assess these differences formally, we use the cross-validation procedure described in Section (ref). 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, AguirregabiriaMira2007 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.\footnote{As before, we compute the TD criterion for ten random samples.} While we view this as further evidence that 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.
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 semi-gradient, 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 Rust1987 engine replacement problem show that the proposed algorithms work well in practice.