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.
116,642 characters · 15 sections · 103 citation commands
Efficient and Convergent Sequential Pseudo-Likelihood Estimation of Dynamic Discrete Games
Estimation of dynamic discrete choice models -- particularly dynamic discrete games of incomplete information -- is a topic of considerable interest in economics. Broadly, likelihood-based estimation of these models takes the form
where $Q_{N}$ is the log-likelihood function based on $N$ independent markets, $\theta$ is a finite-dimensional vector of parameters, $Y$ is a vector of important auxiliary parameters, and $G(\theta,Y)=0$ is a vector equality constraint that represents equilibrium conditions. The parameters $\theta$ usually consist of the structural parameters of the model. Common examples of auxiliary parameters $Y$ include expected/integrated value functions or conditional choice probabilities, since the equality constraint is often derived from an equilibrium fixed point condition of the form $G(\theta,Y)\equiv Y-\Gamma(\theta,Y)=0$.
One approach to estimating these models is to directly impose the fixed point equation for each trial value of $\theta$ visited by the optimization algorithm by solving for $Y_{\theta}$ such that $G(\theta,Y_{\theta})=0$. This approach was pioneered by Miller1984, Wolpin1984, Pakes1986, and RustGMC1987 for single-agent models, where the fixed point is unique and can be computed via standard value function iteration or backwards induction. Solution algorithms are available for dynamic games PakesMcguire1994,PakesMcGuire2001, but it is often infeasible to nest those within an estimation routine because the computational burden can be quite large. Furthermore, in games the model may be incomplete due to multiple solutions $Y_{\theta}$ Tamer2003.
These issues with the nested fixed-point approach led researchers to extend conditional choice probability (CCP) estimators -- first introduced in the seminal work of HotzMil1993 for single-agent models -- to the case of dynamic discrete games. Of particular interest here is the nested pseudo-likelihood approach of Aguirregabiria and Mira AgMir2002,AgMir2007.\footnote{Some other examples of CCP estimators are described in HotzMSS1994,BajBenkLevin2007,PakesOstBerry2007,PesSD2008.} They suggest using a $k$-step nested pseudo-likelihood ($k$-NPL) approach, which defines a sequence of estimators, as an algorithm for computing the nested pseudo-likelihood (NPL) estimator, a fixed point of the sequence. In single-agent models, AgMir2002 show that the $k$-NPL estimator is efficient for $k\geq1$ when initialized with a consistent CCP estimate, in the sense that it has the same limiting distribution as the (partial) maximum likelihood estimator. Furthermore, KasShim2008 showed that the sequence converges to the true parameter values with probability approaching one in large samples. Indeed, the Monte Carlo simulations in AgMir2002 show that single-agent $k$-NPL reliably converges to the maximum likelihood estimate. The combination of computational simplicity, efficiency, and convergence stability make $k$-NPL an attractive alternative to other approaches, such as nested fixed point (computationally burdensome) or standard Newton or Fisher scoring steps on the full maximum likelihood problem (often diverges in finite samples in practice).\footnote{This is a well-known limitation of standard Newton steps, and further regularization is often needed to ensure convergence.}
Unfortunately, these attractive properties of $k$-NPL are lost in dynamic games. AgMir2007 show that $k$-NPL estimates are in general not efficient for $k\leq\infty$, although they show that the $\infty$-NPL (taking $k\to\infty$ or until convergence) estimator outperforms the $1$-NPL estimator in efficiency when both are consistent. But PesSD2010 show that the sequence may fail to converge to the equilibrium that generated the data, even with very good starting values, so that $\infty$-NPL may not be consistent (see also, KasShim2012, EgLaiSu2015, and AgMarcoux2019). KasShim2012 show that inconsistency occurs when the NPL mapping is unstable at the data-generating equilibrium, which is essentially equivalent to best-response instability of the equilibrium.\footnote{The NPL mapping is first-order equivalent to the best-response mapping at the equilibrium.}
Another type of CCP estimator is the minimum-distance estimator (AltMill1998,PesSD2008). This estimator is both consistent and efficient in dynamic games, and BugniBunt2020 develop a sequential version of this estimator, which we refer to as $k$-MD. However, Monte Carlo simulations show that these econometric properties come at the expense of greatly increased computational burden, taking 12 to 26 times longer per iteration than $k$-NPL in even the simplest small-scale dynamic games.\footnote{BugniBunt2020 perform Monte Carlo simulations with a very small-scale game (two players, two actions, four states), and even in that setting they remark, “... computing the optimally weighted 1-MD, 2-MD, and 3-MD estimators takes us roughly 33%, 75%, and 80% more time than computing the 20-{[}NPL{]} estimator, respectively.” These translate to roughly 26, 17.5, and 12 times longer per iteration for each respective case.} This difference in computation time is likely to grow with the size of the game, which would be a serious concern for empirical applications. This leaves the researcher with an undesirable tradeoff between the computationally simple $k$-NPL sequence and the more burdensome $k$-MD sequence with better efficiency properties.
But there is another concern shared by the $k$-NPL and $k$-MD sequences: finite sample performance of successive iterations. The $k$-MD estimator uses the NPL mapping to update choice probabilities between iterations, so there is reason to be concerned that it may mimic $k$-NPL's finite sample properties when the data-generating equilibrium is NPL-unstable. While the first-order asymptotic analysis of BugniBunt2020 implies that both sequences are consistent for finite $k$, this asymptotic consistency does not necessarily lead to good performance in finite samples. Indeed, the finite-sample performance of $k$-NPL deteriorates with $k$ when the equilibrium is NPL-unstable. For example, PesSD2008 consider a finite number of iterations for $k$-NPL and find that it is severely biased in large-but-finite samples when the equilibrium is NPL-unstable. In our own Monte Carlo simulations presented later, we find that substantial bias appears rather quickly, even for low values of $k$. This issue may also be a concern for $k$-MD, since it is also consistent for finite $k$ but has unknown stability properties as $k\to\infty$, which may lead to deterioration in finite-sample performance even for fixed $k$, similar to $k$-NPL. This concern arises because KasShim2012 show that instability/inconsistency of the $k$-NPL sequence arises from instability of the NPL mapping used to update the choice probabilities between iterations.\footnote{However, a rigorous econometric analysis of the the $k$-MD estimator's behavior as $k\to\infty$ is beyond the scope of this paper.}
With these concerns about various sequential methods in mind, an important question arises: is there a CCP-based sequence that achieves a balance between computational simplicity, asymptotic efficiency, and good finite-sample properties with any number of iterations -- including as $k\to\infty$ (iterating to convergence) -- in dynamic games? The primary contribution of this paper is to provide such a method, which we name the $k$-step Efficient Pseudo-Likelihood ($k$-EPL) estimator. We show that $k$-EPL estimates are first-order asymptotically equivalent to the maximum likelihood estimate for any number of iterations, $k\geq1$. Thus, every estimate in the sequence is efficient. Furthermore, we also show that higher-order improvements are achieved with iteration, so the $k$-EPL sequence converges to the finite-sample maximum likelihood estimator almost surely. The convergence rate is fast, approaching super-linear as $N\to\infty$. This convergence result for $k$-EPL holds even when the data-generating equilibrium is best-response unstable, rendering $k$-NPL inconsistent.
One key distinction between $k$-EPL and $k$-NPL lies in how we incorporate simultaneous play by multiple agents. While $k$-NPL focuses on single agents' responses to a combination of an exogenous state transition process and other agents' equilibrium play, $k$-EPL incorporates the simultaneous nature of the game and is based on the joint behavior of multiple players. Incorporating this additional information yields increased asymptotic precision.
Despite this conceptual modification, our $k$-EPL estimator retains a simple computational structure similar to $k$-NPL. When utility is linear in the parameters of interest, both $k$-EPL and $k$-NPL iterations proceed in two stages: (i) solve a set of linear systems to generate pseudo-regressors; and (ii) use the pseudo-regressors in a static logit/probit maximum likelihood problem. The linear systems only need to be solved once per iteration in the sequence, and the static logit/probit problem is a low-dimensional, strictly concave problem that has a unique solution and is easy to solve with out-of-the-box optimization software. Because $k$-EPL incorporates all players simultaneously, the linear systems in $k$-EPL have larger dimension than those in k-NPL. While this increases the relative computation time in large-dimensional problems, the increase appears much less severe than the computation time required for $k$-MD. We also find that iterating $k$-EPL to convergence can be faster than doing so with $k$-NPL when the game is not too large, since the need for fewer iterations results in lower overall computation time.
One interesting implication of our higher-order analysis is that iterating $k$-EPL can also provide a convenient algorithm for computing the maximum likelihood estimator. We explore this in some of our Monte Carlo simulations and find that it performs quite well, even with random starting values. However, our primary focus is on the entire $k$-EPL sequence, beginning with consistent initial estimates. We find that much of the practical improvement from iteration is achieved with just a few iterations in our Monte Carlo experiments, suggesting that low values of $k$ can be quite effective when the initial CCP estimates are consistent.
In recent related work, AgMarcoux2019 studied the finite-sample properties of $\infty$-NPL ($k$-NPL with $k\to\infty$) and introduced a variant of the $k$-NPL algorithm that updates the conditional choice probabilities by applying spectral methods to the CCP updates. The goal of their algorithm is to improve convergence properties of $\infty$-NPL for unstable fixed points.\footnote{In the population or in finite samples, the NPL operator may have spectral radius larger than one for some equilibria, rendering it unstable. Conversely, the spectral radius of the EPL operator is zero in the population and near zero in finite samples.} However, there are several differences between their work and ours. First, they limit their analysis of the spectral algorithm to computing the best fixed point of the $k$-NPL sequence, whereas we provide analysis for all the iterations in the $k$-EPL sequence. Second, upon convergence, the $k$-NPL and spectral $k$-NPL algorithms do not produce the maximum likelihood estimator, and convergence can require many iterations. In contrast, our $k$-EPL estimator has the same limiting distribution as the MLE at each iteration and usually converges locally to the MLE after few iterations in finite samples. We verify these properties in our simulation studies.
Our work is also related to methods leveraging Neyman orthogonalization, which has played a central role in recent advances in the broader econometrics literature (see, e.g., ChernozhukovEtAl_TEJ2018,ChernozhukovEtAl_ECMA2022a,ChernozhukovEtAL_ECMA2022b). $k$-EPL leverages a type of quasi-Newton step in its construction, leading to an important “zero Jacobian” property. Consequently, each estimate in the $k$-EPL sequence is asymptotically Neyman orthogonal to the previous estimate, and which leads to many of $k$-EPL's attractive econometric properties.
We demonstrate the application of the $k$-EPL estimator through an empirical analysis of the U.S. wholesale club industry, with a specific focus on its three major players: Sam's Club, Costco, and BJ's. We construct a structural dynamic model to examine the industry's competitive landscape. Leveraging data on club stores operating across the United States from 2009 to 2021, we employ $k$-EPL to estimate structural parameters such as fixed costs, entry costs, the effect of market size, and the competitive effect. Additionally, we consider a counterfactual experiment designed to identify the key determinants of market entry behavior and to explore their potential influence on the industry's structure.
The remainder of the paper proceeds as follows. Section (ref) describes a generic dynamic discrete choice game of incomplete information. Section (ref) describes the $k$-EPL estimator, its asymptotic and finite-sample properties, and its numerical implementation. Section (ref) provides Monte Carlo simulations, and additional simulation results are included in the appendix. Section (ref) describes our empirical application to the U.S. wholesale club industry. Section (ref) concludes. All proofs appear in the Appendix.
Here we describe a canonical stationary dynamic discrete game of complete information in the style of AgMir2007 and PesSD2008. Time is discrete, indexed by $t=1,2,3.\dots$. In a given market, there are $J$ firms indexed by $j\in\mathcal{J}=\{1,2,\dots,|\mathcal{J}|\}$. Given a vector of state variables observable to all agents and the econometrician, $x_{t}$, and its own private information $\varepsilon_{t}^{j}$, each firm chooses an action, $a_{t}^{j}\in\mathcal{A}=\{0,1,2,\dots,|\mathcal{A}|-1\}$. Action zero is the “outside option” when applicable. All players choose their actions simultaneously.
Agents have flow utilities (profits), $U^{j}(x_{t},a_{t}^{j},a_{t}^{-j},\varepsilon_{t}^{j};\theta_{u})$, where $a_{t}^{-j}$ are the actions of the other players. States transition according to $p(x_{t+1},\varepsilon_{t+1}\mid a_{t},x_{t},\varepsilon_{t};\theta_{f})$, and the discount factor is $\beta\in(0,1)$. Agents choose actions to maximize expected discounted utility, \[ \mathrm{E}\left\{ \sum_{s=0}^{\infty}\beta^{s-t}U^{j}(x_{s},a_{s}^{j},a_{s}^{-j},\varepsilon_{s}^{j};\theta_{u})\,\bigg|\,x_{t},\varepsilon_{t}^{j}\right\} . \] The primary parameter of interest is $\theta=(\theta_{u},\theta_{f})$. Furthermore, we impose the following standard assumptions on the primitives.
Assumptions 1--4 here correspond to Assumptions 1--4 in AgMir2007. In most applications in the literature, the private shocks are assumed to be either i.i.d. Type 1 Extreme Value or normal, both of which satisfy the assumptions.
The operative equilibrium concept here will be that of a Markov Perfect Nash equilibrium. We will consider stationary equilibria only, so from here we drop the time subscript. Because moves are simultaneous, the actions of player $j$ do not depend directly on $a^{-j}\in\mathcal{A}^{|\mathcal{J}|-1}$, but rather on $P^{-j}(x)\in\Delta^{|\mathcal{J}|-1}$, where $P^{-j}(x)$ is player $j$'s belief about the other players' probability of playing the corresponding actions in state $x$ and $\Delta$ is the unit simplex in $\mathbb{R}^{\left|\mathcal{A}\right|-1}$. So, from here on out we will work with the following utility function and transition probabilities: \[ u^{j}(x,a^{j};P^{-j},\theta_{u})=\sum_{a^{-j}\in\mathcal{A}^{|\mathcal{J}|-1}}P^{-j}(x,a^{-j})\bar{u}(x,a^{j},a^{-j};\theta_{u}) \]
\[ f^{j}(x'\mid x,a^{j};P^{-j},\theta_{f})=\sum_{a^{-j}\in\mathcal{A}^{|\mathcal{J}|-1}}P^{-j}(x,a^{-j})f(x'\mid x,a^{j},a^{-j};\theta_{f}) \]
Now consider the vector of player $j$'s (expected) choice-specific value functions, $v^{j}\in\mathbb{R}^{|\mathcal{X}|\times|\mathcal{A}|}$, and define the corresponding choice probabilities as $\Lambda^{j}(x,a^{j};v^{j})$, which is the probability agent $j$ chooses action $a^{j}$ in state $x$, conditional on having choice-specific value function $v^{j}$.\footnote{We note that the choice-specific value functions, $v^{j}$, are also often referred to as conditional value functions.} In equilibrium, the choice probabilities will be $P^{j}(a)=\Lambda^{j}(a;v^{j})$. And let \[ \Lambda^{-j}(v^{-j})=(\Lambda^{1}(v^{1}),\dots,\Lambda^{j-1}(v^{j-1}),\Lambda^{j+1}(v^{j+1}),\dots,\Lambda^{|\mathcal{J}|}(v^{|\mathcal{J}|})), \] so that in equilibrium $P^{-j}=\Lambda^{-j}(v^{-j})$. AgMir2007 show that in equilibrium, the choice-specific value functions are equal to \[ v^{j}(x,a^{j})=u^{j}(x,a^{j};P^{-j},\theta_{u})+\beta\sum_{x'}f^{j}(x'\mid x,a^{j};P^{-j},\theta_{f})\Gamma^{j}(x';\theta,P), \] where \[ \Gamma^{j}(\theta,P)=\left(I-\beta F(\theta_{f},P)\right)^{-1}\sum_{a^{j}}P^{j}(a^{j})*\left(u^{j}(a^{j};P^{-j},\theta_{u})+e(a^{j};P^{j})\right) \] maps into ex-ante (or integrated) value function space. Here, $e\left(a^{j};P^{j}\right)$ stacks the values of $e(x,a^{j};P^{j})\equiv E[\varepsilon^{j}(a^{j})\mid x,a^{j},P^{j}]$, as defined in Aguirregabiria and Mira (AgMir2007, Equation 11), and $F(\theta_{f},P)$ is an unconditional state transition matrix with elements $F(\theta_{f},P)\{k,l\}=\sum_{a\in\mathcal{A}^{J}}\left(\Pi_{j=1}^{J}P^{j}(a^{j}\mid x=k)\right)f(x'=l\mid x,a;\theta_{f})$.\footnote{See footnote 6 in AgMir2007.} They then define the NPL operator, $\Psi(\theta,P)$, such that $\Psi^{j}(\theta,P)=\Lambda^{j}(\Gamma^{j}(\theta,P))$ and combining all players yields the following fixed-point condition that describes any Markov perfect equilibrium (Aguirregabiria and Mira AgMir2007, Lemma 1): \[ P=\Psi(\theta,P). \]
While this equilibrium representation based on the NPL operator is often useful, we will ultimately want to work with an alternative representation when implementing our new estimator. This alternative arises due to a change of variables from $P$ space to $v$ space. (See Section (ref) for a detailed discussion of the importance of this change.) Define the function
where $\Phi:\Theta\times\mathbb{R}^{|\mathcal{J}|\times|\mathcal{X}|\times|\mathcal{A}|}\to\mathbb{R}^{|\mathcal{J}|\times|\mathcal{X}|\times|\mathcal{A}|}$ and $S(\cdot)$ is McFadden's social surplus function.\footnote{For example, $S(v^{j}(x))=\ln(\sum_{a^{j}}\exp(v^{j}(x,a^{j}))+\bar{\gamma}$ when the private values are i.i.d. and follow the type 1 extreme value distribution, where $\bar{\gamma}$ is the Euler-Mascheroni constant.} This $\Phi(\cdot)$ function allows us to characterize the equilibrium with an alternative fixed-point equation, as described in the following Lemma.
This section describes the $k$-EPL estimator and discusses its asymptotic and finite-sample properties, as well as computational aspects of its implementation in dynamic discrete choice games. We begin by discussing maximum likelihood estimation, subject to an equilibrium constraint based on some nuisance parameter, $Y$. The model is parameterized by a finite-dimensional vector, $\theta\in\Theta\subset\mathbb{R}^{|\Theta|}$, and a constraint $G(\theta,Y)=0$ where $Y\in\mathcal{Y}\subset\mathbb{R}^{|\mathcal{Y}|}$ and $G:\Theta\times\mathcal{Y}\to\mathbb{R}^{|\mathcal{Y}|}$. The true parameter values are $\theta^{*}$ and $Y^{*}$, with $G(\theta^{*},Y^{*})=0$. Note that there may be other values of $Y$ satisfying the constraint at $\theta^{*}$, but we will assume that the data are generated from only one such value, a common assumption in the literature.
The alternative statements of the equilibrium conditions in Section (ref) yield two potential choices for $Y$ and the corresponding constraint. Our asymptotic consistency and efficiency results in this section do not depend on the choice of nuisance parameter, but this choice can have serious implications for the computational implementation of the $k$-EPL estimator. So, for computational purposes we ultimately use the choice-specific value functions as the nuisance parameter ($Y\equiv v$), and the constraint comes from the equation in Lemma (ref) ($v-\Phi(\theta,v)=0$). We discuss the computational concerns in more detail in Section (ref). But for now, we present the asymptotic theory with the generic nuisance parameter, $Y$.
Let $w_{i}$ for $i=1,\dots,N$ denote the observations from $N$ independent markets, and define
where $q_{i}(\theta,Y)\equiv\ln Pr(w_{i}\mid\theta,Y)$. Furthermore, let $Q^{*}(\theta,Y)\equiv\mathrm{E}[Q_{N}(\theta,Y)]$.
Assumptions (ref)(a)-(c) echo standard identification assumptions. We note that assuming that $Q^{*}$ has a unique maximum does not rule out games with multiple equilibria. Non-singularity of the Jacobian in (d) is the defining feature of regular Markov perfect equilibria in the sense of DorEscobar2010. Regularity essentially means that the equilibrium is locally isolated and we can apply the implicit function theorem to obtain $Y(\theta)$ locally.\footnote{AgMir2007 directly assume the local existence of $Y(\theta)$, instead of appealing to the implicit function theorem.}
One method to estimate these models is via constrained maximum likelihood:
An equivalent statement is $\hat{\theta}_{\text{MLE}}=\textrm{arg max}_{\theta}\quad Q_{N}(\theta,Y(\theta))$, where $G(\theta,Y(\theta))=0$ and $\hat{Y}_{\text{MLE}}=Y(\hat{\theta}_{\text{MLE}})$. Pseudo-likelihood estimation replaces $Y(\theta)$ with some other mapping. AgMir2007 define $Y\equiv P$ and replace $Y(\theta)\equiv P(\theta)$ with their $\Psi(\theta,\hat{P}_{k-1})$ for the $k$-th iteration in the $k$-NPL sequence. However, this procedure suffers from the issues discussed in Section (ref).
Our $k$-step Efficient Pseudo-Likelihood ($k$-EPL) sequence instead uses a “Newton-like” step, which provides a good approximation to the full Newton step but uses a fixed value of $\nabla_{Y}G(\theta,Y)^{-1}$; this value varies between steps but does not vary as the optimizer searches over different values $\theta$ within each step. Algorithm (ref) below defines our sequential estimation procedure. It uses our Newton-like mapping, $\Upsilon(\cdot)$, which is a function of an initial compound parameter vector, $\gamma=(\theta,Y)$, and an additional (possibly different) value of $\theta$.
This procedure enjoys some nice econometric properties, both asymptotically and in finite samples. These properties arise from some convenient features of the $\Upsilon(\cdot)$ function, which are detailed in the following lemma.
Lemma (ref) is the key to most of the results in this section, which arise from applying the lemma at $(\theta^{*},\gamma^{*})$ and $(\hat{\theta}^{MLE},\hat{\gamma}^{MLE})$. For now, we note that Result 3 of Lemma (ref) is analogous to the “zero Jacobian” property from Proposition 2 of AgMir2002, which was the key to both their efficiency results and the finite-sample convergence results of KasShim2008 for single-agent $k$-NPL. By utilizing Newton-like steps on the equilibrium constraint, the $k$-EPL algorithm restores this zero Jacobian property in dynamic games.
One notable difference between $k$-EPL and some other sequential estimators is that $k$-EPL is initialized from a consistent estimate of both the structural parameter and the nuisance parameter, while other estimators -- including $k$-NPL -- only require an initially consistent estimate of the (nuisance) CCPs. Other examples include the finite-dependence-based estimators (HotzMil1993,ArcidMiller2011) and inequality-based estimators (BajBenkLevin2007). While none of these other estimators offer efficiency or convergence guarantees as $k\to\infty$, they can be used to obtain a consistent parameter estimate for initializing $k$-EPL. If the model exhibits finite-dependence, then finite-dependence-based estimators can be particularly attractive for obtaining initial estimates because of there computational simplicity. However, we note that the example models in our Monte Carlo simulations and empirical application may not have the finite dependence property because exit is not permanent. So, we instead use a $1$-NPL estimate for initialization.\footnote{Finite dependence is more difficult to establish in entry/exit games without a permanent exit decision that leads directly to a terminal state property. See, e.g., the discussion in Section 5.2 of ArcidEllick2011. However, ArcidiaconoMiller2019 show how to derive finite dependence in a class dynamic games that does not require a terminal state.}
One implication of Lemma (ref) is that $k$-EPL then gives a sequence of asymptotically efficient estimators that converges almost surely in large samples. We state this result formally in the following Theorem.
The results of Theorem (ref) for $k$-EPL in games are shared by $k$-NPL in single-agent models, but not games. In short, the zero Jacobian property ensures that $\hat{\gamma}_{k}=(\hat{\theta}_{k},\hat{Y}_{k})$ is asymptotically orthogonal to $\hat{\gamma}_{k-1}$, so that using $\hat{\gamma}_{k-1}$ is asymptotically equivalent to using $\gamma^{*}=(\theta^{*},Y^{*})$ at each step. This drives the consistency (Result 1) and asymptotic equivalence to MLE (Result 2) of each step. Intuitively, an EPL step is similar to a Newton step on the full maximum likelihood problem, although iterating on that procedure is notoriously unstable unless properly regularized.\footnote{Train2009 discusses the need to alter the step size in order to obtain global convergence (pp. 189-191). Nesterov2004 also discuses divergence (Section 1.2).} Our Monte Carlo simulations in Section ((ref)) show that $k$-EPL is stable without further regularization.
While the asymptotic distribution is insensitive to iteration, we can obtain substantial finite-sample improvements by iterating and can even compute the finite-sample MLE by iterating to convergence. In this section, we proceed with a formal econometric analysis of the local convergence rate of the iterations to the maximum likelihood estimator, and we discuss the implications for finite sample performance. Later on, the finite sample properties are illustrated in the Monte Carlo simulations in Section (ref).
Our results for the convergence rate to MLE are similar to those of the single-agent version of $k$-NPL in AgMir2002, which also has the zero Jacobian property. In their Monte Carlo simulations, single-agent $k$-NPL iterations exhibit rapid finite-sample improvements and reliably converged to the finite-sample MLE. KasShim2008 then provided a formal econometric explanation for these results. The analysis of $k$-EPL's finite sample properties in this section is similar but also applies to games.
The only additional requirement for our finite-sample results is that the Jacobian of the equality constraints, $G$, with respect to $Y$ is nonsingular at the finite-sample MLE.
Assumption (ref) guarantees the existence of an implicit function, $Y(\theta)$, around $\hat{\theta}_{MLE}$ and also that the quasi-Newton mapping, $\Upsilon(\theta,\hat{\gamma}_{MLE})$ is valid. Assumption (ref) is enough to guarantee that Assumption (ref) is satisfied almost surely as $N\to\infty$, since it implies that $\det\left(\nabla_{Y}G\left(\hat{\theta}_{MLE},\hat{Y}_{MLE}\right)\right)\overset{a.s.}{\to}\det\left(\nabla_{Y}G\left(\theta^{*},Y^{*}\right)\right)\neq0$ by continuity of $\det(\cdot)$, continuity of $\nabla_{Y}G(\cdot)$ (Assumption 5(d)), and strong consistency of $\hat{\gamma}_{MLE}$. Furthermore, the set of singular matrices has measure zero. So, we view this as a relatively mild (but important) assumption.
The first result of Theorem (ref) establishes that the MLE is a fixed point of the $k$-EPL iterations in a finite sample, similar to Aguirregabiria and Mira (AgMir2002, Proposition 3) for single-agent $k$-NPL. The second result gives an asymptotic analysis of convergence to MLE, which provides a theoretical explanation of why we should expect iteration to yield improvements in finite samples. This result is analogous to Proposition 2 of KasShim2008, although their result only applies in the single-agent case. In short, even though iteration on EPL provides no improvement up to $O_{p}(N^{-1/2})$, it still yields higher-order improvements. To see why, suppose the initial estimates are such that $\hat{\gamma}_{0}-\gamma^{*}=O_{p}(N^{-b})$ for $b\in(1/4,1/2]$, so that $||\hat{\gamma}_{0}-\hat{\gamma}_{\text{MLE}}||=O_{p}(N^{-b})$.\footnote{For $b\in(1/4,1/2]$, we have $\hat{\gamma}_{0}-\hat{\gamma}_{\text{MLE}}=\hat{\gamma}_{0}-\gamma^{*}-(\hat{\gamma}_{\text{MLE}}-\gamma^{*})=O_{p}(N^{-b})+O_{p}(N^{-1/2})=O_{p}(N^{-b})$ . Additionally, this can be used to show higher-order equivalence to the MLE.} Repeated substitution gives $||\hat{\gamma}_{k}-\hat{\gamma}_{\text{MLE}}||=O_{p}(N^{-(k-1)/2-2b})$. In particular, in the case where the state space is finite and frequency or kernel estimates are used, $b=1/2$ and $||\hat{\gamma}_{k}-\hat{\gamma}_{\text{MLE}}||=O_{p}(N^{-(k+1)/2})$, where $N^{-(k+1)/2}\to0$ as $k\to\infty$ for $N>1$. Our own Monte Carlo simulations in Section (ref) exhibit such improvements.
The third result in Theorem (ref) allows us to consider EPL iterations as a computationally attractive algorithm for computing the MLE. It establishes that we can expect the EPL iterations to be a local contraction around the MLE in the finite sample with a very fast convergence rate. The full proof appears in the appendix, but it essentially proceeds by noting that the $k$-EPL sequence satisfies $\hat{\gamma}_{k}=H_{N}(\hat{\gamma}_{k-1})$, where $\hat{\gamma}_{MLE}$ is a fixed point of the function $H_{N}$. And due to the zero-Jacobian property in Lemma (ref) (Result 3), we obtain $\nabla_{\gamma}H_{N}(\hat{\gamma}_{MLE})\overset{a.s.}{\to}0$.\footnote{This drives the Neyman orthogonality discussed in the introduction.}
For the population analogue of the EPL iterations, the convergence rate is then super-linear. However, we only have access to finite samples in practice, so we should expect the convergence rate to be linear with a small Lipschitz constant, implying that we'll need only a few iterations to achieve convergence. We can therefore use EPL iterations to compute the MLE even when a consistent $\hat{\gamma}_{0}$ is unavailable. We can simply use multiple starting values, iterate to convergence, and use the converged estimate that provides the highest log-likelihood.\footnote{AgMir2007 suggest a similar procedure to find $\infty$-NPL when no initial consistent estimate is available, and multiple starting values are often used when computing the maximum likelihood estimate with other methods (SuJudd2012,EgLaiSu2015).} We demonstrate this usage of $k$-EPL with Monte Carlo simulations in the appendix and find that it works well.
Aside from $k$-EPL, there are two potential alternative algorithms for computing the MLE: the nested fixed-point (NFXP) algorithm \'a la RustGMC1987 and the MPEC approach proposed by SuJudd2012 and extended to dynamic games by EgLaiSu2015. The NFXP algorithm searches over $\theta$ in an outer loop and finds $Y_{\theta}$ such that $G(\theta,Y_{\theta})=0$ in an inner loop. MPEC leverages modern optimization software to search over $\theta$ and $Y$ simultaneously, only imposing that $G(\theta,Y)=0$ at the solution. The algorithm of choice may depend on the structure of the model.
While this section discusses the $k$-EPL algorithm in the context of a general constrained maximum-likelihood problem, we are ultimately focused on estimating dynamic discrete choice games of incomplete information. As discussed in the introduction, NFXP is often computationally unattractive---or even infeasible---in such games.\footnote{We note that NFXP still performs well in single-agent dynamic models. See DorJudd2012 and ABBE2016 for details on the computational burden of computing equilibria in discrete-time dynamic discrete games.} MPEC, however, remains feasible and performs well, as demonstrated by EgLaiSu2015. The key difference here between MPEC and $k$-EPL is that $k$-EPL will be able to more heavily exploit the structure of the problem.\footnote{EgLaiSu2015 exploit sparsity patterns in their MPEC implementation but do not further exploit other features of the problem structure.} In Section (ref), we show that --- much like $k$-NPL in single-agent models --- common modeling assumptions lead to EPL iterations composed of two easily-computed parts: solving linear systems to form pseudo-regressors, followed by solving an unconstrained, globally concave maximization problem \'a la static logit/probit with those pseudo-regressors. Neither of these operations require sophisticated commercial optimization software; and repeating them just a few times may ultimately be more computationally attractive than using MPEC to simultaneously solve for all variables in a non-concave, large-scale, constrained maximization problem.
$k$-EPL is particularly useful when we are interested in estimating the flow utility parameters, $\theta_{u}$. In many cases -- including our Monte Carlo experiments -- the transition parameter, $\theta_{f}$, is known in advance. So, we focus on estimating only the flow utility parameters and let $\theta\equiv\theta_{u}$. Alternatively, $\theta_{f}$ can be estimated in a first stage, with $\theta_{u}$ then estimated via partial maximum likelihood. Similarly to single-agent $k$-NPL, $k$-EPL iterations based on this partial MLE problem yield asymptotic equivalence and finite-sample convergence to partial MLE because the zero Jacobian property still applies.
Assumption (ref) requires the flow utilities to be linear in $\theta$, which is a standard assumption for dynamic discrete choice models.\footnote{See, e.g., RustGMC1987,AgMir2002,AgMir2007,BajBenkLevin2007,PakesOstBerry2007,PesSD2008,ArcidMiller2011,EgLaiSu2015,BugniBunt2020.} Assumption (ref) requires log-concavity of the mapping from choice-specific values into CCPs. A sufficient condition for this is that the distribution of private shocks, $g(\cdot)$, is log-concave (CaplinNalebuff1991). Consequently, Assumption (ref) is satisfied in the ubiquitous case of logit shocks, as well as when shocks follow a normal distribution. These two assumptions have important computational implications, which is the focus of the rest of this section.
The choice of nuisance parameter, $Y$, does not affect the asymptotic results in Section (ref), but it does have tremendous implications for computation. In order to provide a computationally simple estimator, we use the choice-specific values as the nuisance parameter and the equilibrium condition from Lemma (ref).
Coupled with Assumptions (ref) and (ref), the choice of nuisance parameter and constraint in Assumption (ref) leads to some convenient computational properties.
By Assumption (ref), we have $u^{j}(\theta,x,a^{j},P^{-j})=h(x,a^{j},P^{-j})'\theta$. Because $P^{-j}=\Lambda^{-j}(v^{-j})$, we can re-write this in terms of $v$: $u^{j}(\theta,x,a^{j},v^{-j})=h(x,a^{j},v^{-j})'\theta$. Inspecting the form of $\Phi(\theta,v)$, we see that it will be linear in $\theta$ and therefore so will $G(\theta,v)=v-\Phi(\theta,v)$: \[ G(\theta,v)=H(v)\theta+z(v), \] where $H(\cdot)$ is a matrix and $z(\cdot)$ is a vector. As a result, $\Upsilon(\theta,\hat{\gamma}_{k-1})$ is also linear in $\theta$:
Additionally, the optimization step in Algorithm (ref) ($k$-EPL) becomes \[ \hat{\theta}_{k,EPL}=\underset{\theta\in\Theta}{\textrm{arg max}}\quad N^{-1}\sum_{i=1}^{N}\sum_{t}\sum_{j}\ln\Lambda\left(\Upsilon(x_{t},a_{t}^{j};\theta,\hat{\gamma}_{k-1})\right). \] It turns out that this is a concave optimization problem, as described in the next proposition.
Proposition (ref) shows how our choice of nuisance parameter and constraint lead to a computationally simple estimation sequence. Computing $A(\hat{\gamma}_{k-1})$ and $b(\hat{\gamma}_{k-1})$ in Proposition (ref) requires computing $\nabla_{v}G(\hat{\theta}_{k-1},\hat{v}_{k-1})^{-1}H(\hat{v}_{k-1})$ and $\nabla_{v}G(\hat{\theta}_{k-1},\hat{v}_{k-1})^{-1}z(\hat{v}_{k-1})$, respectively, which are the solutions to linear systems. Importantly, these linear systems can be solved outside the optimization search in Step 2 of Algorithm (ref). So, the computation procedure alternates between i) computing “pseudo-regressors” by solving linear systems; and ii) maximizing a concave optimization problem that using the pseudo-regressors as inputs. Furthermore, $\nabla_{v}G(\hat{\theta}_{k-1},\hat{v}_{k-1})$ can be computed analytically when $\Lambda(\cdot)$ and its derivative have an analytic form, such as the logit or probit cases.
Notably, this computational simplicity is not available if the CCPs are chosen as the nuisance parameter. That is, when $Y\equiv p$ and $G(\theta,P)\equiv P-\Psi(\theta,P)$. In this case, \[ \Upsilon(\theta,\hat{\gamma}_{k-1})=\hat{P}_{k-1}-\left(I-\nabla_{P}\Psi(\hat{\theta}_{k-1},\hat{P}_{k-1})\right)^{-1}\left(P-\Psi(\theta,\hat{P}_{k-1})\right) \] and the optimization step in the $k$-EPL algorithm solves \[ \hat{\theta}_{k,\text{EPL}}=\underset{\theta\in\Theta}{\textrm{arg max}}\quad N^{-1}\sum_{i=1}^{N}\sum_{t}\sum_{j}\ln\Upsilon(x_{t},a_{t}^{j};\theta,\hat{\gamma}_{k-1}). \] Several computational issues arise. First, instead of solving linear systems once before the optimization step, we must repeatedly solve linear systems throughout the optimization because of the need to compute $(I-\nabla_{P}\Psi(\hat{\theta}_{k-1},\hat{P}_{k-1}))^{-1}\Psi(\theta,\hat{P}_{k-1})$ for each new value of $\theta$ in the search. Second, we will lose the guarantee of concavity of the optimization problem in each step. Even though $\Psi(\theta,\hat{P}_{k-1})$ is log-concave in $\theta$, this does not guarantee log-concavity of $\Upsilon(\theta,\hat{\gamma}_{k-1})$ because affine transformations of log-concave functions are not necessarily log-concave. And third, the Newton-like steps rely on an implicit linearization: even though $\Psi(\theta,\hat{P}_{k-1})$ maps into the probability simplex, $\Upsilon(\theta,\hat{\gamma}_{k-1})$ can arrive at values outside the simplex.\footnote{Technically, $\Psi(\cdot)$ maps into a Cartesian product of the interior of the unit simplex due to each player having their own strategies.} Thus, we would need to add constraints to the optimization problem to ensure $\Upsilon(\theta,\hat{\gamma}_{k-1})$ does not leave the unit simplex. Whereas, the formulation in $v$-space does not require any constraints because $v$ can take any value on a Cartesian product of the real line.
While our EPL iterations with $Y\equiv v$ have a similar computational structure to NPL iterations (with $Y\equiv P$) insofar as both require solving linear systems then a globally concave optimization problem, the dimension of the linear systems in $k$-EPL is larger. AgMir2007 show that $k$-NPL requires solving $|J|(|\Theta|+1)$ different systems of linear equations, each of dimension $|\mathcal{X}|$ , resulting in a worst-case bound of $O((|\Theta|+1)|\mathcal{X}|^{3}|\mathcal{J}|)$ flops.\footnote{There are $|\Theta|+1$ systems for each of the $|\mathcal{J}|$ players.} On the other hand, $k$-EPL requires solving $|\Theta|+1$ different systems of linear equations, each of dimension $|\mathcal{J}||\mathcal{X}||\mathcal{A}|$, resulting in a larger worst-case bound of $O((|\Theta|+1)|\mathcal{X}|^{3}|\mathcal{J}|^{3}|\mathcal{A}|^{3})$ flops. Sparsity of the linear systems -- a common feature in dynamic discrete choice models (EgLaiSu2015) -- can lower these bounds for both $k$-NPL and $k$-EPL. Fortunately, in practice, in our largest Monte Carlo experiment the relative difference is much lower than suggested by the worst-case bounds.
We then see a tradeoff between efficiency and computational burden, a common theme in estimating dynamic games. This theme appears in Bugni and Bunting's BugniBunt2020 analysis comparing their efficient $k$-MD (minimum distance) estimator to $k$-NPL. Even in what is essentially the smallest scale game possible -- 2 players, 2 actions, 4 states -- they report for $k$-MD a large increase in computational burden over $k$-NPL, with individual iterations taking about 12 to 26 times longer on average, depending on the number of iterations. It is perhaps reasonable to expect that this difference will grow with the size of the game, which is a concern for practitioners who must balance econometric efficiency with computational feasibility. Our $k$-EPL estimator, on the other hand, induces a much less severe increase in computational burden while retaining efficiency. In Section (ref), we also explore a 2 player, 2 action, 4 state game and find that the difference between computational time for $k$-EPL and $k$-NPL iterations is negligible in that setting. Additionally, we explore a much larger-scale game -- 5 players, 2 actions, 160 states -- that is more representative of empirically relevant models. In this larger-scale game, we find that there is indeed an increase in time per iteration for $k$-EPL relative to $k$-NPL, but this increase in the larger-scale game (between 4 to 8 times) is not even as large as the 12 to 26-fold increase for $k$-MD in the much smaller-scale game.
Even with an increase in computation time per iteration relative to $k$-NPL, $k$-EPL can still ultimately be more attractive than $k$-NPL. First, its asymptotic efficiency, convergence properties, and rapid finite-sample improvements are attractive features that may be worth the increased computational burden of each iteration. Second, even in cases where both $k$-EPL and $k$-NPL converge to consistent estimates, $k$-EPL enjoys a much faster convergence rate than $k$-NPL, resulting in fewer iterations to convergence. So, iterating to convergence on $k$-EPL to obtain the finite-sample MLE can still be faster than computing the $\infty$-NPL estimator (if it converges), even though each individual iteration takes longer.
In many applications, the dominant source of computational burden for either estimator will often be the size of the state space, $|\mathcal{X}|$, since it can be large when $|\mathcal{A}|$ and $|\mathcal{J}|$ are small and also tends to grow with both $|\mathcal{A}|$ and $|\mathcal{J}|$ in dynamic games. One simple yet salient illustration arises when the state is determined by the previous actions of the players, so that $|\mathcal{X}|=|\mathcal{A}|^{|\mathcal{J}|}$. Thus, as $|\mathcal{J}|$ grows, $|\mathcal{X}|$ ultimately becomes the main source of computational burden for the linear systems in both $k$-NPL and $k$-EPL. To help deal with large state spaces, AgMir2007 show that the linear systems required for $k$-NPL can be solved via an iterative process reminiscent of value function Bellman iteration, so that their worst-case computational burden reduces to $O((|\Theta|+1)|\mathcal{X}|^{2}|\mathcal{J}|)$. Similarly, $k$-EPL can also use alternative iterative methods to solve the linear systems such as Krylov subspace methods, although it cannot use Bellman-style iteration.
We conclude this section by showing that $k$-NPL in a single-agent dynamic discrete choice model (AgMir2002) is equivalent to $k$-EPL with a slightly modified definition of $\Upsilon(\cdot)$. Here, we can work directly in probability space, $Y\equiv P$, and let \[ G(\theta,P)=P-\Psi(\theta,P), \] with $\Psi(\theta,P)$ defined as in Section (ref) but with only a single agent.
We now have \[ \nabla_{P}G(\theta,P)=I-\nabla_{P}\Psi(\theta,P). \]
Proposition 2 from AgMir2002 shows that $\nabla_{P}\Psi(\theta,P_{\theta})=0$, where $P_{\theta}=\Psi(\theta,P_{\theta})$. Thus, $\nabla_{P}G(\theta,P_{\theta})=I$ for all $\theta$. So, we can use a modified definition of $\Upsilon(\cdot)$, where $\nabla_{P}G(\hat{\theta}_{k-1},\hat{P}_{k-1})$ is simply replaced with the identity matrix, $I$, and we obtain
This modified implementation of $k$-EPL is identical to $k$-NPL.
This equivalence of $k$-NPL to $k$-EPL in single agent models is unsurprising for two reasons. First, we stated in the introduction that the motivation for $k$-EPL is to extend the nice properties of $k$-NPL from single-agent models to dynamic games. So there should be, at the very least, substantial conceptual overlap between the methods. Second, Aguirregabiria and Mira (AgMir2002, Proposition 1(c)) show that their policy iterations are equivalent to Newton-like iterations on the (ex-ante) value function in single-agent models. Since $k$-EPL is built around Newton iterations, such an equivalence is again suggestive of the relationship shown here.
In this section, we present Monte Carlo simulation results to illustrate $k$-EPL's finite sample properties. The simulations presented here are based on the dynamic game of entry and exit in Example (ref), parameterized to match a model with five heterogeneous firms from AgMir2007. The appendix includes further Monte Carlo simulations for two other models: (i) a small-scale dynamic model from PesSD2008; and (ii) a static game from PesSD2010. The dynamic model in PesSD2008 exhibits multiple -- possibly best-reply-unstable -- equilibria (with data generated from only one of them), which can be challenging for other iterative methods. The static game in PesSD2010 provides a setting where we can easily compare $\infty$-EPL to the MLE computed via the nested fixed-point algorithm.
The simulations in this section are based on an empirically relevant model that forms the basis of many applications but has also been used (sometimes in simplified forms) in simulation studies by KasShim2012, EgLaiSu2015, BugniBunt2020, BlevKim2019, and AgMarcoux2019. In particular, AgMarcoux2019 discuss in detail how the spectral radius of the NPL operator increases with the strength of competition in the model. As such, we take as our baseline case the parameters as Experiment 2 of AgMir2007. We then follow AgMarcoux2019 and increase the competitive effect parameter, $\theta_{\text{RN}}$ to investigate how $k$-EPL and $k$-NPL behave as the spectral radius of the NPL operator increases to the point of instability and beyond.
The model is a dynamic entry-exit game with $\lvert\mathcal{J\rvert}=5$ firms that operate in $N$ independent markets. The firms have heterogeneous fixed costs. There is a single common market state, the market size, which can take one of five values: $s_{it}\in\lbrace1,2,3,4,5\rbrace$. Market size follows a $5\times5$ transition matrix and we use the same transition matrix as AgMir2007. The other observable states are the incumbency statuses of the five firms, denoted $a_{i,t-1}^{j}$. There are therefore $5\times2^{5}=160$ distinct states in the model. Hence the state in market $i$ at time $t$ can be represented in vector form as $x_{it}=(s_{it},a_{i,t-1}^{1},a_{i,t-1}^{2},a_{i,t-1}^{3},a_{i,t-1}^{4},a_{i,t-1}^{5})$.
Given the state of the model at the beginning of the period, firms simultaneously choose whether to operate in the market, $a_{it}^{j}=1$, or not, $a_{it}^{j}=0.$ They make these decisions in order to maximize expected discounted profits, where the period profit function for an active firm is
\[ \bar{u}^{j}(x_{it},a_{it}^{j}=1,a_{it}^{-j};\theta)=\theta_{\text{FC},j}+\theta_{\text{RS}}s_{it}-\theta_{\text{RN}}\ln\left(1+\sum_{l\neq j}a_{it}^{l}\right)-\theta_{\text{EC}}(1-a_{i,t-1}^{j}) \] and $\bar{u}^{j}(x_{it},a_{it}^{j}=0,a_{it}^{-j};\theta)=0$ for inactive firms. The game is dynamic because firms must pay an entry cost $\theta_{\text{EC}}$ to enter the market and because firms have forward-looking expectations about market size and entry decisions of rival firms. The private information shocks $\varepsilon_{it}^{j}(a_{it}^{j})$ are independent and identically distributed across time, markets, players, and actions and follow the standard type I extreme value distribution.
We choose the model parameters following AgMir2007 and AgMarcoux2019. In particular, the fixed costs for the five firms are $\theta_{\text{FC},1}=-1.9$, $\theta_{\text{FC},2}=-1.8$, $\theta_{\text{FC},3}=-1.7$, $\theta_{\text{FC},4}=-1.6$, and $\theta_{\text{FC},5}=-1.5$. The coefficient on market size is $\theta_{\text{RS}}=1$ and the common firm entry cost is $\theta_{\text{EC}}=1$. The only parameter that differs across our three experiments is the competitive effect $\theta_{\text{RN}}$, which we set to be $\theta_{\text{RN}}=1$ in Experiment 1, $\theta_{\text{RN}}=2.5$ in Experiment 2, and $\theta_{\text{RN}}=4$ in Experiment 3. For easy comparison, these parameter values correspond closely to the “very stable,” “mildly unstable,” and “very unstable” cases of AgMarcoux2019.
For each experiment, we carry out $1000$ replications using two sample sizes, $N=1600$ and $N=6400$, noting that $N=1600$ is the sample size used by AgMir2007. For each replication and each sample size, we draw a sample of size $N$. With this sample we calculate the iterative $k$-NPL and $k$-EPL estimates. For $k$-NPL, we follow the original AgMir2007 implementation. We initialize $k$-NPL with estimated semiparametric logit choice probabilities. We then initialize $1$-EPL using the parameter estimates and value function from the $1$-NPL iteration, which is consistent even in cases where further iteration may lead to inconsistency. Aside from the initialization, and using the same sample, the $k$-EPL iterations proceed independently from the $k$-NPL iterations. For $k$-EPL, as before we represent the equilibrium condition in terms of $v$ as $G(\theta,v)=v-\Phi(\theta,v)$ and use analytical derivatives for the Jacobian $\nabla_{v}G(\theta,v)$.
We report estimates from one, two, and three iterations of each estimator as well as the converged values, which we denote as $\infty$-NPL and $\infty$-EPL. We limit the number of iterations to 100 for both estimators, and if the algorithm has not converged we use the estimate from the final iteration. We use the same convergence criteria for both estimators: at each iteration we check the sup norm of the change in the parameter values and choice probabilities. Based on the criteria used by AgMir2007, if both are below $10^{-2}/K$, where $K$ is the number of parameters, we terminate the iterations and return the final converged estimate.
Figure\ (ref) shows the Monte Carlo distributions of two key parameter estimates in the model: $\hat{\theta}_{\text{RN}}$, the competitive effect, and $\hat{\theta}_{\text{EC}}$, the entry cost. Each panel compares two histograms: the histogram above in dark gray corresponds to $\infty$-EPL and the histogram below in light gray is for $\infty$-NPL. These histograms are based on data for 1,000 replications of each experiment with 6,400 observations.
The left panels of Figure\ (ref) show the distributions of $\hat{\theta}_{\text{RN}}$, which is related to the strength of competition in the market and is closely related with the spectral radius of the NPL operator (while the spectral radius of the EPL operator is always zero). Indeed, we can see that as the competitive effect becomes large the distribution of $\infty$-NPL estimates is truncated at around $\hat{\theta}_{\text{RN}}=2.5$ in Experiment 2. It remains concentrated around $\hat{\theta}_{\text{RN}}=2.6$ in Experiment 3 even though the value used in the data generating process was $\theta_{\text{RN}}=4$. Yet for all experiments the distributions of $\infty$-EPL estimates are centered around the true values.
Although the true entry cost parameter is fixed at $\theta_{\text{EC}}=1$ in all three experiments, the distribution of estimates can be affected when we vary the competitive effect from $\theta_{\text{RN}}=1$ to $\theta_{\text{RN}}=4$. In the right panels of Figure\ (ref), the truncation from above of $\hat{\theta}_{\text{RN}}$ in the $\infty$-NPL distribution induces truncation from below in $\hat{\theta}_{\text{EC}}$, as lower estimated values of competition result in higher estimated entry costs, leading to distortion of both distributions and to parameter estimates that would lead a policy-maker to possibly very different economic implications. In Experiment 3, although the distribution of $\infty$-NPL estimates again appear to be normally distributed for both parameters, they are biased and are not centered around the true values.
Figure\ (ref) shows the distributions of iteration counts and computational times for $\infty$-EPL and $\infty$-NPL across the three experiments.\footnote{Computational times are reported using Matlab R2020b on a 2019 Mac Pro with a 3.5 GHz 8-Core Intel Xeon W processor.} The histograms in the left panels show the distribution of iteration counts required to achieve convergence for both estimators. Recall that the maximum number of iterations allowed was 100. $\infty$-EPL requires fewer iterations for all experiments, especially for Experiments 2 and 3 where $\infty$-NPL sometimes fails to converge in Experiment 2 and always fails to converge in Experiment 3. $\infty$-EPL converged for all replications in all experiments.
The histograms in the right panels of Figure\ (ref) show the distribution of total computational time, in seconds, across the 1000 replications. In Experiment 1, where $k$-NPL is stable, $\infty$-NPL is faster even though it requires more iterations on average. In other words, each $k$-EPL iteration is more expensive on average, but fewer are required. This is in line with the analysis in Section (ref). However, in Experiments 2 and 3 the computational times for $\infty$-NPL increase as it requires more iterations, eventually overtaking the time required for $\infty$-EPL yet still frequently (Experiment 2) or always (Experiment 3) failing to converge.
Summary statistics for all parameter estimates and all three experiments can be found in Tables (ref)--(ref). Each table reports summary statistics for 1-, 2-, 3-, and $\infty$-NPL and 1-, 2-, 3-, and $\infty$-EPL over the 1000 Monte Carlo replications for $N=1600$ or $N=6400$. The upper two panels report the mean bias and mean square error (MSE) for each parameter and each sequential estimator. The next panel reports summary statistics for the number of iterations completed until either convergence or failure at 100 iterations. This includes the median, maximum, and inter-quartile range (IQR) of iteration counts across the replications as well as the number of replications for which $\infty$-NPL failed to converge. Finally, the bottom panel reports the computational times for each estimator including the mean and median total time (for all completed iterations) across the 1000 replications and the median time per iteration.
For Experiment 1, both the $k$-NPL and $k$-EPL estimators perform equally well with low bias, as can be seen in Tables (ref) and (ref). The parameter with the most finite-sample variation is also the main parameter of interest in our study: $\theta_{\text{RN}}.$ Note that in this model, in order to obtain estimates with performance similar to the converged estimates it would suffice to stop at 3 iterations with either estimator. Typically, both $\infty$-NPL and $\infty$-EPL converge in 4 to 5 iterations. However, even in this specification where $k$-NPL performs well, in one case out of 1000, $\infty$-NPL fails to converge in 100 iterations or less while $\infty$-EPL always converges in at most 7 iterations. In terms of computational time, in this model due to the computational complexity, in the median experiment one iteration of $k$-EPL is more expensive (0.3 seconds) than one iteration of $k$-NPL (0.045 seconds). Because roughly the same number of iterations are required in this model, the overall times for $\infty$-EPL are also longer than for $\infty$-NPL. However, even with the increased complexity each replication of $\infty$-EPL takes only around 1.3 seconds to estimate, so it remains quite feasible.
In Experiment 2, we begin to see a divergence between the two methods. As reported by AgMarcoux2019, the spectral radius of the population NPL operator is slightly larger than one for this specification. In finite random samples, sometimes the sample counterpart is stable and sometimes it is unstable. This leads to the situation illustrated by Table (ref), where $\infty$-NPL fails to converge in 612 of 1000 replications. $\infty$-EPL, on the other hand, is stable and converges for all 1000 replications. Importantly, the 1-NPL estimator obtained without further iterations is always consistent. However, there is substantial bias in the $\infty$-NPL estimates. In an apparent contradiction, the MSE for $\theta_{\text{RN}}$ is actually lower for $\infty$-NPL than for $\infty$-EPL. This pattern of larger bias and lower MSE is seen again with the large sample size in Table (ref), however it can be understood simply by recalling the histogram of the $\theta_{\text{RN}}$ estimates (panel (b) of Figure (ref)). The $\infty$-NPL sampling distribution appears to be truncated near 2.4, which perhaps not coincidentally is near the value where the spectral radius exceeds one (AgMarcoux2019). Since this happens to be close to the true parameter value, the MSE is artificially low. However, the sampling distribution is neither normally distributed nor centered at the true value. A K-S test for normality of the $\infty$-NPL estimates has a $p$-value equal to zero up to three decimal places, while for $\infty$-EPL the $p$-value is 0.622.
In Experiment 2 we also see a reversal of the order of computational times: the non-convergent $\infty$-NPL cases require more iterations and more time per iteration, while $\infty$-EPL always converges in 8 or fewer iterations. Thus, in thinking about the trade-off between robustness and computational time we should also consider convergence. A non-convergent estimator may take longer and yield worse results in the end. $k$-EPL does require more time per iteration in this model, but it is more robust to the strength of competition in the model.
Turning to Experiment 3, we increase the effect of competition even further to where $\theta_{\text{RN}}=4$, well beyond the point where $k$-NPL becomes unstable. In this case, both with small samples and large samples, $\infty$-NPL fails to converge in all 1000 replications but $\infty$-EPL converges in all 1000 replications. In these cases, the bias in $\infty$-NPL estimates is larger and in this case and so is the MSE, since the value of $\theta_{\text{RN}}$ is farther from the point of truncation than in Experiment 2. The $\infty$-NPL estimator systematically underestimates the competitive effect $\theta_{\text{RN}}$ and overestimates the entry cost $\theta_{\text{EC}}$.
Overall, across the three experiments the performance of $k$-EPL is stable and of similar quality despite the increasing competitive effect. This result agrees with our theoretical analysis showing that $k$-EPL is stable, convergent, and efficient.
Finally, we note that with the $k$-EPL estimator there is very little performance improvement after the first three iterations. The performance of $\infty$-EPL is achieved already, up to two decimal places, by 3-EPL. Thus, for this model one can reduce the computational time required while retaining efficiency and robustness by carrying out only a few iterations of $k$-EPL.
The U.S. wholesale club store industry is a retail segment that offers members a wide range of merchandise at discounted prices. The industry is dominated by three major players: Sam's Club, Costco, and BJ's Wholesale Club. These companies operate as membership-based clubs, with a focus on bulk purchasing and high-volume sales.
The modern wholesale club store industry emerged in the 1970s with the founding of Price Club in San Diego, California. Price Club was a pioneer in the industry, offering its members deep discounts on products, including groceries, electronics, and household goods. In 1983, Costco was founded in Seattle, Washington, and quickly became a major player in the industry. Sam's Club, a subsidiary of Walmart, was also founded in 1983. BJ's Wholesale Club was founded in 1984 in Massachusetts. In 1993, Price Club merged with Costco and adopted the Costco name (costco-about-us).
Today, the wholesale club store industry is a major force in retail, with the three main players generating over \$260 billion in annual revenue, according to data from their 2021 annual reports. Costco is the largest company in the industry, with over 800 stores worldwide and annual revenue of over \$190 billion (costco-2021). Sam's Club is the second-largest, with around 600 stores and annual revenue of over \$60 billion (sams-2021). BJ's Wholesale Club is the smallest of the three, with over 200 stores and annual revenue of over \$16 billion (bj-2021).
Although we focus on Costco, Sam's Club, and BJ's, which are by far the largest firms in the industry, there are other smaller players. An example is DirectBuy, which was the 4th largest firm (by total markets served) in our sample. While it was a significant player in the home furnishings and home improvement space, it was not a direct competitor of Costco, Sam's Club, and BJ's and was significantly smaller. DirectBuy faced financial difficulties and filed for bankruptcy in 2018 (furniture-today).
Our data come from the Data Axle Business Database which contains information on businesses across the United States. The company uses a variety of sources to gather information on businesses, including public records, government filings, and proprietary data sources. Our data begins in 2009 and ends in 2021. From this database we first extract records for each Sam's Club, Costco, and BJ's location. We augment this with ZIP code level population data obtained from NHGIS (manson2022ipums). We then aggregate to county-level markets using the 2021 ZIP code to county crosswalk files provided by the U.S. Department of Housing and Urban Development (HUD). We consider all counties in the 50 states of the United States and the District of Columbia that have between 20,000 and 600,000 residents. This serves to exclude very small counties that would clearly not be considered for entry as well as some very large, atypical markets.
Overall, our final sample consists of $N=1,600$ counties observed over $T=12$ years.\footnote{Although we have 13 years of data, we require lagged actions to construct the incumbency status state variable. As a result, the time dimension of our sample is reduced to $T=12$.} Among these markets, the average peak population during our sample was 104,841 with standard deviation 112,759.\footnote{We compute peak population for a market as the maximum over the years in the sample.} In the model, our market size variable, $s_{it}$, is the logarithm of population discretized into 5 equal bins.\footnote{We also tried 10 bins without any major substantive changes in the results.} Table (ref) presents summary statistics for our sample. On average, there are 0.348 active firms in each market with a standard deviation of 0.622. The autoregressive coefficient for the number of active firms is 0.987, indicating a strong positive correlation between the number of active firms in the current period and the previous period.\footnote{This refers to the estimated autoregressive coefficient in an AR(1) regression of the number of current active firms on the number of active firms in previous period.} The average number of entrants is 0.010, while the average number of exits is 0.006. Excess turnover, defined as (\#entrants + \#exits) - |\#entrants - \#exits|, is effectively zero. The correlation between entries and exits is -0.007. The probability of being active is highest for Sam's Club at 0.201, followed by Costco at 0.093 and BJ's at 0.054. The distribution of market size is such that there are relatively more small-to-medium size markets markets and relatively fewer large markets.
Our model of wholesale club competition follows the dynamic oligopoly model with heterogeneous firms described in Example (ref) and also used in our Monte Carlo experiments in Section (ref). In our application the firms are denoted $\mathcal{J}=\{\text{SC},\text{CC},\text{BJ}\}$. In each market $i=1,\dots,N$ and time period $t=1,\dots,T$, firms decide whether to operate in a market ($a_{it}^{j}=1$) or not ($a_{it}^{j}=0$). The profit function has the following form for an active firm $j$: \[ \bar{u}^{j}(x_{it},a_{it}^{j}=1,a_{it}^{-j};\theta)=\theta_{\text{FC},j}+\theta_{\text{RS}}s_{it}-\theta_{\text{RN}}\ln\left(1+\sum_{l\neq j}a_{it}^{l}\right)-\theta_{\text{EC}}(1-a_{i,t-1}^{j}). \] Here, $\theta_{FC,j}$ is the fixed cost of operation for firm j, $\theta_{EC}$ is the cost incurred by a new entrant, $\theta_{RS}$ represents the effect of market size $s_{it}$ (the discretized logarithm of population, defined above), and $\theta_{RN}$ captures the effect of competition. When firm $j$ is inactive, $\bar{u}^{j}(x_{it},a_{it}^{j}=0,a_{it}^{-j};\theta)=0$.
Using $k$-EPL, we estimate the heterogeneous fixed costs $\theta_{\text{FC},\text{SC}}$, $\theta_{\text{FC},\text{CC}}$, and $\theta_{\text{FC},\text{BJ}}$ as well as the entry cost $\theta_{\text{EC}}$, the coefficient on market size $\theta_{\text{RS}}$, and the competitive effect $\theta_{\text{EC}}$ parameters of the model.
Table (ref) reports the point estimates from the observed sample along with standard errors and 95% confidence intervals estimated using 250 cross-sectional bootstrap replications. The estimates all have the expected sign and are all significantly different from zero at the 5% level. Fixed costs are lowest for Sam's Club and highest for BJ's. Firms are more profitable in larger markets and entry by competitors reduces profits. The entry cost is large relative to fixed costs, as expected.
We note that for this application $k$-NPL yields very similar results to $k$-EPL. The estimated competitive effect is small, implying that $k$-NPL is likely to be stable. However, we could not know this a priori (recall that in our Monte Carlo Experiments, the $k$-NPL estimates of the competitive effect are biased towards zero). However, ex post, with the stable $k$-EPL estimates in hand, we can understand the performance of $k$-NPL this setting.
The industry under investigation is characterized by a significant number of monopoly markets. In the latest year of our sample, 2021, we observed a mere 14 triopoly markets (less than 1% of our sample). Only 7% of the markets had a duopoly, with two firms present, while 20% of the markets were monopolies, where a single firm operated. Interestingly, 72% of counties in our sample had no wholesale club stores at all in 2021. Our counterfactual exercise aims to explore the reasons behind this relative scarcity of duopoly and triopoly markets. Are strong competitive effects or high costs responsible?
To address this question, we conduct a counterfactual simulation in which we entirely eliminate the competitive effect in the model, allowing firms to operate as independent agents without considering their competitors' actions. If we observe a substantial increase in new entries, we may deduce that strong competition is deterring other firms from entering the market. Conversely, if we notice minimal change in entry behavior, it may suggest that competitive effects are insignificant, and costs are the primary factor driving firms' entry decisions.
To investigate this, we take the estimated structural parameters from the observed sample, set $\theta_{\text{RN}}=0$, and compute the counterfactual equilibrium. We solve the nonlinear system of equilibrium equations using the estimated equilibrium as the starting value. Subsequently, we compare the results of simulations conducted under this counterfactual scenario to those obtained using the estimated model parameters. We perform this procedure for each bootstrap replication to calculate standard errors for the counterfactual quantities of interest. For every simulation, we employ the observed market configuration at the beginning of the sample in 2009 and simulate until the end of our sample period in 2021.
In Table (ref), we present aggregate statistics from simulations using both the estimated and counterfactual parameters, where the competitive effect has been set to zero. The reported figures represent the estimated means and standard errors derived from simulated sample paths using the parameters from the 250 bootstrap replications. Before delving into the counterfactual analysis, it is worth noting that this can also serve as a test of the model's fit. We observe that the average number of active firms, entries, and exits, as well as the distribution of the number of firms present, align reasonably well with the observed values from the data. Upon examining the counterfactual, we discover that removing the effect of competition does not significantly reduce the number of unserved markets (from 1164 to 1156). However, it does shift the distribution away from monopoly markets (from 341 down to 300) towards a higher prevalence of duopoly and triopoly markets (from 93 up to 118, and from 11 up to 34, respectively). These findings indicate that competition in this industry is relatively weak, and costs play a more significant role in determining entry behavior.
Table (ref) provides a more detailed breakdown of our simulations by firm and market size state $s$. We estimate that, on average, Sam's Club would enter approximately 3 additional small markets ($s=1,2,3$) and around 15 additional large markets ($s=4,5$). Costco would enter, on average, 4 additional small markets and 28 large markets. BJ's would enter, on average, 3 additional small markets and 27 large markets. Consequently, eliminating the competitive effect appears to have a proportionally much larger impact on BJ's compared to Costco and Sam's Club. This effect is also considerably more pronounced in larger markets than in smaller ones.
Our analysis also incorporates average profit simulations using both the estimated parameters and the counterfactual parameters. To achieve this, for each simulation, we calculate $v^{j}(x,a^{j})+e(x,a^{j};P^{j})$ for each firm and average it over each simulated sample path of states and choices.\footnote{We simulate actions by sampling from the distribution defined by $P^{j}$, so $e(x,a^{j};P^{j})\equiv E[\varepsilon^{j}(a^{j})\mid x,a^{j},P^{j}]$ must be included to capture the contribution of $\varepsilon^{j}$ to average profit. With Type 1 Extreme Value shocks, $e(x,a^{j};P^{j})=-\ln P^{j}(x,a^{j})+\bar{\gamma}$, where $\bar{\gamma}$ is the Euler-Mascheroni constant.} We perform this separately for the simulations based on the model estimates and those under the counterfactual parameters. Comparing these two unit-less measures of profitability offers another perspective on our research question.
The results of this exercise are reported in Table (ref). The resulting changes in average profits are relatively minor: a 2.0% increase for Sam's Club, a 2.5% increase for Costco, and a 1.9% increase for BJ's. These results provide further insights regarding the relative significance of competition and costs in determining market entry behavior. The modest changes in average profits under the counterfactual scenario, where competitive effects are removed, suggest that competition does not play a dominant role in shaping firms' entry decisions. Instead, this reinforces our earlier findings that costs are a more substantial determinant of entry behavior in this industry.
We proposed an iterative $k$-step Efficient Pseudo-Likelihood ($k$-EPL) estimation sequence that extends the attractive econometric and computational properties of the single-agent $k$-NPL sequence to games. The nice econometric properties arise because $k$-EPL uses Newton-like steps on the fixed point constraint at each iteration. As a result, $k$-EPL is stable for all regular Markov perfect equilibria, each EPL iteration has the same limiting distribution as the MLE, and further iterations achieve higher-order equivalence and quickly converge to the finite-sample MLE almost surely. Computational advantages follow from defining the equilibrium conditions with choice-specific value functions, with standard modeling assumptions reducing each EPL iteration to two steps: (i) solving linear systems to generate pseudo-regressors, followed by (ii) solving a globally concave static logit/probit maximum likelihood problem using the pseudo-regressors.
In a real-world application, we use $k$-EPL to investigate the effect of competition on entry and exit of U.S. wholesale club stores. Our estimated model indicates that competition among wholesale club stores has a relatively mild effect on their entry and exit.Our Monte Carlo simulations show that $k$-EPL performs favorably in finite samples, is robust to data-generating processes where standard $k$-NPL encounters serious problems, and scales better than other iterative alternatives to $k$-NPL.
One limitation of our analysis is that we did not consider time-invariant unobserved heterogeneity in estimating dynamic discrete games. $k$-EPL can easily accommodate a proxy variable approach (e.g., CollardWexler2013), where an observed time-invariant variable is used to proxy for the the time-invariant unobserved heterogeneity. Without a proxy variable, it may also be possible to modify the $k$-EPL algorithm to incorporate time-invariant unobserved heterogeneity while preserving computational convenience and econometric efficiency. However, we leave such a substantial and challenging extension as an avenue for future research.
We thank J\'er\^ome Adda (editor) and three anonymous referees for numerous comments which greatly improved the paper. We also thank Dan Ackerberg, Victor Aguirregabiria, Lanier Benkard, Chris Conlon, Arvind Magesan, Mathieu Marcoux, Robert Miller, Salvador Navarro, John Rust, and Eduardo Souza-Rodrigues for comments and insightful discussions. The paper has benefited from comments by seminar participants at Boston College, Cornell University (Johnson), Georgetown University, University of Notre Dame, University of Pennsylvania (Wharton), University of Toronto, the 2019 International Industrial Organization Conference (IIOC, Boston), the 2019 University of Calgary Empirical Microeconomics Workshop (Banff), the 2020 ASSA Annual Meeting (San Diego), the Spring 2020 $\textrm{IO}^{2}$ Seminar, the 2020 Econometric Society World Congress (Milan), the 2021 North American Meeting of the Econometric Society, and the 2023 International Association for Applied Econometrics conference (Oslo).
The data and code underlying this research is available on Zenodo at \url{https://dx.doi.org/10.5281/zenodo.10582004}.