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.
78,614 characters · 16 sections · 40 citation commands
Efficient Estimation of Structural Models via Sieves
Keywords: Efficient Estimation, Sieves, Empirical Games, Joint Algorithm, Nested Algorithm
A structural model builds on economic theory and describes how a set of endogenous variables are related to a set of explanatory variables. This relation is often in the form of an implicit function. In particular, it generates endogenous function $p$ determined by an equation system:
where $\theta$ is the parameter of interest and $\Psi$ is a representation of the structural model.\footnote{In discrete choice models, the parameter captures consumer preferences and the observable is consumer choice; in auctions, the parameter captures the value distribution and the observable is the bid distribution; in dynamic models, the parameter describes the agent's intertemporal tradeoff and the observable is intertemporal choice.} While $\Psi$ is explicit, solving for $p$ could be difficult or costly. Such computational burden limits the use of standard estimators. For instance, the maximum likelihood estimator (MLE) repeatedly guesses $\theta$ and evaluates data likelihood using the solution of (ref), $p^*(\theta)$. However, finding the solution can be computationally intensive and is often hindered by a lack of robust algorithms, particularly in empirical games.
We introduce a new class of estimators: sieve-based efficient estimators for structural models (SEES). Our approach applies to a broad class of models, including empirical games, and avoids solving the model. Our SEES\ is motivated by two popular approaches to infinite-dimensional optimization problems: approximation and penalization. See shen1997methods, shen1998method, and chen2007large. Because the likelihood function $\ell(p^*,\theta; \text{data})$ involves an unknown function $p^*$, maximizing the likelihood with respect to $p^*$ and $\theta$ may lead to an asymptotically inefficient estimator for the parameter, and the resulting estimator may not necessarily be close to the solution of (ref). To address these issues, prior studies utilize sieves that are less complex but dense to approximate the original function space, and regularization that assumes smoothness of this function.
In this paper, we estimate structural models by approximating the solution to avoid solving the model and regularizing with the equilibrium conditions that are built into the model itself. By combining the data fitting and model fitting criteria, we formulate our penalized log-likelihood criterion, \[ \underbrace{\ell(\beta,\theta;\text{data})}_{\text{data likelihood}} - \underbrace{\omega \times \rho(\beta,\theta)}_{\text{penalization}}, \] where $\ell$ and $\rho$ measure the data fitting and model fitting, respectively.\footnote{While our idea extends to other types of estimators, we focus on likelihood-based ones here.} Moreover, $\beta \in \mathbb{R}^K$ and $\omega \in \mathbb{R}_+$ govern the approximation and the weighting, respectively. Instead of imposing stronger smoothness assumptions than typically implied by theory, our approach relies solely on the model to regularize the sieve approximation. The smoothing parameter $\omega$ explicitly captures the weighting of the data likelihood and the equilibrium condition, and the dimension of the approximation parameter $\beta$, denoted by $K$, balances computational cost and solution accuracy.
Allowing these tuning parameters to diverge at appropriate rates, the proposed parameter estimator of $\theta$ is consistent, asymptotically normal, and asymptotically efficient. Intuitively, by gradually updating the smoothing parameter, we shift the weight from the data to the equilibrium condition. At the minimum, a preliminary nonparametric estimate of $p$ (by letting $\omega =0$) constitutes a good starting value but is subject to issues with nonparametric estimates. When the smoothing parameter increases, more weight is given to the equilibrium condition. By forcing model restrictions more strongly, the estimator converges to the MLE.
We prescribe several algorithms to implement SEES. The first is a joint algorithm that finds the combination of sieve approximations and model parameters that best explains the data and satisfies the equilibrium conditions. That is, we maximize the penalized log-likelihood function with respect to $(\beta,\theta)$. The second is a nested algorithm that consists of two main parts. First, for each model parameter $\theta$, we find the sieve approximation of $p$ that best explains the data and satisfies the equilibrium conditions. Second, based on the approximated solution, we find the model parameter that best fits the data. While the joint algorithm is attractive because it results in a single-level optimization problem, the nested algorithm is quite intuitive, resembling MLE.
Our estimator allows for discrete and continuous state/heterogeneity in the model to be estimated. The standard practice of discretizing continuous state variables or covariates leads to efficiency loss. Under mild regularity conditions, we show that our estimator has the same asymptotic distribution as MLE in both cases. To our knowledge, we are the first to combine approximation and penalization in estimating structural models. While some studies have adopted approximation approaches, none combines them with penalization. Another important advantage of our method is that it produces standard errors in the same way as the standard MLE using the Fisher information matrix, which is of considerable convenience in empirical work. As a side product, we also derive a similar approach for the mathematical program with equilibrium constraints (MPEC) estimators, which provides a faster alternative than the bootstrapping method previously proposed by su2012constrained.
We acknowledge several limitations inherent in our methodology. First, our sieve-based approach presupposes the smoothness of the solution within the state variables or covariates, leaving the treatment of discontinuities as a subject for future investigation. Second, our approach provides a robust solution that works with minimal assumptions on the solution, which is particularly valuable in models with unfavorable or unknown properties. However, it may not always be the most expedient choice in scenarios where the solution exhibits favorable properties, such as contraction mappings. For instance, as demonstrated in the simple model outlined in Section (ref), it exhibits a relatively slower performance when compared to a nested-fixed point algorithm. Throughout this paper, we refrain from comparing computational time across different estimators, as it is often model-specific and, hence, more relevant in richer empirical models.
The remainder of the paper is organized as follows. Section (ref) explains the idea using a simple example. Section (ref) proposes the class of sieve-based efficient estimators for structural models and derives its asymptotic properties. Section (ref) demonstrates the performance of our estimators in estimating an empirical game. Section (ref) concludes. The Appendix contains all omitted proofs and details.
Our new method differs from existing methods by how we leverage data and model restrictions. We now compare it with popular existing methods, such as maximum likelihood estimation (MLE), two-step approaches, and nested pseudo-likelihood (NPL), through a motivating example.
Consider a monopolist $j$, facing logit demand, sells a product at a price $P_j$. That is, consumer $i$ gets utility of \[ u_{ij} = \xi_j - \alpha P_j + \epsilon_{ij} , \] where $\xi_j$ is continuous product quality, $\alpha$ is the price coefficient, and $\epsilon_{ij}$ represents the standard Type 1 extreme value (T1EV) taste shock. The firm's profit maximization problem is \[ \max_{P_j} \quad (\underbrace{P_j - c_j}_{\text{profit margin}}) \times \underbrace{\frac{\exp(\xi_j - \alpha P_j)}{1+\exp(\xi_j - \alpha P_j)}}_{\text{market share}} , \] where $c_j$ represents the constant marginal cost. The optimal price is determined by the FOC, \[ \alpha(P_j - c_j) = 1+\exp(\xi_j - \alpha P_j), \] where $P_j$ appears both inside and outside an exponential function. As a result, the mapping from the parameters to the optimal price is implicit.
For simplicity, we focus on estimating the parameter $\theta$ that governs consumer preferences over product feature $x_j \in \mathbb{R}$ using observed prices. Specifically, we treat it as known that $c_j = 0,~\alpha = 1, ~\xi_j = \log x_j + \log \theta + 1$, where “1” is quality normalization for simplicity. Appendix (ref) shows that the optimal price satisfies
where $y_j^* = P_j^* -1 $ represents the normalized price and $p(x_j; \theta)$ is defined by \[ p(x_j; \theta) e^{p(x_j; \theta)} = \theta x_j \quad \text{ or } \quad p(x_j; \theta) = \theta x_j e^{-p(x_j; \theta)} , \] the first of which has the standard form of the Lambert W function\footnote{The Lambert W function $W(x)$ is defined by $W(x) e^{W(x)} = x$.} and the second of which has the same form as Equation (ref). We denote this function as $p(\cdot; \theta)$ to indicate its dependence on the parameter.
Consider a data generating process (DGP) that is a noisy measurement of the optimal price $y_j = y_j^* + e_j$, where $e_j$'s are measurement errors that are i.i.d. draws from the standard normal distribution. Therefore, the observed (normalized) price $y_j$ is from the standard normal distribution with a location $p(x_j;\theta_0)$: \[ y_j \sim \mathcal{N}(p(x_j;\theta_0), 1), \text{ where } j=1,\ldots, N . \] The data contain the product characteristics $x_j$ and the prices $P_j= P_j^* + e_j$ (equivalently, the normalized prices $y_j$). The parameter of interest is $\theta$.
Maximum Likelihood Estimation: The standard MLE solves the following problem \[ \max_{\theta} \quad \sum_{j=1}^N \log \phi(y_j - p(x_j;\theta)), \] where $\phi(\cdot)$ represents the density function of the standard normal distribution. Because $p(\cdot;\theta)$ is implicitly defined, this estimator is computationally costly. For each trial of $\theta$, we need to find $p(x_j;\theta)$ for each data point $x_j$. The number of equations that need to be solved equals the sample size multiplied by the number of likelihood function evaluations.
Despite its asymptotic efficiency, the standard MLE requires solving the model for each parameter and thus solution algorithms that are sufficiently efficient and robust. When a contraction mapping solution for the model is available, it is often referred to as the nested fixed point algorithm (NFXP). See, e.g., rust1987optimal in dynamic discrete choice models and berry1995automobile in demand models.
A Two-Step Approach: We can “invert the FOC” and obtain a representation of the “unknown” in terms of the optimal prices: \[ \theta = \frac{p(x_j; \theta) e^{p(x_j; \theta)} }{x} , \] where the normalized price $y_j^*= p(x_j; \theta)$ is unobserved. In principle, this FOC inversion allows estimating the parameter using the optimal price in any market.
Due to measurement errors, the rewritten FOC suggests a simple two-step approach that avoids solving the model repeatedly in estimation. In the first step, we consider $y_j^* = p(x_j; \theta_0)$ and estimate the optimal price as a function of the covariate. Although the true parameter $\theta_0$, the endogenous variable $p$ and thus the RHS $p(x_j; \theta_0)$ are all unobserved, we can estimate the LHS $y_j^*(x_j)$ nonparametrically using the observed price and covariate pairs $\{x_j, y_j\}_{j=1}^N$. In particular, we run a nonparametric regression\footnote{We apply kernel regression using the optimal bandwidth estimated by cross-validation, as provided on Professor Yingying Dong's webpage: http://yingyingdong.com.}, \[ y_j = y_j^*(x_j) +e , \] and obtain an estimate of the normalized price $\widehat{y_j^*(x_j)}$. In the second step, we have a simple plug-in estimator, \[ \widehat{\theta } = \text{median}\left\{\frac{\widehat{y_j^*(x_j)} \times \exp{\widehat{y_j^*(x_j)} } }{x_j}: j = 1, \ldots, N\right\}. \]
Two-step approaches avoid repeatedly solving the economic model at the expense of efficiency. In the first step, the analyst obtains a nonparametric estimate of the endogenous variable $\widehat{p}$. In the second step, the estimate is obtained from $\widehat{p} = \Psi(\widehat{p}, \theta)$ in various ways. In auction models, guerre2000optimal use the estimated bid distribution to construct pseudo values, which are then used to estimate the underlying value distribution. In dynamic discrete choice models, the conditional choice probability (CCP) approach of hotz1993conditional plugs the estimated CCPs into the optimal decision rules. In dynamic games, one can obtain a nonlinear least squares estimate of $\theta$ by replacing $p$ with the estimated CCP $\widehat{p}$ in the function; see pesendorfer2008asymptotic.
Nested Pseudo-Likelihood Algorithm: In each iteration, the NPL algorithm solves the following problem: \[ \max_{\theta} \quad \sum_{j=1}^N \log \phi(y_j -\theta x_j \exp(-\widehat{p}_j^\tau)), \] where $\widehat{p}_j^\tau$ represents some estimate of the optimal price in market $j$. Denote the solution as $\widehat{\theta}^\tau$. We can then update the price estimates $\widehat{p}_j^{\tau+1} = \widehat{\theta}^\tau x_j \exp(-\widehat{p}_j^\tau)$. We iterate the process till the parameter estimate converges.
Given some estimates $\widehat{\theta}$ and $\widehat{p}$, the NPL algorithm obtains new estimates of the choice probabilities by applying the mapping $\tilde{p} = \Psi(\widehat{p},\widehat{\theta})$ and then updates the parameter estimate by maximizing the pseudo-likelihood function $\ell(\tilde{p}, \theta)$.
A Sieve-Based Efficient Estimator: In this paper, we propose a method that obviates solving the model repeatedly. In particular, we approximate the “solution” function $p^{\star}(\cdot)$ by B-spline basis functions: \[ p^\beta(\cdot) = \sum_{k=1}^K \beta_k s_k(\cdot), \] where $s_k(\cdot)$ is a cubic spline basis function, and $K$ denotes the number of basis functions. Our sieve-based estimator of $\theta$ maximizes the likelihood \[ \sum_{j=1}^N \log \left\{\phi\left(y_j - p^{\widehat{\beta}(\theta,\omega)}(x_j)\right)\right\} , \] where $\widehat{\beta}(\theta,\omega)$ is defined by \[ \arg\max_{\beta \in \mathbb{R}^K} \quad \sum_{j=1}^N \log \left\{\phi(y_j - p^\beta(x_j)) \right\} - \underbrace{\omega \int_{\mathcal{X}} \left[p^\beta(x) e^{p^\beta(x)} - \theta x\right]^2 dx}_{\text{penalization}}, \] where $\omega >0$. Because $p^\beta(\cdot)$ is an approximation, the second term penalizes the likelihood by the amount of deviation by definition of the Lambert W function.
We now compare the above-mentioned estimators. First, SEES\ and MLE are asymptotically equivalent and almost identical in finite samples. However, NFXP algorithms may converge slowly, and such mappings may not even exist in important models. For instance, empirical games, such as asymmetric auctions and dynamic games, are notoriously difficult to solve, making MLE difficult to apply. In contrast, we avoid solving the model repeatedly by approximating the solution flexibly.
Second, two-step approaches are limited by the first-step nonparametric estimation of the endogenous variable and may suffer from the “curse of dimensionality” when $x$ has multiple dimensions stone1980optimal. As a result, the finite-sample estimation error can be substantial. In contrast, our approximated solution avoids this issue, as its final version relies almost entirely on the model.
Third, our estimator is also related to the nested pseudo-likelihood algorithm proposed by aguirregabiria2002swapping,aguirregabiria2007sequential. Exploiting the unique feature of dynamic discrete choice models that the Jacobian matrix of $\Psi_\theta$ is always zero, their iterative refinement converges to MLE. However, it requires some discretion in applying it to empirical games. See pesendorfer2010sequential. Both algorithms bridge the gap between the standard MLE and two-step methods, and are asymptotically equivalent to MLE. However, they are based on very different ideas. Our estimator is flexible to accommodate different estimation algorithms, including one that resembles NPL, and robust in applications to various models, including empirical games.
To our knowledge, we are the first to combine approximation and penalization in estimating structural models. While some studies have adopted approximation approaches, none combines them with penalization. For instance, keane1994solution,keane1997career use sieves to approximate solutions in dynamic structural models, combining approximation and NFXP. In estimating dynamic games, sweeting2013dynamic uses parametric approximations to the value function, combining approximation and NPL. Most related, barwick2015costs approximates the value function using sieves and imposes the Bellman equation as an equilibrium constraint.
Another related algorithm is MPEC, which is an alternative computational algorithm to MLE. See, e.g., su2012constrained and dube2012improving.\footnote{Several papers have compared MPEC with the original estimators for various models. For a comparison of NFXP and MPEC, see, e.g., lee2015computationally for the berry1995automobile model and iskhakov2016comment for dynamic models.} It avoids solving the model repeatedly by augmenting the unknown to $(\theta,p)$ and imposing the equilibrium condition as a constraint:
We will show that our SEES's dual problem is a natural extension of MPEC in discrete state settings. Our SEES\ nests MPEC as a limiting case when the number of basis functions is the same as the dimension of $p$ and the regularization parameter equals infinity. MPEC forms the Lagrangian function using Lagrange multipliers $\Lambda$ that are of the same size as $p$: $\max_{\theta,p,\Lambda} \ell(p,\theta;\text{data}) - \Lambda^\prime (p - \Psi(p,\theta))$. As a result, it solves $2 \times \text{dim}(p)+\text{dim}(\theta)$ equations in the same number of unknowns. Our SEES\ approximates $p$ by $\beta$ and introduces a scalar regularization parameter $\omega$, which reduces the problem to an unconstrained optimization problem with $K +1+\text{dim}(\theta) $ unknowns. Therefore, SEES\ is computationally convenient because $K+1 \ll 2 \times \text{dim}(p)$.
Consider $x_j \sim \text{Uniform}[0,\overline x]$ and $\theta_0 = 1$. For MLE, we use the bisection method to solve for the Lambert W function. We omit MPEC here for two reasons. First, we focus on the statistical properties rather than the computational time of the estimators. Second, MPEC is an alternative computational algorithm to MLE with identical statistical properties. For the two-step approach, we use local linear kernel regression and apply the optimal bandwidth chosen via cross validation. For the proposed method, we use the cubic spline explained in luo2018structural and let $K=6$; the choice of $\omega$ follows the method that we propose later. We also provide the analytic gradient of the outer loop and the analytic gradient and Hessian of the inner loop maximization problem; see Appendix (ref).
Table (ref) shows the simulation results of 1,000 replications with a sample size of 1,000. SEES, MLE, and NPL perform very similarly. In particular, SEES\ and MLE are almost identical in each replication. Figure (ref) compares MLE and the iterations of NPL and SEES\ in a typical replication.\footnote{SEES\ usually converges in 2-4 iterations using our proposed choice of $\omega$. For better visualization in this figure, we increase $\log \omega$ in 7 equal steps to match the number of iterations of NPL.} While their earlier iterations could differ from MLE significantly, NPL and SEES\ both converge to a close neighbourhood of MLE.
In contrast, the two-step approach generates a larger bias and standard error. Alternatively, we can take the sample average in the second step. However, the noise in nonparametric estimates near boundaries deteriorates the estimates substantially. The median performs much better the mean. It is clear that the performance of the two-step estimator is affected by the first-stage nonparametric regression.
Now, we describe our estimator in detail. Instead of solving for $p^*(\theta)$ for each parameter $\theta$ in the likelihood evaluation, we approximate the true solution by $p^\beta$. The choice of approximation infrastructure depends on its approximation properties and computational convenience. A popular one often adopted in empirical studies is the method of sieves; that is, $p^\beta(\cdot)=\sum_{k = 1}^K \beta_{k}s_{k}(\cdot)$, where $\{s_1, \ldots, s_K\}$ represent the basis functions of the finite-dimensional sieve space $\mathcal{B}$. Typical choices of $s_k$ include B-spline basis functions and Bernstein polynomials. Such methods are flexible in accounting for shape restrictions imposed by the structural model, such as nonnegativity and monotonicity. For instance, if $p$ represents choice probabilities, we can use $p^\beta(\cdot)=[1+\exp\{\sum_{k = 1}^K \beta_{k}s_{k}(\cdot)\}]^{-1}$ to ensure that $p\in (0,1)$.
Moreover, our SEES\ imposes the model constraints by penalizing the difference between $p^\beta$ and $\Psi(p^\beta,\theta)$. This difference is independent of the data sample in measuring the fidelity of approximation to the equilibrium conditions. The smaller the difference, the better the approximate $p^\beta$ satisfies equilibrium conditions.
We formulate the penalized log-likelihood criterion by combining the data fitting and model fitting criteria,
where $\omega >0$ is the smoothing parameter, and $\rho$ is a metric that measures the difference between $p^\beta$ and $\Psi(p^\beta,\theta)$. For instance, $\rho(p,\Psi(p,\theta)) = \|p - \Psi(p, \theta)\|_2^2$, where $\|\cdot\|_2$ is the Euclidean norm. For simplicity, we shall use as shorthand $\ell(\beta,\theta)$ and $\rho(\beta,\theta)$.
We develop two algorithms to implement our estimator given each smoothing parameter $\omega$: a joint algorithm and a nested algorithm.
Joint Algorithm This algorithm is attractive because it involves a single-level optimization problem. We augment the unknown to $(\beta,\theta)$ and solve the following problem:
which leads to $\widehat{\beta}(\omega)$ and $\widehat{\theta}(\omega)$. We recommend supplying an analytic gradient and Hessian to reduce computational cost and to increase precision.
Nested Algorithm This algorithm is intuitive, resembling MLE. There are two layers of optimization problems to be solved. In the inner layer, given $(\theta, \omega)$, we find the best approximation parameter $\widehat{\beta}(\theta; \omega)$ that solves the following problem:
Solving (ref) indicates that the maximizer $\widehat{\beta}(\theta; \omega)$ is an implicit function of $\theta$.
In the outer layer, applying the best fitting approximation parameter, we search for the structural parameter $\widehat{\theta}(\omega)$ that maximizes the following likelihood:
Note that we have considered the structural equation (ref) in the inner layer. Therefore, the equilibrium conditions are embedded in $\widehat{\beta}(\theta, \omega)$. As illustrated above, the optimizer of (ref), $\widehat{\beta}(\theta; \omega)$, is an implicit function of $\theta$. Therefore, $\ell \left(\widehat{\beta}(\theta,\omega),\theta \right) $ in (ref) is a function of $\theta$. We obtain the final estimator of $\theta$, denoted by $\widehat{\theta}(\omega)$, by directly maximizing (ref) with respect to $\theta$.
We propose a new method that selects the smoothing parameter $\omega$ based on the performance of parameter estimation. Intuitively, we choose a sufficiently large $\omega$ to ensure fidelity in approximating the equilibrium conditions.\footnote{Cross-validation is commonly used to balance bias and variance in estimation or prediction when the function of interest is unknown, by estimating prediction error or evaluating the likelihood function on held-out data. However, in our case, this trade-off is not a concern because the function $p(\cdot)$ is fully specified by the structural model (ref).} In particular, we start with a moderate $\omega_1$ and update it till the estimates converge. For each $\omega = \omega_\tau$, we can conduct the joint or nested algorithm and obtain an estimate $\widehat{\theta}(\omega_{\tau})$. Consider a significance level of $\alpha$. We obtain its standard error $\widehat{\sigma}(\omega_{\tau})$ and confidence interval $\mathcal{I}(\omega_{\tau}) = [\widehat{\theta}(\omega_{\tau}) +z_{\alpha/2} \widehat{\sigma}(\omega_{\tau}), \widehat{\theta}(\omega_{\tau}) +z_{1-\alpha/2} \widehat{\sigma}(\omega_{\tau})]$, applying the standard formula of standard error calculation for MLE.
We multiply the smoothing parameter by $L$ each time, i.e. $\omega_{\tau+1} = L \times \omega_{\tau}$ and obtain a new estimate and its confidence interval $\mathcal{I}(\omega_{\tau+1})$. Continue this process till the overlapping portion of the two intervals accounts for more than a threshold percentage of both of the two. That is, our final choice of the smoothing parameter is $\widehat{\omega} = \omega_\tau$ if \[ \min \left\{ \frac{|\mathcal{I}(\omega_{\tau}) \cap \mathcal{I}(\omega_{\tau-1})|}{|\mathcal{I}(\omega_{\tau-1})|}, \frac{|\mathcal{I}(\omega_{\tau}) \cap \mathcal{I}(\omega_{\tau-1})|}{|\mathcal{I}(\omega_{\tau})|} \right\} \geq \mathfrak{c} , \] where $|\cdot|$ represents the length of the interval. In case the parameter of interest $\theta$ is multi-dimensional, we check this condition element by element. The final estimate of the model parameter follows $\widehat{\theta}(\widehat{\omega})$.
All the simulations reported in this paper adopt $\alpha = 0.05$, $L = 10$, and $\mathfrak{c} = 95\%$. Therefore, $z_{\alpha/2} = -1.96$ and $z_{1-\alpha/2} = 1.96$. To illustrate how the proposed method works in terms of selecting the tuning parameter, consider the monopoly pricing example with $x_j = 1$. Figure (ref) reports the parameter estimate and confidence intervals when the smoothing parameter varies. The x-axis represents $\log \omega$. While the bias seems small for small values of $\omega$, the confidence interval is large. As $\omega$ increases, it shrinks to the MLE confidence interval.
The approximation parameter $K$ (i.e., the number of basis functions) should be chosen properly: large enough to approximate the equilibrium well. Exactly how many is sufficient depends on the complexity of the equilibrium solution. In our motivating example, the solution is simple; as a result, we find that four cubic basis functions are adequate to approximate the solution well. When the solution is complex and the number of basis functions needs to be large, the analyst should start with a large smoothing parameter to avoid over-fitting in the inner loop; i.e., the likelihood function dominates. Sometimes, there are a finite number of states in the structural model, which is often assumed in estimating dynamic models. See, e.g., aguirregabiria2002swapping and pesendorfer2008asymptotic. Such finite states often come from discretization of covariates. In this case, the approximation can be perfect, i.e., $\beta = p$. That is, $s_k(x) = \mathbbm{1}(x = p_k)$, where $\mathbbm{1}(\cdot)$ is the indicator function and $p_k$ is the $k$th element in the endogenous variable $p$. Note in this scenario the number of basis functions is identical to the dimensionality of $p$.
To capture empirically relevant covariates without losing much efficiency, any approximation methods would suffer from a computational curse of dimensionality --- the total number of basis functions has to grow fast as the dimensionality of $x$ increases. We propose to resolve this issue in several ways. First, more advanced approximation methods are often preferable to simple ones. See, e.g., chen2021efficient compare neural networks-based estimators. Additional shape constraints, sparsity patterns, and better grid choices are useful in reducing computational burden. See, e.g., chen2007large discusses various sieve-based methods, and kristensen2021solving discuss various approximation architectures for approximating value functions in dynamic models. Second, there are also many model-specific techniques for approximating functions using a small number of basis functions. The model may generate multiple endogenous objects, some directly observable while others intermediate. A well-chosen $p$ simplifies its approximation and evaluating $\Psi$ and data likelihood. For instance, in static games of asymmetric information, if the deterministic component in the payoff function is linear in the parameters, see, e.g., bresnahan1991empirical and bajari2010estimating, how covariates determine the endogenous variable becomes a multiple-index model. In empirical auctions, chen2023identification approximate the bid-stage primitives by flexible Bernstein polynomial sieves. One can borrow techniques from the existing literature in estimating such a model.
As mentioned above, our proposed method can handle both continuous states, which result in an infinite-dimensional endogenous variable $p$, and discrete states, which lead to a finite-dimensional $p$. We begin by examining the discrete state settings, as this will provide insights into the continuous state cases.
Let $p = (p_1, \ldots, p_{d_1})^{\prime} \in \mathbb{R}^{d_1}$ be the endogenous variable in (ref). Without loss of generality, we assume that $\theta \in \mathbb{R}^{d}$ and $\Theta$ denotes the space of $\theta$. Under certain conditions that are specified below, we establish consistency and asymptotic normality of the joint estimator in the main text. The nested estimator is also consistent and asymptotically normal; see Theorems (ref) and (ref) in Appendix (ref). In fact, the two estimators have the same asymptotic distribution.
For any given $\theta \in \Theta$, we aim to maximize the following function with respect to $p$
where $\ell_n$ denotes the log-likelihood corresponding to $n$ i.i.d. observations and $\|\cdot\|_2$ is the Euclidean norm of a vector. Suppose $Y_1, \ldots, Y_n$ are i.i.d. observations, taking values in $\mathbb{R}^{d_2}$. We assume that the likelihood function in (ref) can be written as $$ \ell_n(p) = \frac{1}{n} \sum_{i = 1}^n f(Y_i, p), $$ where $f$ is a function defined on $\mathbb{R}^{d_2} \times \mathbb{R}^{d_1}$. For simplicity, we assume that $d_2 = 1$. In this context, $f(y, p)$ is actually the log density function of $Y_i$.
By Brouwer's fixed-point theorem, there must exist a solution to (ref) for any $\theta \in \Theta$. For instance, $p\in [0,1]$ represents CCP in dynamic games. Define $g(p, \theta) = \Psi(p, \theta) - p$. Obviously, the solution to Equation (ref), denoted as $p^{*}(\theta)$, satisfies $g(p^{*}(\theta), \theta) = 0$. We impose the following regularity condition on $g(p, \theta)$.
Assumption (ref) can be understood as a local inverse Lipschitz condition. Consider the Lambert function $p = W(\theta)$, which is defined implicitly by $pe^{p} = \theta$. In correspondence, $g(p, \theta) =\theta e^{-p} - p$. Thus, $(p + g)e^{p + g} = \theta e^{g}$, which implies that $p + g = W(\theta e^{g})$ or $$ p(g) = -g + W\left(\theta e^{g}\right). $$ Since $W$ is a continuously differentiable function by the implicit function theorem, for any $g_1, g_2$ satisfying $|g_1| \vee |g_2| \leq r$ with some constant $r$, we have $|p(g_1) - p(g_2)| \leq C|g_1 - g_2|$ for some constant $C$.
By the implicit function theorem, this assumption ensures that the solution to Equation (ref), $p = p^{*}(\theta)$, is a continuously differentiable function of $\theta$. Let $\theta_0$ denote the true value of $\theta$. Define $$ M(\theta) = {\rm E}\,_{\theta_0}[f(Y_i, p^{*}(\theta))] \quad\quad\text{for}~ \theta \in \Theta, $$ where the expectation is taken with respect to $\mathbb{P}_{\theta_0}$.
This assumption ensures that the true parameter $\theta_0$ is identifiable.
This assumption holds for square loss functions, i.e., $f(y, p) = -(y - p)^2$. Given any $\theta \in \Theta$ and a positive $\omega$, recall the sieve estimate of $p$ is given by
The following theorem indicates the approximate solution to the structural equation (ref) is uniformly close to the exact solution $p^{*}(\theta)$.
To establish the consistency of the estimator $\hat{\theta}_n$, we need a stronger version of Assumption (ref).
Let $\tilde{\theta}_n$ denote the estimator of $\theta$ obtained from the joint algorithm. Actually, $\tilde{\theta}_n$ is defined by
where $\hat{p}(\theta)$ is given by (ref).
To derive asymptotic normality, we need a stronger condition than Assumptions (ref) and (ref).
The joint algorithm is attractive because it involves a single-level optimization problem and computes the Hessian matrix with respect to $(\beta,\theta)$ at the solution directly. The following corollary provides a natural way to calculate the standard error of $\tilde{\theta}_n$ using the Hessian matrix generated from the joint algorithm.
In this section, we consider continuous states in the structural model and examine the large-sample properties of SEES.
Let $x$ denote the continuous states. Without loss of generality, we assume that the dimension of $x$ is 1. Given any $\theta \in \Theta \subset \mathbb{R}^d$, $p^*(x; \theta) \in \mathbb{R}^{d_1}$ denotes the solution to the following the structure equation:
where $x \in [0, T]$ with a fixed $T > 0$. Without loss of generality, we assume $T = 1$. We approximate the solution to (ref) using a sieve method. In particular, we take the sieve space, denoted as $\mathcal{B}_n$, to be the space of cubic B-spline functions equipped with knots $\tau^{(n)} = \left\{0 = t_1^{(n)} < \cdots < t_{M_n}^{(n)} = T\right\}$. Let $|\tau^{(n)}| = \max_{1 \leq i \leq M_n - 1} |t_{i + 1}^{(n)} - t_{i}^{(n)}|$ be the largest distance of adjacent knots in $\tau^{(n)}$. For any element $\eta \in \mathcal{B}_n$, there exists a $\beta = (\beta_1, \ldots, \beta_{K_n})^{\top} \in \mathbb{R}^{K_n}$ such that $\eta(x) = \sum_{j = 1}^{K_n} \beta_j s_j(x)$, where $s_j$'s are cubic B-spline basis functions and $K_n = M_n + 3$.
Suppose that $(Y_1, X_1), \ldots, (Y_n, X_n)$ are i.i.d observations, where $X_i$'s are independently sampled from a distribution $Q$ on $[0, T]$ and $Y_i$'s take values in $\mathbb{R}^{d_2}$. For simplicity, we assume that $d_1 = d_2 = 1$. Following the notations defined in the last subsection, given $\theta \in \Theta$, we have
where the likelihood $\ell_n$ can be written $$ \ell_n(p(\cdot)) = \frac{1}{n} \sum_{i = 1}^n f(Y_i, p(X_i)), $$ for any function $p$ defined on $[0, T]$, and the penalty function $\rho$ is given by $$ \rho[p(\cdot), \Psi(p(\cdot), \theta)] = \int_0^T \{p(x) - \Psi(p(x), \theta)\}^2 \mathop{}\!\mathrm{d} x. $$ Then the nested estimator of $\theta$ is given by
We next study the asymptotic properties of $\hat{\theta}_n$ in this model. Similar to Section (ref), we first establish consistency for the estimator and then develop the asymptotic normality of $\hat{\theta}_n$. To this end, we need the following regularity conditions.
Assumption (ref) ensures that $X_i$'s are evenly distributed over $[0, T]$, which is entailed by a good estimation of $p$ over the entire domain. Moreover, this assumption is commonly adopted in the literature of nonparametric smoothing; see stone1985additive and chen2023identification for example.
Define $g(p(x), \theta) = \Psi(p(x), \theta) - p(t)$ for any function $p$ defined on $[0, T]$. Let $h^{(k)}$ denote the $k$th order derivative of function $h$ for any integer $k \geq 0$, and define $$ C^k([0, T]) = \{h: h^{(k)}~\text{is continuous on}~[0, T] \}. $$ The following two assumptions extend their counterparts in the parametric model to the semi-parametric one.
By Lemma (ref) in Appendix (ref), we define the sieve space to be $$ \mathcal{B}_n(r) = \left\{\eta(x): \eta(x) = \sum_{j = 1}^{K_n} \beta_j s_j(x), \|\eta\|_{\infty} \leq r \right\}. $$ with equally spaced knots for some sufficiently large constant $r$. Therefore, for any $\theta \in \Theta$, there exits an $p_{\theta, n} \in \mathcal{B}_n(r)$ such that $\|p^*(\cdot; \theta) - p_{\theta, n}\|_{\infty} = O(K_n^{-4})$. Let
Similar to Section (ref), for the sieve estimator $\hat{p}(\cdot; \theta)$ defined in (ref), we establish an important approximation error bound, which will be used to develop the asymptotic normality for $\hat{\theta}_n$ later.
Let $\theta_0$ denote the true value of $\theta$ and $G_{\theta_0}$ denote the joint distribution of $(X_i, Y_i)$ under this true value. Define $$ M(\theta) = {\rm E}\,_{\theta_0}[f(Y_i, p^*(X_i; \theta)], \quad\quad \theta \in \Theta, $$ where the expectation is taken with respect to $G_{\theta_0}$. The follow assumptions are essentially the same as those in Section (ref).
To establish consistency and asymptotic normality for $\hat{\theta}_n$, a stronger version of Assumption (ref) is entailed.
MPEC: When the state space is discrete and $p$ is finite, our method could incorporate each element of the endogenous variable $p$ as a basis function in sieve approximation and put all of the weight on the equilibrium conditions in the inner loop. In this case, our estimator becomes the MPEC estimator.
Moreover, it is easy to show that for a given $\omega >0$, there exists an $\epsilon>0$ such that the optimization problem of the joint algorithm (ref) has the same solution as its dual optimization problem
The dual problem is a natural generalization of the MPEC estimator. However, solving it numerically is challenging.
Our above-mentioned derivation also suggests a natural way to calculate standard errors for MPEC estimators.
To the best of our knowledge, this result is new in the literature. su2012constrained suggest obtaining standard errors through bootstrapping. We derive the general result in Section (ref). Here, we consider $p,\beta,\theta \in \mathbb{R}$, as in the simple example with $x_j=1$, to explain the idea. When $\omega = \infty$ and $\beta = p$, our estimator is effectively an MPEC estimator, \[ \max_{g(\beta,\theta) =0} \quad \ell(\beta, \theta). \] The MPEC approach forms the Lagrangian function $h(\beta,\theta,\omega) = \ell (\beta, \theta ) + \lambda g(\beta, \theta)$. Note that this multiplier $\lambda$ should not to be confused with the smoothing parameter $\omega$ for general PSE. By definition, we have $g(\widehat{\beta}(\theta),\theta ) = 0$. Its first-order and second-order derivatives are
which allow for expressing $\widehat{\beta}^{\prime}(\theta)$ and $\widehat{\beta}^{\prime\prime}(\theta)$ in the gradient of $g$.
On the other hand, the second-order derivative of the likelihood is
where $\lambda$ denotes the associated Lagrange multiplier reported by a constrained maximization algorithm. The last equation follows from the Lagrange multiplier theorem that $\ell_\beta + \lambda g\lo\beta =0$ at the optimum $(\beta = \widehat{\beta}(\widehat{\theta}) , \theta = \widehat{\theta})$ and the first-order and second-order derivatives of the equilibrium constraints. All terms on the RHS are readily available if MPEC converges. We recommend supplying the analytic gradient and Hessian, as the numerical one can be inaccurate.
To the best of our knowledge, the theoretical properties of the MPEC estimator for structural models with continuous states remain unexplored. In contrast, our proposed method offers a rigorous framework for conducting statistical inference for $\theta$ with either discrete, continuous, or both types of states. We believe this represents a critical advancement for practical applications.
Approximate MLE: We now discuss the extreme case when we let $\omega = \infty$. That is, for each guess of the model parameter $\theta$, we find the best approximation to minimize any deviation from the equilibrium condition and then evaluate the likelihood by plugging in this best approximation. Specifically, our estimator becomes equivalent to
which looks similar to MLE, with an important difference that we only search for the best approximation in the inner loop. We call this special case of our estimator the approximate MLE (AMLE). Such approximate solution approaches have appeared in the literature. See, e.g., keane1994solution,keane1997career use sieves to approximate solutions in dynamic structural models.
One may wonder about the advantages of gradually changing $\omega$ instead of directly considering the limiting case. AMLE ignores the data when finding the best approximation of the solution for each $\theta$. Because the data are informative about the true strategies $p$, our general sieve-based efficient estimator may perform better than AMLE. By gradually updating the smoothing parameter, we shift the weight from the data to the equilibrium condition. At the minimum, a preliminary nonparametric estimate of $\widehat{p}$ (by letting $\omega =0$) constitutes a good starting value for the inner loop but is subject to issues with nonparametric estimates. When the smoothing parameter increases, more weight is given to the equilibrium condition. By forcing model restrictions more strongly, the estimates converge to MLE estimates.
In this section, we apply our methodology to an entry game between Walmart and Kmart, using a dataset published by jia2008happens. A detailed description of the industry and data is available in the original paper.
The original dataset includes 2,065 markets, each representing a county with an average population ranging from 5,000 to 64,000, covering the years 1988 to 1997. For our analysis, we focus on the year 1997. The market-level variables include the log of county population (pop), the log of retail sales per capita (spc), and the percentage of urban population (urban). Walmart-specific variables include an intercept, the log of distance to Bentonville (dbenton), and an indicator for the southern region. Kmart-specific variables include an intercept and an indicator for the Midwest region (midwest). These variables capture key variations in the data. For instance, a simple scatter plot of the total number of firms shows that neither firm enters the market when SPC is too low.
Denote the data as $\{d_{\mathcal{W} m}, d_{\mathcal{K} m}, x_{\mathcal{W} m}, x_{\mathcal{K} m}, z_m \}_{m=1}^M$, where $\mathcal{W}$ and $\mathcal{K}$ represent Walmart and Kmart, respectively. Here, $d_{jm}$ is firm $j$'s entry decision in market $m$, $x_{jm}$ includes firm-specific covariates, including a constant, and $z_m$ contains market-specific covariates. Table (ref) provides summary statistics for the sample used in our analysis.
For the purpose of illustrating our method, we model the entry game between Walmart and Kmart as a static game with incomplete information. Two players, Walmart ($\mathcal{W}$) and Kmart ($\mathcal{K}$), decide whether to enter a market. We assume that they make independent decisions across markets. Let $d_{j}=1$ if firm $j$ is active and $0$ otherwise. The payoff function of firm $j$ depends on its own productivity, whether its competitor enters or not, market- and firm-specific covariates, and private information: \[ u_{j}(d_{j},d_{-j})=\underbrace{X_{j}^\prime \beta - Z^\prime \gamma}_{\xi_{j}} - \Delta d_{-j} + \epsilon_{j1} , \] if $d_{j}=1$ and $=\epsilon_{j0}$ otherwise, where $X =(X_{\mathcal{W}},X_{\mathcal{K}})^{\prime}$ is firm characteristics that affect only the focal firm's profit and $Z$ is market characteristics common to both firms. For convenience, we denote $\xi_{j} = X_{j}^\prime \beta - Z^\prime \gamma$.
Firm $j$'s profit is $\xi_j$ under monopoly and $\xi_j - \Delta$ under duopoly. Note we allow asymmetry in monopoly profit by including a constant in firm-specific covariates. The term $Z^\prime \gamma$ is common among all firms. Denote $\theta = (\beta, \gamma, \Delta)^\prime$, market- and firm-specific characteristics $(x,z)$ are common knowledge, and firm $j$'s private information $\epsilon_{j}$ is type-1 extreme value distributed and independent of $\epsilon_{-j}$.
Therefore, the probability that firm $j$ chooses to enter is \[ p_{j} = \frac{1}{1+\exp \{- \xi_{j} + p_{-j} \Delta\}} , \] where $p_{-j}$ is its competitor's entry probability. Denote the CCPs as $p=(p_{\mathcal{W}},p_{\mathcal{K}})^{\prime}$. Define the best response mapping from CCP to CCP $\Psi: p\rightarrow p$. In equilibrium, we must have \[ p=\Psi(p,\theta). \]
We define the likelihood function as \[ \ell(p^\beta, \theta)= \sum_{j=\mathcal{W},\mathcal{K}}\sum_{m=1}^{M}\left\{ d_{jm} \log \Big[p_{j}^\beta(\chi_{m}) \Big] +(1-d_{jm} )\log \Big[1-p_{j}^\beta(\chi_{m}) \Big] \right\}, \] where $\chi_m = (x_{\mathcal{W} m}, x_{\mathcal{K} m}, z_m)^\prime$, the approximation structure $p_j^\beta$ will be introduced below, and the penalization as \[ \rho(\beta, \theta) = \sum_{j=\mathcal{W},\mathcal{K}}\sum_{m=1}^{M} \left[p_j^\beta(\chi_m) - \Psi_j\left(p_{-j}^\beta(\chi_m) , \theta\right)\right]^2 , \] which accounts for the equilibrium conditions for the set of observed market-specific covariates. In addition, we supply analytic gradient; see Appendix (ref).
Approximation Structure: We now consider the approximation of the CCPs. A naive approach is to approximate them as a flexible function of all market- and firm-specific covariates $p_j(x_j, x_{-j},z)$, which is of six dimensions in our empirical setting. To ensure that the approximation error disappears in first order asymptotics, the dimension of the approximation parameter $K$ needs to be large, leading to substantial computational challenges.
We propose a novel approximation structure that leverages the model structure: the deterministic component in the payoff function is linear in the parameters. As a result, how covariates $(x,z)$ determine the endogenous variable $p$ becomes a two-index model $p^*(\xi_j, \xi_{-j})$, which is much easier to approximate than a six-dimensional function. Using cubic basis functions following luo2018structural, we propose to approximate the CCPs in our empirical model by
where $s_\imath(\cdot)$ and $s_\jmath(\cdot)$ are cubic spline basis functions on $[0,1]$, $\beta = (\beta_{11},\ldots, \beta_{K K})^\prime$ and $\sigma(\cdot) = (1+e^{-\cdot})^{-1}$ representing the logistic function. In principle, we can use a different number of basis functions in the two dimensions. For convenience, we will use the same number $K$ and refer to it as the approximation parameter.
Note that the logistic function appears three times but for different reasons. First, because $s_\imath(\cdot)$ and $s_\jmath(\cdot)$ are cubic spline basis functions on $[0,1]$, the inner ones $\sigma(\xi_j)$ and $\sigma(\xi_{-j})$ transforms unbounded payoff indices $\xi_j$ and $\xi_{-j}$ into bounded ones on $[0,1]$. Interestingly, there are just the stand-alone entry probabilities when firms ignore competition. Second, the outer one transforms an approximation of the ex-ante value of entry $\sum_{\imath=1}^K\sum_{\jmath=1}^K\beta_{\imath \jmath} s_\imath \big( \cdot \big) s_{\jmath} \big( \cdot \big)$, before observing T1EV errors, into CCPs. Altogether, our approximation structure is a hybrid of a simple neural network and a tensor product linear sieve space. It leverages the index structure in the payoff function and hence reduces the dimension of the approximation parameters needed.
The algorithms proposed in Section (ref) share the same asymptotic properties and perform similarly in simple settings. For practical, real-world applications, we recommend breaking the search process into more manageable steps. Specifically, we suggest using the joint algorithm with a small smoothing parameter to identify good starting values, followed by the alternating iterative version of the nested algorithm for the main estimation. Table (ref) shows the estimated parameters when the number of basis functions, $K$, varies from 10 to 30. The second last column reports the maximum likelihood estimates assuming equilibrium uniqueness.\footnote{We conduct fine grid search to check equilibrium uniqueness and find that the equilibrium is unique in each market.}
The last column reports the two-step estimates. Two-step methods are generally less efficient than MLE and also rely on consistent first-stage estimates of the CCPs. Ideally, this first-stage estimation should be nonparametric, as the functional form of the solution is unknown, even when the profit and best response functions are known. However, this leads to the well-known curse of dimensionality. To address this, we propose a novel two-step approach that leverages the single-index structure, thereby avoiding the curse of dimensionality.\footnote{Similar strategies can be found in the econometrics literature on single-index regression models, such as stoker1986consistent and powell1989semiparametric.} Specifically, we first obtain the sieve MLE of the CCPs, $\widehat{p}_{j}$, using the same approximation structure as in Equation (ref), and then estimate the parameters by maximizing the pseudo-likelihood function \[ \max_{\theta} \quad \sum_{j=\mathcal{W},\mathcal{K}}\sum_{m=1}^{M}\left[d_{jm} \log\left\{p_{j}^\theta(\chi_{m})\right\} +(1-d_{jm} )\log\left\{1-p_{j}^\theta(\chi_{m})\right\}\right], \] where $p_{j}^\theta = \frac{1}{1+\exp \{- \xi_{j}(\theta) + \widehat{p}_{-j} \Delta\}}$.
The estimates are quite similar across different estimators. All estimates are significant at the 5% level and their signs are consistent with jia2008happens. More populated areas tend to have more stores, and higher retail sales per capita predict increased entry. Urbanized areas also attract more entry. The southern region dummy variable and the log of the distance to Walmart's headquarters in Bentonville, Arkansas, both significantly predict Walmart's entry decisions. Similarly, because Kmart’s headquarters are located in Troy, Michigan, the dummy variable for the Midwest region is predictive of Kmart's entry decisions.
The maximum likelihood estimates imply that the mean and standard deviation of $\xi_{\mathcal{W}}$ are -0.11 and 3.43, respectively, while the mean and standard deviation of $\xi_{\mathcal{K}}$ are -2.20 and 3.30, respectively. The large coefficients on firm dummies suggest substantial entry costs. Note that Walmart is a dominant firm with a penetration rate of 48%, while Kmart is relatively weak with a penetration rate of 19%. This explains the much lower coefficient on the Kmart firm dummy. The proposed sieve estimator performs well across different $K$. The larger the approximation parameter $K$ is, the closer the estimates become to the maximum likelihood estimates. The estimates from the proposed two-step estimator have larger biases.
A structural model is based on economic theory and describes how endogenous variables relate to a set of explanatory variables. This relationship is often expressed as an implicit function dependent on unknown parameters, which can be costly to solve. Two-step methods avoid solving the model but rely heavily on the accuracy of the first-step nonparametric estimation. We introduce SEES\ as a new class of estimators that use a sieve to approximate the solution while penalizing deviations from the equilibrium condition. SEES\ are straightforward to apply, at least as fast as alternative approaches like MLE, and more robust across various models. We believe our method will become a valuable tool in structural estimation.