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.
117,522 characters · 31 sections · 101 citation commands
Sequential algorithm for structural estimations with equilibrium constraints
Structural estimations with equilibrium constraints are prevalent in economics. Examples include dynamic discrete choice models (henceforth DDC; e.g., rust1987optimal), demand estimations (e.g., berry1995automobile, henceforth BLP), widely applied in industrial organization, labor economics, and other fields. Structural estimations with equilibrium constraints are also applied to general equilibrium models (e.g., ahlfeldt2015economics).
The traditional approach for estimating structural models with fixed point constraints is the Nested-Fixed Point (NFXP) algorithm. Generally, NFXP aims at exactly solving the following constrained continuous optimization problem:
In the problem above, $\theta$ denotes the structural parameters of the model, and $Y$ denotes the economic variables of the model.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {$Y$ is also known as nuisance parameters.} $G(Y;\theta)=0$ denotes the equilibrium constraint. The idea of the NFXP algorithm is to solve $G(Y;\theta)=0$ for $Y$ given candidate structural parameters $\theta$ in the inner loop, and search for the optimal $\theta$ that minimizes the objective function $Q\left(\theta;Y(\theta)\right)$ in the outer loop. However, as noted by su2012constrained, the computational burden of NFXP can be a large obstacle.
The current study investigates sequential algorithms possessing a property called “Zero-Jacobian Property” (we call ZJP hereinafter) in the previous studies, such as aguirregabiria2002swapping, kasahara2008pseudo (NPL for single-agent DDC estimations), and dearing2024efficient (EPL for dynamic game estimations).\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {NPL and EPL denote Nested Pseudo-Likelihood and Efficient Pseudo-Likelihood, respectively. Although not explicitly discussed, ABLP algorithm proposed by lee2015computationally for estimating static BLP (random-coefficient logit demand) model is also the sequential algorithm satisfying the ZJP.} In sequential algorithms, $\theta$ and $Y$ are iteratively updated. These studies discussed convergence properties and showed good performance of their algorithms. The current study generalizes these algorithms by developing a more unified framework to theoretically analyze the convergence properties of them. It then proposes a novel algorithm termed Sequential Linearly Constrained (SLC) algorithm, which can be easily applied to more general structural models than previous studies and often much faster than the NFXP algorithm. By conducting numerical experiments on dynamic discrete games with time\nobreakdash-varying unobserved heterogeneity and on dynamic BLP models, the present study demonstrates that the SLC can reduce computational time by several times relative to the NFXP.
Concerning the convergence properties of the sequential algorithms with the ZJP, this study demonstrates that their local convergence rate is nearly quadratic, provided the sample size is sufficiently large---even when a consistent initial estimate of $Y$ is unavailable. As it is unclear whether general methods to obtain initial consistent estimate of $Y$ for non-standard structural models are available, the non-dependence on the availability of initial consistent estimates is a notable strength. Furthermore, the current study also shows that the discussions hold not only for the Maximum Likelihood Estimation (MLE) but also for the Generalized Method of Moments (GMM), whereas the previous studies focused on the MLE. Both MLE and GMM share the property $\nabla_{Y}Q\left(\widehat{\theta},\widehat{Y}\right)=O_{p}\left(N^{-\frac{1}{2}}\right)$, where $\left(\widehat{\theta},\widehat{Y}\right)$ denotes the solution of ((ref)) and $N$ denotes the number of samples, which allows for a unified discussion of local convergence rates. Given that not a few structural estimation applications employ GMM, these results significantly broaden the applicability of sequential algorithms. The current study also shows that the estimator attains the same (statistical) efficiency as the NFXP estimator even after one iteration when the iteration is started from consistent initial values, in line with the previous methods utilizing ZJP.
Building on the unified framework, the present study proposes a novel SLC algorithm as an alternative to previously developed sequential methods. Each iteration is essentially equivalent to solving a surrogate constrained optimization problem in which the constraint is replaced by a linearization of the original constraint $G\left(Y;\theta\right)=0$.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {We can draw an analogy between the SLC update for constrained optimization and the classical one\nobreakdash-step Newton-Raphson update for unconstrained optimization. Both updates can yield asymptotically efficient estimates when supplied with consistent initial estimates for MLE. While the Newton--Raphson update solves a surrogate problem based on a quadratic approximation of the objective function, the SLC update instead solves a surrogate problem based on a linear approximation of the constraint, while continuing to use the original objective function.} The current study further discusses that the SLC algorithm can be implemented without explicitly computing the Jacobian of the equilibrium constraint function $G\left(Y;\theta\right)$, $\nabla_{Y}G\left(Y;\theta\right)$, which is very beneficial because computing Jacobians is very cumbersome and sometimes prone to coding mistakes.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {In recent years, automatic differentiation (AD) packages have become increasingly prevalent across many programming languages, allowing researchers to avoid manually coding the Jacobians of functions. However, this study finds that AD packages can sometimes be significantly slower than manually coded Jacobians, as illustrated by the examples in the Supplemental Appendix (ref). Although such performance drawbacks may be mitigated by future technological advancements, practitioners should remain aware of the current limitations of AD-based implementations.} In addition, the matrix storing $\nabla_{Y}G\left(Y;\theta\right)$ itself may require substantial memory, especially when the dimension of $Y$ is large. Thus, a Jacobian-free approach is highly desirable.
A major advantage of the SLC is that it does not rely on specific structures of the structural models, such as the existence of a fixed-point mapping with the ZJP (as in NPL) or the linearity of the equilibrium constraint $G(Y;\theta)$ concerning the structural parameters $\theta$ (as in EPL). Consequently, SLC is in principle applicable to a wider range of structural models. Note that we can also achieve global convergence of the SLC algorithm by introducing a line search procedure under moderate conditions, as demonstrated in the Supplemental Appendix (ref), although numerical results show that the SLC algorithm typically converges even without the procedure. SLC can be combined with fixed\nobreakdash-point iteration acceleration methods (e.g., Spectral), whose effectiveness has been highlighted in recent economic studies. Note that, under the linearity of $G\left(Y;\theta\right)$ concerning $\theta$, the SLC coincides with the EPL, whose effectiveness has been demonstrated in dynamic games without unobserved heterogeneity by dearing2024efficient.
Other than the NFXP approach and sequential algorithms, the “MPEC” (Mathematical Programming with Equilibrium Constraints) approach---typically relying on Lagrangian-based constrained optimization solvers---is also the alternative for structural estimation. The current study shows that the sequential algorithms with ZJP, especially the SLC algorithm, have similarities with the Sequential Quadratic Programming (SQP) algorithm, which is the basis of many solvers applied in the previous studies.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Most of the previous studies proposing the MPEC approach (e.g., su2012constrained; dube2012improving; egesdal2015estimating) applied commercial KNITRO solver, which builds on the Interior Point Method, in their numerical experiments.} Both the SQP and the SLC sequentially solve constrained optimization problems where the constraint equation corresponds to the linearization of the equilibrium constraint $G(Y;\theta)=0$. The key difference lies in the objective function of the constrained objective function: the SLC uses the original objective function $Q\left(\theta;Y\right)$, but the SQP uses a quadratic objective function using the candidate value of the Lagrange multiplier.
Compared to NFXP and “MPEC” (Lagrangian-based) approaches, the SLC algorithm offers several advantages. Table (ref) compares the NFXP, SLC, non-Lagrange-based sequential algorithms (e.g., EPL), and Lagrange-based sequential algorithms (e.g., SQP). Notably, SLC is more robust in the presence of multiple solutions to $G\left(Y;\theta\right)=0$, and it is often significantly faster than NFXP, as demonstrated in numerical experiments. Additionally, unlike Lagrangian-based methods, SLC can be implemented without explicitly computing the Jacobian $\nabla_{Y}G\left(Y;\theta\right)$.
The contribution of the current paper is twofold. First, it develops a unified theoretical framework for analyzing the convergence properties of sequential algorithms that possess the ZJP, extending prior work in this area. Second, it proposes the SLC algorithm, which is computationally more efficient, applicable to a broader class of structural models than previously proposed methods, and implementable with only minor modifications to the steps of the widely used NFXP algorithm.
The rest of this paper is organized as follows: Section (ref) explains the related literature. Section (ref) describes sequential algorithms using ZJP. Note that the overall discussions in Section (ref) does not rely on any statistical properties in the applications of structural estimations. Instead, it identifies key factors that contribute to the good behavior of the sequential algorithms, which are asymptotically satisfied in the applications of structural estimations. Section (ref) then examines the properties of the sequential algorithms in the context of structural estimation, incorporating formal statistical asymptotic analysis. Section (ref) presents numerical results, and Section (ref) finally concludes. Appendix (ref) provides additional discussions. Appendix (ref) shows the proofs of the statements in the main part of the current paper. The Supplemental Appendix of the current paper provides further discussions, including the introduction of inequality constraints on $\theta$ such as nonnegative constraints, and the global convergence results concerning the SLC.
Sequential algorithms with ZJP, including the SLC, are closely related to the NFXP approach, as both approaches aim at solving the constrained optimization problem. It is important to emphasize that the current study does not intend to claim that the NFXP approach is always inferior to sequential algorithms. Prior research has developed techniques to accelerate the convergence of the inner-loop iterations in specific models,\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {In the context of the BLP estimation, conlon2020best and fukasawa2024fast highlighted that introducing acceleration methods (e.g., SQUAREM, Spectral, Anderson acceleration) significantly improves the convergence speed of BLP contraction mapping iterations. fukasawa2024fast further proposed fast and simple fixed-point mappings for the inner loop of static and dynamic BLP estimation procedures. For dynamic models, bray2019strong introduced a relative value function iteration algorithm. iskhakov2016comment employed Newton's method to efficiently implement the inner loop of NFXP algorithms in single-agent DDC models.} and should consider these methods where applicable. However, such techniques are not always available or well-investigated, especially in the case of non-standard structural models. The SLC algorithm proposed in the current study is in principle applicable to any models, and is worth considering as the alternative to the NFXP approach. In Appendix (ref), we discuss in what cases the SLC would yield better performance than the NFXP in terms of the computational cost.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Previous studies have highlighted the computationally inferior performance of the “MPEC” approach relative to NFXP in certain models. Notable examples include lee2015computationally, pal2023comparing for static BLP models, and sun2019computationally for dynamic BLP models. iskhakov2016comment demonstrated that the NFXP algorithm can be highly competitive with the MPEC approach in single-agent DDC models, provided that the inner and outer loops are appropriately designed.}
Other than the NFXP, MPEC, and sequential algorithms, two-step estimation methods have been proposed for certain standard models. Examples include hotz1993conditional and others' Conditional Choice Probability (CCP)-based methods for DDC models, and lu2023semi for static BLP models.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Bayesian-type methods (e.g., imai2009bayesian; gallant2022constrained) are also proposed to mitigate computational problems in the literature.} Although convenient to use and usually computationally less demanding, they typically require large sample sizes to consistently estimate first-stage parameters. In small samples, bias can be substantial. In addition, in structural models where such methods have yet to be developed (e.g., dynamic BLP, to my knowledge), they are not applicable. Note that such methods can be used as the initial values in the sequential algorithms with ZJP to obtain an statistically efficient estimator.
The current study also relates to studies on sequential estimation methods that do not necessarily satisfy the ZJP (e.g., NPL for dynamic games; cf. aguirregabiria2007sequential). \ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {bugni2021iterated also proposed a sequential estimation method using a minimum distance criterion function, and showed that the estimator is efficient if an optimal sequence of weight matrices is taken.} Compared to these, sequential algorithms with ZJP offer advantages in terms of statistical efficiency and theoretical guarantees on local convergence speed.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Computing standard errors for sequential estimators that do not exhibit the ZJP can be analytically complex in practice (cf. yamaguchi2019effects), often requiring reliance on computationally demanding bootstrap-based methods. In contrast, sequential estimators with the ZJP attain the same asymptotic efficiency as the NFXP estimator, and their standard errors can be computed in the same manner as those for the NFXP, without resorting to bootstrap-based procedures.}\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {As discussed in kasahara2012sequential, NPL algorithm for dynamic games may not be locally convergent even in large samples. To tackle the problem, kasahara2012sequential introduced some techniques such as the relaxation method for fixed-point iterations, and aguirregabiria2021imposing proposed the use of the spectral algorithm.} While this paper focuses on ZJP-based algorithms, the insights may also be useful for algorithms without ZJP, which are sometimes simpler to implement.
The idea of the sequential algorithms with ZJP is closely related to the Neyman orthogonality, which has been extensively studied in the recent econometrics and machine learning literature. By constructing Neyman orthogonal scores or moment conditions for MLE/GMM, we can obtain debiased estimates of $\theta$ when the appropriate consistent estimate of $Y$ is available.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {chernozhukov2018double presents a general procedure to construct Neyman orthogonal scores/moments. chernozhukov2022locally and sawadogo2025efficient propose methods for estimating DDC models without unobserved heterogeneity using Neyman orthogonality.} Concerning the sequential algorithms with ZJP, we can easily show that the gradient of the surrogate objective function used in the algorithm satisfies the Neyman orthogonality condition in the MLE case. Even in the GMM setting, the Neyman orthogonality condition asymptotically holds.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {dearing2024efficient showed that the EPL asymptotically satisfies the Neyman orthogonality condition for MLE. The current study further advances the discussions by demonstrating that the Neyman orthogonality holds even under non-asymptotic settings for MLE.} This suggests that the first order condition of the optimization problem to be solved in the sequential algorithm is less likely to be affected by the value of the nuisance parameters around true parameter values. Hence, we can obtain updated parameters which would be close to the true value, even if the candidate nuisance parameter values $\gamma_{k}$ are not very close to the true values.
Although exploring Neyman orthogonality-based methods is a promising direction, they require consistent initial estimates of the nuisance parameter $Y$ which may be difficult to obtain in complex structural models. Additionally, constructing Neyman orthogonal functions often involves computing second derivatives of the objective function $Q$, which may not be straightforward. In contrast, the sequential algorithms proposed here---especially the SLC algorithm---can be implemented without consistent initial estimates or explicit derivative computations. Thus, they offer a practical alternative to Neyman orthogonality-based approaches for structural estimations with equilibrium constraints.
The current section discusses sequential algorithms using a mapping with the zero-Jacobian property (ZJP). Note that the overall discussions in the current section does not rely on any statistical discussions specific to the application of structural estimations. Instead, we briefly discuss factors which contribute to the good behavior of the sequential algorithms, which are asymptotically satisfied in the application of structural estimations, as described in detail in Section (ref).
Structural estimation with equilibrium constraints basically solves the following continuous constrained optimization problem:
where $\theta$ denotes the structural parameters of the model, and $Y$ denotes economic variables of the model. $n_{\theta}$ and $n_{Y}$ denotes the dimension of the variables. Here, we impose the following assumptions:
Let $\widehat{\gamma}\equiv\left(\widehat{\theta},\widehat{Y}\right)$ be the solution of the problem. Then, the following statement holds:
Here, we define \[ \Gamma_{FOC}\equiv\left\{ \gamma\equiv\left(\theta,Y\right):G\left(Y;\theta\right)=0\ \text{and}\ \nabla_{\theta}Q\left(\theta,Y\right)-\left(\nabla_{Y}Q\left(\theta,Y\right)\right)\left(\nabla_{Y}G\left(Y;\theta\right)\right)^{-1}\left(\nabla_{\theta}G\left(Y;\theta\right)\right)=0\right\} . \] Then, $\widehat{\gamma}\in\Gamma_{FOC}$ holds under Assumptions (ref) and (ref).
The original constrained optimization problem $Q\left(\theta,Y\right)\ s.t.\ G\left(Y,\theta\right)=0$ can be reformulated as the following problem, as discussed in su2012constrained:
Let $\mathcal{Y}\left(\theta\right)=\arg\min_{Y(\theta)\ s.t.\ G(Y(\theta);\theta)=0}Q\left(\theta,Y(\theta)\right)$. If we define $Q_{NFXP}\left(\theta\right)\equiv Q\left(\theta,\mathcal{Y}(\theta)\right)$,\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {We cannot exclude the possibility that $\mathcal{Y}\left(\theta\right)$ is empty, especially when $\theta$ is far from $\widehat{\theta}$. In this case, we cannot define $Q_{NFXP}$$\left(\theta\right)$ at the parameter value.} the problem can be written as the unconstrained optimization problem $\min_{\theta\in\mathbb{R}^{n_{\theta}}}Q_{NFXP}\left(\theta\right)$. Note that, in the neighborhood of $\widehat{\theta}$, the function $\mathcal{Y}\left(\theta\right)$ is $C^{2}$ class by the implicit function theorem and Assumption (ref), and it implies that $Q_{NFXP}\left(\theta\right)$ is $C^{2}$ class. Here, let $\theta_{NFXP}$ be the solution of $\min_{\theta}Q_{NFXP}\left(\theta\right)$. In addition, the following statement shows the relationship between the solution of the original constrained optimization problem and the NFXP problem, as shown by su2012constrained:
Proposition (ref) and Assumption (ref) imply $\theta_{NFXP}=\widehat{\theta}$, and $\theta_{NFXP}$ is the unique minimizer of $Q_{NFXP}\left(\theta\right)$.
In this study, we consider sequential algorithms which can solve the constrained optimization problem ((ref)). Algorithm (ref) shows the general steps of the algorithms.
Here, $\Upsilon\left(\theta;\gamma_{k}\right)$ denotes a mapping which is a function of $\theta$ and $\gamma_{k}=(\theta_{k},Y_{k})$. For general forms of $\Upsilon\left(\theta;\gamma_{k}\right)$, the set of the fixed points of Algorithm (ref) may not include the solution of ((ref)). Hence, we impose the following assumption:
Let $\Gamma_{seq}$ be the set of fixed points of Algorithm (ref). Then, the following statement holds:
So far, the following forms of $\Upsilon$ have been proposed in the literature:\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {dearing2024efficient also mention the specification $\Upsilon\left(\theta;\gamma_{k}\right)=Y_{k}-\left(Z\left(Y_{k};\theta_{k}\right)\right)^{-1}\left(G\left(Y_{k};\theta\right)\right)$, where $Z$ is a continuously differentiable function and $Z\left(\theta,Y\right)=\nabla_{Y}G\left(\theta,Y\right)$ for all $(\theta,Y)$ such that $G\left(Y;\theta\right)=0$, in Appendix C of their paper.}
dearing2024efficient showed that the EPL mapping $\Upsilon\left(\theta;\gamma_{k}\right)\equiv Y_{k}-\left(\nabla_{Y}G\left(Y_{k};\theta_{k}\right)\right)^{-1}\left(G\left(Y_{k};\theta\right)\right)$ satisfies the Assumption (ref)'s three conditions in Lemma 2 of their paper. They also mentioned in Appendix C of their paper that the mapping $\Upsilon\left(\theta;\gamma_{k}\right)\equiv Y_{k}-\left(\nabla_{Y}G\left(Y_{k};\theta\right)\right)^{-1}\left(G\left(Y_{k};\theta\right)\right)$ also satisfies the three conditions. We can also easily show that the NPL mapping satisfies the conditions, as shown in Appendix (ref).
The current study proposes the following mapping:
As shown in Appendix (ref), the mapping satisfies the three conditions in Assumption (ref). The mapping can be considered as a slight extension of the EPL mapping. In the EPL mapping $\Upsilon\left(\theta;\gamma_{k}\right)\equiv Y_{k}-\left(\nabla_{Y}G\left(Y_{k};\theta\right)\right)^{-1}\left(G\left(Y_{k};\theta\right)\right)$, the Taylor approximation of $G\left(Y_{k};\theta\right)$ around $\theta_{k}$ is $G\left(Y_{k};\theta_{k}\right)+\left(\nabla_{\theta}G\left(Y_{k};\theta_{k}\right)\right)\left(\theta-\theta_{k}\right)$, and we can obtain the representation ((ref)) by using $G\left(Y_{k};\theta_{k}\right)+\left(\nabla_{\theta}G\left(Y_{k};\theta_{k}\right)\right)\left(\theta-\theta_{k}\right)$ instead of $G\left(Y_{k};\theta\right)$ in the EPL mapping. Note that the EPL mapping and the mapping ((ref)) coincide when $G\left(Y;\theta\right)$ is linear concerning $\theta$.
The mapping ((ref)) has a computational advantage over the EPL mapping, in that $\left(\nabla_{Y}G\left(Y_{k};\theta_{k}\right)\right)^{-1}\left(G\left(Y_{k};\theta_{k}\right)\right)$ and $\left(\nabla_{Y}G\left(Y_{k};\theta_{k}\right)\right)^{-1}\left(\nabla_{\theta}G\left(Y_{k};\theta_{k}\right)\right)$ can be precomputed before solving the optimization problem \[ \min_{\theta\in\mathbb{R}^{n_{\theta}}}Q\left(\theta,Y_{k}-\left(\left(\nabla_{Y}G\left(Y_{k};\theta_{k}\right)\right)^{-1}\left(G\left(Y_{k};\theta_{k}\right)\right)\right)-\left(\left(\nabla_{Y}G\left(Y_{k};\theta_{k}\right)\right)^{-1}\left(\nabla_{\theta}G\left(Y_{k};\theta_{k}\right)\right)\right)\left(\theta-\theta_{k}\right)\right). \] In contrast, the EPL mapping in general requires repeated evaluations of $\left(\nabla_{Y}G\left(Y_{k};\theta\right)\right)^{-1}\left(G\left(Y_{k};\theta\right)\right)$ and $\left(\nabla_{Y}G\left(Y_{k};\theta_{k}\right)\right)^{-1}\left(\nabla_{\theta}G\left(Y_{k};\theta_{k}\right)\right)$, which necessitate solving a linear equation given candidate $\theta$.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Computational cost in the optimization step using the ABLP mapping is higher because of the additional need to compute $\nabla_{Y}G\left(Y_{k};\theta\right)$ for each candidate value of $\theta$.}
We call the proposed mapping ((ref)) SLC (Sequential Linearly Constrained) mapping, because we can easily show that the sequential algorithm (ref) using the SLC mapping is equivalent to sequentially solving the following optimization problem with linear constraints:
The constraint is the linearization of the original constraint $G\left(Y;\theta\right)=0$ around $(\theta_{k},Y_{k})$. The representation is convenient for comparing the SLC algorithm with the Lagrangian-based algorithms, such as the SQP algorithm. More specifically, the difference with the SLC and the SQP is the objective function in each iteration: SLC uses the original objective $Q\left(\theta,Y\right)$, but the SQP uses the quadratic approximation of the Lagrangian function. For details, see Appendix (ref).
To make the algorithm well-defined, we impose the following assumption:
Note that $\nabla_{Y}G\left(Y;\theta\right)$ is nonsingular at least near the solution $\left(\widehat{Y},\widehat{\theta}\right)$ under Assumption (ref). In models where $G$ is derived from a contractive fixed-point constraint (e.g., single-agent DDC, Static BLP (random coefficient logit) model), we can theoretically show the non-singularity at any points.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {If $G\left(Y;\theta\right)=Y-\Phi\left(Y;\theta\right)$, $\nabla_{Y}G\left(Y;\theta\right)=I-\nabla_{Y}\Phi\left(Y;\theta\right)$ holds. If $\Phi$ is a contraction on $Y$, the spectral radius of $\nabla_{Y}\Phi\left(Y;\theta\right)$ should be smaller than 1, and $\nabla_{Y}G\left(Y;\theta\right)=I-\nabla_{Y}\Phi\left(Y;\theta\right)$ is nonsingular (cf. Neumann series).} Furthermore, the assumption can be weaker than the assumption that at least one solution of $G\left(Y;\theta\right)=0$ concerning $Y$ exists given $\theta$, which the NFXP algorithm relies on.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Nonsingular Jacobian assumption is sometimes imposed in globally convergent nonlinear equation solution algorithms (e.g., li2000derivative).}
Other than the advantage of the smaller computational costs in the optimization step among the sequential algorithms with the ZJP, the proposed SLC algorithm has an advantage in that the updating direction is descending given an appropriate merit function.
Here, we define a $l_{1}$-norm merit function $\phi_{1}\left(\theta,Y;\mu\right)\equiv Q\left(\theta,Y\right)+\mu\left\Vert G\left(Y;\theta\right)\right\Vert _{1}$ given a tuning parameter $\mu>0$. As shown in Theorem 17.3 of nocedal2006numerical, $\widehat{\gamma}$ is a local minimizer of $\phi_{1}\left(\gamma;\mu\right)$ if $\mu\geq\left\Vert \lambda\right\Vert _{\infty}$, where $\lambda$ denotes the vector of Lagrange multipliers of the constrained optimization problem ((ref)). Concerning the merit function, the following statement holds:
The statement implies that $\phi_{1}\left(\theta_{k}+\alpha d_{\theta,k},Y_{k}+\alpha d_{Y,k};\mu\right)<\phi_{1}\left(\theta_{k},Y_{k};\mu\right)$ if $\alpha>0$ is sufficiently small. The property can be utilized to show the global convergence of the SLC algorithm by introducing a line search step, as discussed in the Supplemental Appendix (ref).\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {The strategy of the proof is similar to the one for the SQP algorithm (cf. bonnans2006numerical), which utilizes the property $\left(\nabla_{\gamma}G\left(\gamma_{k}\right)\right)\left(\gamma-\gamma_{k}\right)+G\left(\gamma_{k}\right)=0$ in each iteration.}
Here, we consider the convexity of $\widetilde{Q}$ in relations to the NFXP. The following lemma shows the statement.
Because $\widehat{\theta}$ is the unique minimizer of $Q_{NFXP}\left(\theta\right)$, $\nabla_{\theta\theta^{\prime}}Q_{NFXP}\left(\widehat{\theta}\right)$ is positive definite. As discussed in Section (ref) in the applications to structural estimations, $\frac{\partial Q\left(\widehat{\theta},\widehat{Y}\right)}{\partial Y}\approx0$ holds in large samples. Consequently, $\nabla_{\theta\theta^{\prime}}\widetilde{Q}\left(\widehat{\theta};\widehat{\gamma}\right)$ is close to positive definite if the number of samples is large.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Generally, $\nabla_{\theta\theta^{\prime}}\widetilde{Q}\left(\widehat{\theta};\widehat{\gamma}\right)$ may not be positive definite. Supplemental Appendix (ref) investigates the conditions where the function $\widetilde{Q}\left(\theta;\gamma\right)$ is convex. }
Let $H\left(\gamma\right)\equiv\left(
\right)\equiv\left(
\right)$, where $\widetilde{\theta}\left(\gamma\right)\equiv\arg\min_{\theta\in\Theta}\widetilde{Q}\left(\theta;\gamma\right)=Q\left(\theta,\Upsilon\left(\theta;\gamma\right)\right)$. Then,
holds. Here, let $\widetilde{\gamma}\in\Gamma_{seq}$. Note that we allow for $\widetilde{\gamma}\neq\widehat{\gamma}\in\Gamma_{seq}$. By $\nabla_{\theta}\widetilde{Q}\left(\widetilde{\theta}(\gamma),\gamma\right)=0$, $\nabla_{\gamma}\widetilde{\theta}\left(\gamma\right)=\left(\nabla_{\theta\theta^{\prime}}\widetilde{Q}\left(\widetilde{\theta}\left(\gamma\right),\gamma\right)\right)^{-1}\left(\nabla_{\theta\gamma^{\prime}}\widetilde{Q}\left(\widetilde{\theta}\left(\gamma\right),\gamma\right)\right)$ holds by the implicit function theorem if the Hessian matrix $\nabla_{\theta\theta^{\prime}}\widetilde{Q}\left(\widetilde{\theta}\left(\gamma\right),\gamma\right)$ is invertible.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {$\nabla_{\theta\theta^{\prime}}\widetilde{Q}\left(\widetilde{\theta}\left(\gamma\right),\gamma\right)$ is invertible if the matrix is positive definite.} By $\widetilde{Q}\left(\theta,\gamma\right)=Q\left(\theta,\Upsilon\left(\theta,\gamma\right)\right)$, $\nabla_{\gamma}\widetilde{Q}\left(\theta,\gamma\right)=\left(\nabla_{Y}Q\left(\theta,\Upsilon\left(\theta,\gamma\right)\right)\right)\left(\nabla_{\gamma}\Upsilon\left(\theta,\gamma\right)\right)$ holds, and hence:
By Taylor's theorem, we have $H\left(\gamma_{k-1}\right)-H\left(\widetilde{\gamma}\right)=\left(\nabla_{\gamma}H\left(\widetilde{\gamma}\right)\right)\left(\gamma_{k-1}-\widetilde{\gamma}\right)+O\left(\left\Vert \gamma_{k-1}-\widetilde{\gamma}\right\Vert ^{2}\right)$. By $H\left(\gamma_{k-1}\right)=\gamma_{k}$ and $H\left(\widetilde{\gamma}\right)=\widetilde{\gamma}$,
Because $\nabla_{\gamma}\Upsilon\left(\widetilde{\theta};\widetilde{\gamma}\right)=0$ holds under Assumption (ref), ((ref)) and $\Upsilon\left(\widetilde{\theta};\widetilde{\gamma}\right)=\widetilde{Y}$ imply:
Hence, we obtain the following statement:
The proposition implies that $\gamma_{k}-\widetilde{\gamma}\approx O\left(\left\Vert \gamma_{k-1}-\widetilde{\gamma}\right\Vert ^{2}\right)$ holds, which implies $\gamma_{k}$ locally converges to $\widetilde{\gamma}$ near-quadratically under $\nabla_{Y}Q\left(\widehat{\theta},\widehat{Y}\right)\approx0$.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {The definition of quadratic convergence is: $\gamma_{k}$ converges to $\widetilde{\gamma}$ quadratically if $\lim_{k\rightarrow\infty}\frac{\left\Vert \gamma_{k+1}-\widetilde{\gamma}\right\Vert }{\left\Vert \gamma_{k}-\widetilde{\gamma}\right\Vert ^{2}}<\infty$.} As we will discuss in Section (ref), the condition approximately holds at the global minimum $\widehat{\theta}$ in large samples in the application of structural estimations.
SLC algorithm is often computationally more attractive compared to the NFXP approach. Here, we clarify the point by comparing the computational costs of the SLC and the NFXP using the Newton's method as the inner-loop solution method and a gradient-based optimization method (e.g., BFGS) as the outer-loop algorithm. Note that the Newton's method (Newton-Kantorovich iteration) is the fundamental method to solve nonlinear equations, which is known to be locally quadratically convergent.
In the case of the SLC, we need to compute the values of $\left(\nabla_{Y}G\left(Y_{k};\theta_{k}\right)\right)^{-1}\left(G\left(Y_{k};\theta_{k}\right)\right)$ and $\left(\nabla_{Y}G\left(Y_{k};\theta_{k}\right)\right)^{-1}\left(G\left(Y_{k};\theta_{k}\right),\nabla_{\theta}G\left(Y_{k};\theta_{k}\right)\right)$ only once in each iteration.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {The computational cost of computing $\left(\nabla_{Y}G\left(Y_{k};\theta_{k}\right)\right)^{-1}\left(G\left(Y_{k};\theta_{k}\right),\nabla_{\theta}G\left(Y_{k};\theta_{k}\right)\right)$ is not largely different from the one for computing $\left(\nabla_{Y}G\left(Y_{k};\theta_{k}\right)\right)^{-1}G\left(Y_{k};\theta_{k}\right)$ for moderate $n_{\theta}$ when relying on the exact solution method for solving linear equations. Since the $O(n^{3})$ cost of matrix factorization (e.g., LU decomposition) is the dominant bottleneck in solving linear equations, where $n$ is the size of the matrix, adding additional right-hand side terms does not significantly increase the overall computational burden.} In contrast, in the case of the NFXP with the Newton's method,\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Newton's method solves for the solution of $G(Y;\theta)=0$ by repeatedly applying $Y_{m+1}\leftarrow Y_{m}-\left(\nabla_{Y}G\left(Y_{m};\theta_{m}\right)\right)^{-1}\left(G\left(Y_{m};\theta_{m}\right)\right)$ until convergence. It attains quadratic local convergence in general model settings. Note that sometimes utilizing good fixed-point mappings and fixed-point mapping acceleration methods (e.g., Anderson acceleration methods) lead to superior performance than the Newton's method, because the computation of the Jacobian $\nabla_{Y}G\left(Y;\theta\right)$ in the the Newton's method is sometimes costly.} we need to compute the values of $\left(\nabla_{Y}G\left(Y;\theta\right)\right)^{-1}\left(G\left(Y;\theta\right)\right)$ many times to solve for $Y$ such that $G\left(Y;\theta\right)=0$ in each outer-loop iteration. In addition, we need to compute the values of $\left(\nabla_{Y}G\left(Y_{k};\theta_{k}\right)\right)^{-1}\left(\nabla_{\theta}G\left(Y_{k};\theta_{k}\right)\right)$ once to compute the gradient of the objective function $Q_{NFXP}\left(\theta\right)$. Hence, concerning the NFXP algorithm with the Newton's method for the inner-loop, the repeated evaluations of $\left(\nabla_{Y}G\left(Y;\theta\right)\right)^{-1}\left(G\left(Y;\theta\right)\right)$ can lead to larger computational cost than the SLC algorithm in each iteration.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {I thank Victor Aguirregabiria for suggesting the insight.}
Furthermore, as demonstrated in the numerical experiments, the current study finds that the SLC often require fewer number of iterations than the NFXP with first-order gradient-based algorithms (e.g., BFGS, Interior Point method). It can also contribute to the computationally superior performance of the SLC.
Note that the SLC involves the optimization step $\arg\min_{\theta\in\mathbb{R}^{n_{\theta}}}\widetilde{Q}\left(\theta;\gamma_{k}\right)$ in each iteration, which does not appear in the NFXP algorithm. However, as long as the optimization step is computationally much less costly, the SLC attains computationally superior performance than the NFXP. Concerning the inner-loop of the NFXP, Newton's method is sometimes inferior to good fixed-point iteration with acceleration methods (e.g., Anderson acceleration), because the computations of $\nabla_{Y}G\left(Y;\theta\right)$ and the solution of a linear equation $\left(\nabla_{Y}G\left(Y;\theta\right)\right)^{-1}\left(G\left(Y;\theta\right)\right)$ are sometimes costly. Still, numerical experiments in Section (ref) suggest that the SLC can outperform the NFXP with fixed-point iteration-based inner-loop. In Appendix (ref), we discuss more about the comparison of the computational costs. Note that the Newton's iteration is known to be unstable without regularization when the initial value is far from the solution. However, the SLC is usually stable, as demonstrated by numerical experiments. Supplemental Appendix (ref) provides further theoretical discussions with additional numerical experiments.
One of the advantages the SLC and EPL algorithms possess is that they can be implemented without explicitly computing the Jacobians of the equilibrium constraint $\nabla_{Y}G\left(Y;\theta\right)$.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {As shown in Table (ref), it is not clear whether the Lagrangian-based algorithms (e.g. SQP) possess the property.} Generally, coding the Jacobian $\nabla_{Y}G\left(Y;\theta\right)$ is prone to mistakes especially when the function $G\left(Y;\theta\right)$ is complex. In addition, $\nabla_{Y}G\left(Y;\theta\right)$ is a $n_{Y}\times n_{Y}$-dimensional matrix, which implies the memory requirement for storing the matrix $\nabla_{Y}G\left(Y;\theta\right)$ increases quadratically, and the algorithm may not be applicable for high-dimensional problems when the the memory of the computers practitioners use is not necessarily sufficient. In this subsection, we discuss the Jacobian-free implementation of the SLC and EPL algorithms, based on the discussion in fukasawa2025jacobian that proposed the Jacobian-free implementation of the EPL algorithm.
In the SLC and EPL algorithms, we need to solve a linear equation $\left(\nabla_{Y}G\left(Y;\theta\right)\right)x=b$ for $x\in\mathbb{\mathbb{R}}^{n_{Y}}$ to compute the value of $\left(\nabla_{Y}G\left(Y;\theta\right)\right)^{-1}b$ for a vector $b\in\mathbb{R}^{n_{Y}}$. Then, we can avoid the computation of $\nabla_{Y}G\left(Y;\theta\right)$ by utilizing the ideas of numerical derivatives and the Krylov-based algorithms for solving a linear equation. In the Krylov-based algorithms (e.g, GMRES), it is sufficient to compute the values of $\left(\nabla_{Y}G\left(Y;\theta\right)\right)v$ for any $v\in\mathbb{R}^{n_{\theta}}$ to solve a linear equation $\left(\nabla_{Y}G\left(Y;\theta\right)\right)x=b$. Concerning the Jacobian-vector product $\left(\nabla_{Y}G\left(Y;\theta\right)\right)v$, it can be approximated by $\left(\nabla_{Y}G\left(Y;\theta\right)\right)v\approx\frac{\left(\nabla_{Y}G\left(Y+\epsilon v;\theta\right)\right)-\left(\nabla_{Y}G\left(Y-\epsilon v;\theta\right)\right)}{2\epsilon}$ for a sufficiently small $\epsilon>0$ using the idea of numerical derivatives. Consequently, we can solve $\left(\nabla_{Y}G\left(Y;\theta\right)\right)x=b$ for $x\in\mathbb{R}^{n_{Y}}$ without explicitly computing the Jacobian $\nabla_{Y}G\left(Y;\theta\right)$. Note that the Krylov-based algorithms for solving linear equations, such as the GMRES algorithm, are available in many programming languages (e.g., MATLAB, Julia, Python), and the implementation is relatively straightforward as long as we can prepare the function $G\left(Y;\theta\right)$. For details, see also fukasawa2025jacobian, proposing the Jacobian-free approach for the EPL.
In this section, we consider sequential algorithms for solving constrained optimization problems in the context of statistical problems, especially the structural estimations in economics. In structural estimation applications, $\theta$ denotes structural parameters such as utility parameters, and $Y$ denotes economic variables, such as value functions or choice-specific value functions in DDC estimations. Let $\theta^{*}$ be the true parameter value.
To simplify discussions, we consider the setting where $Y$ is a deterministic variable, and it does not depend on the observed data. We also assume the equilibrium constraint $G\left(Y;\theta\right)=0$ does not directly depend on the observed data.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {The setting applies to standard dynamic discrete choice models (single-agent/games). Strictly speaking, the setting does not directly apply to BLP models, because some variables such as unobserved product characteristics generally depend on the observed data, and they are stochastic. However, considering the setting where $Y$ is a stochastic variable further complicates the discussions. Hence, we focus our attention on the setting with deterministic $Y$. Note that aguirregabiria2026nested recently demonstrated a proof of a consistency of a sequential estimation-based estimator in the context of the static BLP allowing for stochastic $Y$, although the ZJP may not hold.} Let $Y^{*}$ be the true value of $Y$ consistent with $\theta^{*}$, which satisfies $G\left(Y^{*};\theta^{*}\right)=0$. We assume the observed data is generated based on $\left(\theta^{*},Y^{*}\right)$, and the objective function $Q\left(\theta,Y\right)$ depends on the observed data $\left\{ w_{i};i=1,\cdots,N\right\} $.
In this section, we consider two types of statistical specifications widely applied in the literature of economics: Maximum Likelihood Estimation (MLE) and Generalized Method of Moments (GMM):\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Although not explicitly discussed, Minimum Distance (MD) Estimation can be treated analogous to the GMM estimation.}
We assume $w_{i}$ is generated from $f^{*}\left(w|\theta^{*},Y^{*}\right)=f^{*}\left(w\right)$. Let $Q^{*}\left(\theta,Y\right)\equiv E\left[Q\left(\theta,Y\right)\right]$.
where $\overline{m}\left(\theta,Y\right)=\frac{1}{N}\sum_{i=1}^{N}m\left(w_{i};\theta,Y\right)$ and $m\left(w_{i};\theta,Y\right)$ depends on the data of sample $i$. Let $Q^{*}\left(\theta,Y\right)\equiv\left(E\left[m\left(w;\theta,Y\right)\right]\right)^{\prime}W\left(E\left[m\left(w;\theta,Y\right)\right]\right)$.
Here, we impose the following assumption:
As discussed in Section (ref), the value of $\nabla_{Y}Q\left(\widehat{\theta},\widehat{Y}\right)$ largely affects the local convergence speed and the local convexity of the objective function. Hence, this subsection investigates the value of $\nabla_{Y}Q\left(\widehat{\theta},\widehat{Y}\right)$. The following discussions show that $\nabla_{Y}Q\left(\widehat{\theta},\widehat{Y}\right)\rightarrow_{p}0$, and more specifically $\nabla_{Y}Q\left(\widehat{\theta},\widehat{Y}\right)=O_{p}\left(N^{-\frac{1}{2}}\right)$, holds if the model is correctly specified and we apply the MLE or GMM.
In the large sample environment, $\nabla_{Y}Q^{*}\left(\theta^{*};Y^{*}\right)=0$ holds. In the MLE and GMM, we can choose $Q^{*}$ such that $Q^{*}\left(\theta,Y\right)\geq0\ \forall\left(\theta,Y\right)$ and $Q^{*}\left(\theta^{*},Y^{*}\right)=0$. Concerning the MLE, it is well known that the MLE is equivalent to minimizing the Kullback-Leibler Information Criterion (KLIC),\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {The log likelihood $\int\ln f\left(w|\theta,Y\right)f^{*}\left(w\right)dw$ satisfies:
Here, we defined $KLIC\left(f^{*};f\left(\cdot;\theta,Y\right)\right)\equiv\int\left[\ln\frac{f^{*}\left(w|\theta^{*},Y^{*}\right)}{f\left(w|\theta,Y\right)}f^{*}\left(w|\theta^{*},Y^{*}\right)\right]dw$. By the property of the KLIC, $KLIC\left(f^{*};f\left(\cdot;\theta,Y\right)\right)\geq0$ holds for all $\left(\theta,Y\right)$, and the equality holds when $\left(\theta,Y\right)=\left(\theta^{*},Y^{*}\right)$. Then, because the second term in the right hand side does not depend on $(\theta,Y)$, the minimizer of the KLIC maximizes the likelihood, as long as $G\left(Y^{*};\theta^{*}\right)=0$ holds.} and the conditions hold if we alternatively choose $Q^{*}$ as the KLIC. Regarding the GMM, the conditions hold because $Q^{*}\left(\theta,Y\right)=\left(E\left[m\left(w;\theta^{*},Y^{*}\right)\right]\right)^{T}W\left(E\left[m\left(w;\theta^{*},Y^{*}\right)\right]\right)=0$, where $E\left[m\left(w;\theta^{*},Y^{*}\right)\right]=0$ and $W$ is a positive definite matrix. Under the conditions, we can easily show that $\nabla_{Y}Q^{*}\left(\theta^{*};Y^{*}\right)=0$ holds because of the first-order condition, as long as $G\left(\theta^{*},Y^{*}\right)=0$ holds.
In finite sample settings, the equality may not hold, but $\nabla_{Y}Q\left(\widehat{\theta},\widehat{Y}\right)$ would be close to 0 when the number of samples $N$ is large.
By differentiating the both sides of $1=\int f\left(w|\gamma\right)dw$ with regard to $Y$ at $\gamma^{*}\equiv\left(\theta^{*},Y^{*}\right)$, we have $0=\int\left(\nabla_{Y}f\left(w|\gamma^{*}\right)\right)dw=\int\left(\nabla_{Y}\ln f\left(w|\gamma^{*}\right)\right)f\left(w|\gamma^{*}\right)dw$. Then,
{
}
Concerning the first term, $\nabla_{Y}\ln f\left(w_{i}|\widehat{\gamma}\right)-\nabla_{Y}\ln f\left(w_{i}|\gamma^{*}\right)=\left(\nabla_{Y\gamma^{\prime}}\ln f\left(w_{i}|\overline{\gamma}\right)\right)\left(\gamma^{*}-\gamma^{*}\right)$ holds by the mean value theorem, where there exists $c\in[0,1]$ such that $\overline{\gamma}=c\widehat{\gamma}+(1-c)\gamma^{*}$. Then, by $\gamma^{*}-\gamma^{*}=O_{p}\left(N^{-\frac{1}{2}}\right)$ and the boundedness of $\nabla_{Y\gamma^{\prime}}\ln f\left(w_{i}|\overline{\gamma}\right)$, the first term is $O_{p}\left(N^{-\frac{1}{2}}\right)$. The second term is also $O_{p}\left(N^{-\frac{1}{2}}\right)$ by the central limit theorem and the assumption that $Var\left(\nabla_{Y}\left(\ln f\left(w_{i}|\gamma^{*}\right)\right)\right)$ is finite. The third term is 0 if the distribution $f$ is correctly specified. Consequently, $\nabla_{Y}Q\left(\widehat{\theta},\widehat{Y}\right)=O_{p}\left(N^{-\frac{1}{2}}\right)$ holds if the distribution $f$ is correctly specified.
Note that $E\left[\nabla_{Y}Q\left(\theta^{*},Y^{*}\right)\right]=0$ also holds. This is because of $E\left[\frac{1}{N}\sum_{i=1}^{N}\left(\nabla_{Y}\ln f\left(w_{i}|\gamma^{*}\right)-\nabla_{Y}\ln f\left(w_{i}|\gamma^{*}\right)\right)\right]=0$ and {$E\left[\frac{1}{N}\sum_{i=1}^{N}\nabla_{Y}\left(\ln f\left(w_{i}|\gamma^{*}\right)\right)-\int\left(\nabla_{Y}\ln f\left(w|\gamma^{*}\right)\right)f^{*}\left(w\right)dw\right]=0$}. We can utilize the property in the discussions related to Neyman orthogonality (Section (ref)).
holds. Here, $\overline{m}\left(\widehat{\theta},\widehat{Y}\right)=\frac{1}{N}\sum_{i=1}^{N}m\left(w_{i};\widehat{\theta},\widehat{Y}\right)=O_{p}\left(N^{-\frac{1}{2}}\right)$ holds if the model is correctly specified.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {More precisely, $\frac{1}{N}\sum_{i=1}^{N}m\left(w_{i};\widehat{\gamma}\right)=\frac{1}{N}\sum_{i=1}^{N}m\left(w_{i};\gamma^{*}\right)+\frac{1}{N}\sum_{i=1}^{N}\left(\nabla_{\gamma}m\left(w_{i};\exists\overline{\gamma}\right)\right)\left(\widehat{\gamma}-\gamma^{*}\right)$ holds by the mean value theorem. Here, $\frac{1}{N}\sum_{i=1}^{N}m\left(w_{i};\gamma^{*}\right)=O_{p}\left(N^{-\frac{1}{2}}\right)$ holds under the assumption of finite $Var\left(m\left(w;\gamma^{*}\right)\right)$ and the central limit theorem. In addition, $\widehat{\gamma}-\gamma^{*}=O_{p}\left(N^{-\frac{1}{2}}\right)$ and the boundedness of $\nabla_{\gamma}m\left(w_{i};\overline{\gamma}\right)$ imply $\frac{1}{N}\sum_{i=1}^{N}\left(\nabla_{\gamma}m\left(w_{i};\overline{\gamma}\right)\right)\left(\widehat{\gamma}-\gamma^{*}\right)=O_{p}\left(N^{-\frac{1}{2}}\right)$.} Hence, by the boundedness of $\nabla_{Y}\overline{m}\left(\theta,Y\right)$, $\nabla_{Y}Q\left(\widehat{\theta},\widehat{Y}\right)=O_{p}\left(N^{-\frac{1}{2}}\right)$ holds.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {In contrast to the case of the MLE, $E\left[\nabla_{Y}Q\left(Y^{*},\theta^{*}\right)\right]=0$ may not hold, because $E\left[\frac{1}{N}\sum_{i}m\left(w_{i};\theta^{*},\gamma^{*}\right)\right]=0$ may not necessarily imply $E\left[\nabla_{Y}Q\left(\theta^{*},Y^{*}\right)\right]=2\left(\nabla_{Y}\frac{1}{N}\sum_{i}m\left(w_{i};\theta^{*},Y^{*}\right)\right)^{\prime}\widehat{W}\left(\frac{1}{N}\sum_{i}m\left(w_{i};\theta^{*},\gamma^{*}\right)\right)$.} Note that, in the exactly identified case ($n_{m}=n_{\theta}$), $\overline{m}\left(\widehat{\theta},\widehat{Y}\right)=0$ holds. This suggests $\nabla_{Y}Q\left(\widehat{Y},\widehat{\theta}\right)=0$, and the local convergence speed is quadratic, regardless of the number of samples, in the setting.
By $\nabla_{Y}Q\left(\widehat{Y},\widehat{\theta}\right)=O_{p}\left(N^{-\frac{1}{2}}\right)$ and Proposition (ref), the following statement holds, which is the counterpart of the ones in dearing2024efficient.
$\gamma_{k}-\widehat{\gamma}=O_{p}\left(N^{-1/2}\left\Vert \gamma_{k-1}-\widehat{\gamma}\right\Vert +\left\Vert \gamma_{k-1}-\widehat{\gamma}\right\Vert ^{2}\right)$ implies that $N^{-1/2}\left\Vert \gamma_{k-1}-\widehat{\gamma}\right\Vert $ is mostly negligible when the number of samples $N$ is large, and $\gamma_{k}-\widehat{\gamma}\approx O_{p}\left(\left\Vert \gamma_{k-1}-\widehat{\gamma}\right\Vert ^{2}\right)$. It implies that the local convergence speed of the algorithm is nearly quadratic in large sample.
Note that the statement does not impose the assumption that the initial value $\gamma_{0}$ is statistically consistent, in contrast to the previous studies (kasahara2008pseudo for the NPL (single-agent DDC) using MLE; dearing2024efficient for the EPL using MLE). The new statement is meaningful, because initial consistent estimates of $\gamma$ is not straightforward to obtain for general models.
Besides, Proposition (ref) in Section (ref) indicates that $\widetilde{\theta}\in\Gamma_{seq}$, which are local minima of the solution, satisfies $\gamma_{k}-\widetilde{\gamma}=A\left(\widetilde{\gamma}\right)\left(\nabla_{Y}Q\left(\widetilde{\theta},\widetilde{Y}\right)\right)\left(\nabla_{\theta\gamma^{\prime}}\Upsilon\left(\widetilde{\theta},\widetilde{Y}\right)\right)\left(\gamma_{k-1}-\widetilde{\gamma}\right)+O\left(\left\Vert \gamma_{k-1}-\widetilde{\gamma}\right\Vert ^{2}\right)$. Nevertheless, the discussion above implies that $\nabla_{Y}Q\left(\widetilde{Y},\widetilde{\theta}\right)\approx0$ may not hold even in large samples if $\left(\widetilde{Y},\widetilde{\theta}\right)\neq\left(\widehat{Y},\widehat{\theta}\right)$. Consequently, the modulus of the sequential algorithm iterations near the local optima other than the global optima may not be close to 0, and less likely to converge to the points, especially when the modulus is larger than 1. This is an unique feature of the sequential algorithms with ZJP discussed in the current paper, and the SQP iteration does not possess the property.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {The SQP obtains the solution of the constrained optimization problem by basically applying Newton's method to the first-order conditions of the constrained optimization problem. Consequently, the SQP iteration attains quadratic local convergence even at local optima. See Chapter 18 of nocedal2006numerical for details of the SQP algorithm.} Although not discussed in the previous studies applying sequential algorithms (NPL for single-agent DDC, EPL for dynamic discrete games), the feature can be thought of as one advantage of the sequential algorithms over the SQP.
In some models (e.g., DDC models), previous studies have developed methods to estimate initial $\gamma$ consistently. In cases where initial consistent estimates of $\gamma_{0}$ is available, we can derive stronger results. The followings are the statements, which corresponds to Theorem 1 of dearing2024efficient in the context of the EPL for MLE:
Note that we do not assume the $\sqrt{N}$-consistency of the initial values $\gamma_{0}$ to derive $\sqrt{N}\left(\theta_{k}-\theta^{*}\right)\rightarrow_{d}N\left(0,\Sigma_{\theta}\right)$, unlike dearing2024efficient. Hence, the efficiency result $\sqrt{N}\left(\theta_{k}-\theta^{*}\right)\rightarrow_{d}N\left(0,\Sigma_{\theta}\right)$ for sequential algorithms applies even when kernel-based estimator is used to consistently estimate $\gamma_{0}$.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {kasahara2008pseudo of Proposition 2 showed that $\sqrt{N}$-consistency assumption imposed in aguirregabiria2002swapping is not necessary to derive the efficiency of the NPL estimator in the single-agent DDC model. The current study generalizes the result to models beyond single-agent DDC models.}
The sequential algorithms with ZJP are closely related to Neyman orthogonality, which is recognized in recent studies (e.g., chernozhukov2018double) as important properties for efficient estimation of parameters under the existence of nuisance parameters . The current section briefly clarifies the point.
First, the derivative of the objective function $\widetilde{Q}\left(\theta;\gamma\right)$ used in the sequential algorithm is $\nabla_{\theta}\widetilde{Q}\left(\theta;\gamma\right)$. By ((ref)) and $G\left(Y;\theta\right)=0\Rightarrow\nabla_{\gamma}\Upsilon\left(\theta^{*};\gamma^{*}\right)=0$ under Assumption (ref), $\left(\theta,Y\right)$ with $G\left(Y;\theta\right)=0$ satisfies $\nabla_{\gamma}\left[\nabla_{\theta}\widetilde{Q}\left(\theta,\gamma\right)\right]=\left(\nabla_{Y}Q\left(\theta,Y\right)\right)\left(\nabla_{\theta\gamma^{\prime}}\Upsilon\left(\theta,\gamma\right)\right)$. Hence, $\nabla_{\gamma}\left[\nabla_{\theta}\widetilde{Q}\left(\theta^{*},\gamma^{*}\right)\right]=\left(\nabla_{Y}Q\left(\theta^{*},Y^{*}\right)\right)\left(\nabla_{\theta\gamma^{\prime}}\Upsilon\left(\theta^{*},\gamma^{*}\right)\right)$ holds.
As discussed in Section (ref), $E\left[\nabla_{Y}Q\left(\theta^{*},Y^{*}\right)\right]=0$ holds for MLE, and it implies $\nabla_{\gamma}\left[E\left[\nabla_{\theta}\widetilde{Q}\left(\theta^{*},\gamma^{*}\right)\right]\right]=0$, i.e., $\nabla_{\gamma}\left[E\left[\nabla_{\theta}\widetilde{Q}\left(\theta^{*},\gamma^{*}\right)\right]\right]\left(\gamma-\gamma^{*}\right)=0\ \forall\gamma$.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {The condition implies that $\widetilde{Q}\left(\theta,\gamma\right)$ is an orthogonal loss function (foster2023orthogonal).} Consequently, the score function $\nabla_{\theta}\widetilde{Q}\left(\theta;\gamma\right)$ is Neyman orthogonal concerning the nuisance parameters $\gamma$ in the MLE setting. Note that $\nabla_{\gamma}\left[\nabla_{\theta}\widetilde{Q}\left(\theta^{*},\gamma^{*}\right)\right]\rightarrow_{p}0$ holds for the GMM setting, as discussed in Section (ref), although $\nabla_{\gamma}\left[E\left[\nabla_{\theta}\widetilde{Q}\left(\theta^{*},\gamma^{*}\right)\right]\right]=0$ may not hold in general.
Intuitively, if $\nabla_{\gamma}\left[E\left[\nabla_{\theta}\widetilde{Q}\left(\theta^{*},\gamma^{*}\right)\right]\right]=0$ or $\nabla_{\gamma}\left[\nabla_{\theta}\widetilde{Q}\left(\theta^{*},\gamma^{*}\right)\right]\rightarrow_{p}0$ holds, the first order condition of the optimization problem $\theta_{k+1}=\arg\min_{\theta}\widetilde{Q}\left(\theta;\gamma_{k}\right)\equiv Q\left(\theta,\Upsilon\left(\theta;\gamma_{k}\right)\right)$ is less likely to be affected by the value of the nuisance parameters $\gamma_{k}$ around $\left(\theta^{*},\gamma^{*}\right)$.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Concerning the SLC using $\Upsilon\left(\theta;\gamma_{k}\right)\equiv Y_{k}-\left(\left(\nabla_{Y}G\left(Y_{k};\theta_{k}\right)\right)^{-1}\left(G\left(Y_{k};\theta_{k}\right)\right)\right)-\left(\left(\nabla_{Y}G\left(Y_{k};\theta_{k}\right)\right)^{-1}\left(\nabla_{\theta}G\left(Y_{k};\theta_{k}\right)\right)\right)\left(\theta-\theta_{k}\right)$, we can treat not only $\theta$ and $Y$ but also $\left(\nabla_{Y}G\left(Y_{k};\theta_{k}\right)\right)^{-1}\left(\nabla_{\theta}G\left(Y_{k};\theta_{k}\right)\right)$ as nuisance parameters $\gamma$, and the ZJP $\nabla_{\gamma}\Upsilon\left(\theta,\gamma\right)=0$ at points satisfying $G\left(Y;\theta\right)=0$ holds even under the setting (See also the proof in Appendix (ref)). This property suggests that locally $\theta_{k+1}$ is less likely to be affected by a slight numerical error in $\left(\nabla_{Y}G\left(Y_{k};\theta_{k}\right)\right)^{-1}\left(\nabla_{\theta}G\left(Y_{k};\theta_{k}\right)\right)$, which can be computed by iterative Krylov-based methods whose solution may not exactly coincide with the exact solution of a linear equation.} Hence, we would be able to obtain updated parameters $\theta_{k+1}$ which would be close to the true value, even if the candidate nuisance parameter values $\gamma_{k}$ are not very close to the true values.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {If we apply the SLC or EPL, we require initial consistent values not only for $Y$ but also for $\theta$ to obtain an efficient estimate of $\theta$ after one iteration. Although this requirement may seem restrictive, studies that rely on Neyman orthogonality, such as sawadogo2025efficient and belloni2018high, also require consistent initial values for both the main parameter and the nuisance parameter to construct Neyman-orthogonal score functions. }
One potential concern with sequential algorithms with ZJP is that, although they have a theoretical guarantee of near\nobreakdash-quadratic local convergence in large samples---as shown in Section (ref) using $\nabla_{Y}Q\left(\widehat{\theta},\widehat{Y}\right)=O_{p}\left(N^{-\frac{1}{2}}\right)$---the convergence can be slow when the pseudo\nobreakdash-objective function $\text{\ensuremath{\widetilde{Q}\left(\theta;\gamma\right)}}$ is relatively flat near the solution. Proposition (ref) indicates that the local convergence rate can be far from quadratic when the Hessian matrix $\nabla_{\theta\theta^{\prime}}\widetilde{Q}\left(\widetilde{\theta},\widetilde{\gamma}\right)$ is close to singular in the finite sample environment. In addition, the iteration may not be contractive when far from the solution.
A practical way to mitigate these issues is to employ fixed\nobreakdash-point iteration acceleration methods (e.g., Spectral, Anderson acceleration, SQUAREM). These methods inherit the spirit of Newton\nobreakdash-type algorithms for solving nonlinear equations but do not require computing additional derivatives. Recent studies in economics have shown that such methods can substantially accelerate and stabilize fixed\nobreakdash-point iterations (e.g., Spectral for dynamic game NPL; SQUAREM and Anderson acceleration for static BLP NFXP\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {See aguirregabiria2021imposing, conlon2020best, and fukasawa2024fast.}).\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {As discussed by pollock2021anderson in the case of the Anderson acceleration method for fixed-point iterations, the use of the acceleration method sometimes lead to the convergence of non-contractive fixed-point mappings. The phenomenon is known in the literature of computational physics, and pollock2021anderson analytically investigates the mechanism behind the phenomenon motivated by the previous studies.} Because sequential algorithms with ZJP repeatedly apply a fixed\nobreakdash-point mapping $H$, these acceleration techniques can be incorporated with minimal modification. For details on fixed\nobreakdash-point acceleration methods, see the Supplemental Appendix (ref). Regarding the SLC, although the original algorithm typically performs well, incorporating acceleration can further improve convergence speed, as demonstrated in Section (ref).
To demonstrate the performance of the SLC, this study conducts a set of numerical experiments. We consider two models: a dynamic discrete game with time-varying unobserved heterogeneity, and a dynamic demand model (dynamic BLP). dearing2024efficient showed that the EPL performs well in dynamic games without unobserved heterogeneity, where the mapping $G$ is linear in $\theta$ and the EPL coincides with the SLC. Accordingly, the present study focuses on settings in which this linearity does not hold.
Because explicit computation of Jacobians is often not straightforward in practice due to the complicated model structure or the memory size, the current study shows results based on the Jacobian-free SLC method, and compare the results with the ones for NFXP algorithms using numerical derivatives. All the experiments were run on a laptop computer with the CPU AMD Ryzen 5 6600H 3.30 GHz, 16.0 GB of RAM, Windows 11 64-bit, and MATLAB 2022b.\cprotect\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Replication code for the numerical experiments is available at \url{https://github.com/takeshi-fukasawa/sequential_SLC}.}
In both models, the tolerance level for the inner-loop of the NFXP algorithm is set to 1E-12. Concerning numerical derivative-based (Jacobian-free) algorithms, we use central finite differences. SLC is implemented in conjunction with the Krylov-based method (GMRES) using numerical derivatives, and the NFXP is implemented by computing the gradient of the objective function concerning statistical parameters by numerical derivatives. Concerning the inner loop of the NFXP, we use fixed-point iterations with Anderson acceleration. The maximum number of SLC iterations is set to 50. We perform the estimation using 20 independent datasets. To mitigate the risk of convergence to local optima, we employ five different starting values.
The dynamic game model we consider introduces time-varying unobserved heterogeneity in the model considered in aguirregabiria2007sequential, egesdal2015estimating, and dearing2024efficient. In addition, the model is similar to the one considered in arcidiacono2011conditional.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {arcidiacono2011conditional considered a setting in which firms exit the market forever once they choose to exit. In contrast, the current study consider the setting where firms can re-enter the market even after having exited, as in aguirregabiria2007sequential, egesdal2015estimating, and dearing2024efficient.} The numerical experiments are done by modifying the replication code of dearing2024efficient.
Time is discrete, and indexed by $t=1,2,\cdots$. In each market, there are $|\mathcal{J}|$ firms, and they are indexed by $j\in\{1,2,\cdots,|\mathcal{J}|\}$. Given states $\omega_{t}$ and private information only observed to each firm $\epsilon_{t}^{j}$, each firm simultaneously chooses its action $a_{t}^{j}\in\mathcal{A}=\{0,1,\cdots,|\mathcal{A}|-1\}$.
Agents maximize the expected discounted utility $E\left\{ \sum_{s=t}^{\infty}\beta^{s-t}\left[\overline{u}^{j}\left(\omega_{s},a_{s}^{j},a_{s}^{-j};\theta_{u}\right)+\epsilon_{s}^{j}(a_{s}^{j})\right]\left|\omega_{t},\epsilon_{t}^{j}\right.\right\} $, where $\beta\in(0,1)$ denotes the discount factor. $\theta_{u}$ and $\theta_{f}$ represent parameters concerning utility and state transitions respectively, and let $\theta\equiv(\theta_{u},\theta_{f})$. Under the assumptions of conditional independence, independent private values, finite state space $\left(\omega\in\Omega=\{1,2,\cdots,|\Omega|\}\right)$, and Markov perfect equilibrium, choice-specific value function $v^{j}\in\mathbb{R}^{|\Omega|\times|\mathcal{A}|}$ satisfies the following equation:\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {See Lemma 1 of dearing2024efficient for the proof.}
Here, we define a function $\Phi_{v}:\Theta\times\mathbb{R}^{|\mathcal{J}|\times|\Omega|\times|\mathcal{A}|}\rightarrow\mathbb{R}^{|\mathcal{J}|\times|\Omega|\times|\mathcal{A}|}$. $\Lambda^{j}(\omega,a^{j};v^{j})$ is the choice probability that agent $j$ chooses action $j$ in state $x$, conditional on having choice-specific value function $v^{j}$. The function $S(\cdot)$ is McFadden's surplus function. We also define $\Lambda^{-j}(v^{-j})\equiv\left(\Lambda^{1}(v^{1}),\cdots,\Lambda^{j-1}(v^{j-1}),\Lambda^{j+1}(v^{j+1}),\cdots,\Lambda^{|\mathcal{J}|}(v^{|\mathcal{J}|})\right)$, and:
Choice-specific value function $v$ is a solution of a nonlinear equation $G(\theta,v)\equiv v-\Phi_{v}(\theta,v)=0$ given $\theta$.
We consider the setting where firms make entry/exit decisions $(a^{j}=1,0)$. Concerning $\overline{u}^{j}$, we assume:
where $\omega_{t}\equiv\left(z_{t},\left\{ a_{t-1}^{j}\right\} _{j\in\mathcal{J}}\right)$ and $z_{t}^{(1)}\in\{1,\cdots,|\mathcal{Z}^{(1)}|\}$ denotes the market size. $z_{t}^{(2)}\in\{1,\cdots,|\mathcal{Z}^{(2)}|\}$ denotes another state affecting the market.
As in dearing2024efficient, we choose parameter values $\left|\mathcal{Z}^{(1)}\right|=5,\theta_{FC,j}=-2+0.1\times j,\theta_{RS_{1}}=\theta_{RN}=\theta_{EC}=1.0$, and $\beta=0.95$. We consider a setting in which data from $N=640$ or $160$ markets with $|\mathcal{J}|=3$ firms are available for $T=10$ consecutive periods. Regarding $\theta_{RS_{2}}$, we assume $\theta_{RS_{2}}=1.0$. Concerning $\theta_{RN}$, we assume $\theta_{RN}=4.0$.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Under $\theta_{RN}=4.0$ and without the additional state $z_{t}^{(2)}$, the NPL algorithm is unstable, unless introducing additional strategies (e.g., spectral algorithm with a line search), as shown in aguirregabiria2021imposing and dearing2024efficient. Under the parameter setting, the present study did not find evidence of multiple equilibria.} Regarding the states $z_{t}^{(1)}$ and $z_{t}^{(2)}$, we assume the state transition probabilities $\left(
\right)$ $(i=1,2)$, where $\pi^{(1)}=\pi^{(2)}=0.8$. In the estimation, we assume the values of $\pi^{(1)}$ and $\beta$ are known.
We assume that the state $z_{t}^{(2)}$ is observed by the firms but not observed by econometricians. Concerning $\left|\mathcal{Z}^{(2)}\right|$, we let $\left|\mathcal{Z}^{(2)}\right|=3$. We estimate utility parameters $\theta$ and the state transition parameter $\pi^{(2)}$. In both the SLC and NFXP algorithms, we treat $\left(v,\overline{p_{1}}\right)$ as the nuisance parameters $Y$, where $\overline{p_{1}}$ denotes the vector of the initial distribution of $\omega$. We assume that the initial period is in the stationary environment and satisfies $\overline{p_{1}}=\Phi_{\overline{p_{1}}}\left(\overline{p_{1}}\right)=\Pi\overline{p_{1}}$, where $\Pi$ denotes the transition probability matrix of $\omega$ in the initial period. Concerning the NFXP, we solve for $\left(v,\overline{p_{1}}\right)$ by applying the fixed-point mappings $\Phi_{v}$ and $\Phi_{\overline{p_{1}}}$ with Anderson acceleration.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Supplemental Appendix (ref) additionally shows numerical results of the SLC and the NFXP where the idea of the Endogenous Value Function Iteration (EVFI), which can accelerate convergence, is applied to these algorithms.} To uniquely identify the parameters, we impose an inequality constraint $\theta_{RN_{2}}\geq0$ in the estimation. We also impose $\pi\in[0,1]$. Note that we can easily introduce inequality constraints on the structural parameters, as discussed in the Supplemental Appendix (ref). The parameters $\left(\theta_{u},\pi^{(2)}\right)$ are estimated by the MLE given the constraints on $\left(v,\overline{p_{1}}\right)$. The current study applies the interior-point algorithm in fmincon function of MATLAB for the outer loop of the NFXP.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {The present study also experimented with the Berndt--Hall--Hall--Hausman (BHHH) algorithm (hall1974estimation) as an outer\nobreakdash-loop optimization method for the NFXP, whose effectiveness is noted in iskhakov2016comment. However, in the current setting, the algorithm occasionally moved to parameter values for which the inner\nobreakdash-loop iterations failed to converge. The current study also experimented with the commercial solver KNITRO, but it did not yield substantial speed improvements.} The tolerance level for the convergence of the SLC and the outer loop of the NFXP is set to 1E-6. Supplemental Appendix (ref) describes the procedure for computing the likelihood; see that section for further details.
Table (ref) shows the computational performance, and Table (ref) shows the estimated parameters. Table (ref) suggests that the SLC is roughly 4-8 times faster than the NFXP, regardless of the use of analytical derivatives and the number of independent markets $N$. Such a speedup mainly comes from the fewer number of main (outer-loop) iterations. Note that the last column of Table (ref) indicates that the differences between the SLC and NFXP estimator values are mostly negligible.
Furthermore, the SLC converges in most cases even absent any stabilization strategies, suggesting that the SLC is practical as the solution method for estimating dynamic games even with time-varying unobserved heterogeneity.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {dearing2024efficient showed numerical results where the convergence of the EPL, which is equivalent to the SLC under the linearity of $G$ concerning $\theta$, is stable even when starting from random initial values in a dynamic game model without unobserved heterogeneity.}
We consider a dynamic BLP model where forward-looking decision-makings of consumers are explicitly incorporated into the standard static BLP model (berry1995automobile). Here, we focus on the model of perfectly durable goods considered in sun2019computationally.
Consumer $i$'s utility when purchasing product $j$ at time $t$ is:
\[ U_{ijt}=\chi_{jt}\theta_{i}+\xi_{jt}+\epsilon_{ijt} \] where $\chi_{jt}\equiv\left(1,\chi_{1jt},\chi_{2jt},\chi_{3jt},p_{jt}\right)$ and $\xi_{jt}$ denote the observed and unobserved product characteristics of product $j$ at time $t$.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Here, $p_{jt}$ denotes the product price.} The coefficient $\theta_{i}\sim N(\overline{\theta},\Sigma)$ represents random coefficients of consumer $i$.
Consumer $i$'s utility when purchasing nothing is:
\[ U_{i0t}=\beta E_{t}[V_{it+1}(\Omega_{t+1})|\Omega_{t}]+\epsilon_{i0t}, \] where $\beta$ denotes the consumers' discount factor, and $\epsilon$ denotes idiosyncratic utility shocks following Gumbel distribution. The variable $V_{it}$ denotes consumer $i$'s (integrated) value function at time $t$, and $\Omega$ denotes state variables.
We estimate demand parameters $(\overline{\theta},\Sigma)$ by GMM by imposing moment conditions $E[\xi|Z]=0$, where $Z$ denotes instrumental variables. We introduce consumer heterogeneity in preferences for the characteristics $\chi_{1jt},\chi_{2jt},p_{jt}$, and let $\sigma_{\chi_{1}},\sigma_{\chi_{2}},\sigma_{p}$ denote the standard deviation parameters of the normal distribution governing the corresponding random coefficients. We further assume that $S_{jt}^{(data)}=s_{jt}\left(\xi,V,\theta\right)\ \forall j,t$, namely, observed market shares are equal to predicted market shares. Concerning consumer expectations, we assume inclusive value sufficiency, as in the previous studies (e.g., gowrisankaran2012dynamics). We use 50 grid points to interpolate the values of inclusive values. We consider a setting where $J=25$ products exist in a market, and the market data is available for $T=25$ consecutive periods. For details of the model setting, see Section 5.2 of sun2019computationally or fukasawa2024fast. Supplemental Appendix (ref) also describes the details.
For the inner-loop of the NFXP, we use the computationally efficient fixed-point mapping with the Anderson acceleration method developed in fukasawa2024fast. fukasawa2024fast considers a fixed-point mapping of the value functions $V$, which does not depend on the mean product utilities, and shows that the new inner-loop algorithm with Anderson acceleration is roughly ten times faster than the traditional fixed-point iteration algorithm that combines the BLP contraction mapping and value function iterations.
In the present application, we rely on the numerical derivative-based Jacobian-free approach. In dynamic BLP models with random coefficients, deriving Jacobian of function $G$ is not straightforward due to the complexity of the model structure.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {See also Supplemental Appendix (ref) for potential limitations of automatic differentiations.} In addition, the Jacobian matrix is typically high-dimensional, implying that a substantial amount of computer memory is required merely to compute and store it.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {In the current setting where the number of discretized consumer types is 50, the Jacobian is a 3750$\times$3750 dimensional matrix, which consumes roughly 107MB of memory. If we allow for more flexible consumer heterogeneity and let the number of discretized consumer types be 1,000, which is typical in the static BLP applications, the Jacobian is a 75000$\times$75000 dimensional matrix, and it requires roughly 43GB of memory. Although the present model includes only two consumer-level state variables (owning a durable good or not), the dimensionality would increase further as additional states are introduced. It is worth noting that, for the MPEC approach (interior-point methods that explicitly compute Hessians), sun2019computationally reported in Appendix C of their paper that more than 32 GB of memory is required even when the number of discretized consumer types is only 50.}
Table (ref) shows the computational performance, and Table (ref) shows results on estimated parameter values. Results suggest that the original SLC usually performs well, and it is roughly 5 times faster than the NFXP. Such a large speedup mainly arises from from the much smaller number of the evaluations of the function $G$. the number of $G$ evaluations required by the SLC is roughly five times smaller than that required by the NFXP. Although the SLC requires a larger number of evaluations of the objective function $Q$, code profiling indicates that the computational cost of evaluating $Q$ is roughly 10 times smaller than that of evaluating $G$. Consequently, the substantially fewer evaluations of $G$ dominate the overall computational savings. Note that the SLC and the NFXP estimators take mostly the same values, as shown in the latter table.
It is worth noting that the original SLC is occasionally slow, and sometimes does not converge within 50 iterations.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {Occasional non\nobreakdash-convergence is also not uncommon even under the “MPEC” approach (su2012constrained; egesdal2015estimating).} As discussed in Section (ref), sequential algorithms with ZJP may require many iterations when the objective function becomes nearly flat around the solution. In addition, local convergence in large samples does not necessarily imply convergence from far initial points. To mitigate the issues, the current study also experimented the SLC with the spectral algorithm, which is one of the fixed-point acceleration methods and utilized in the NPL algorithm for dynamic games (aguirregabiria2021imposing). The results suggest that the spectral algorithm significantly mitigates the problem and enhancing the robustness of the SLC--- SLC-Spectral is roughly 7 times faster than the NFXP, and always converges.\ifx\@typeset@protect \expandafter\SF@@footnote \else \expandafter\SF@gobble@opt \fi {In the current setting, the SLC-Spectral and the NFXP estimators take different values in 5 out of the 20 randomly generated datasets, although the absolute value difference is at most less than 0.05. In the current numerical experiment, the standard deviation parameters of the random coefficients are not restricted to be positive, which may increase the number of potential local minima. Supplemental Appendix (ref) further discusses the issue.}
The current study has investigated sequential algorithms possessing the zero Jacobian property (ZJP) for estimating structural models with equilibrium constraints, which can also be regarded as numerical algorithms for solving constrained optimization problems. By developing a more general framework than those in previous studies, the study shows that they achieve near-quadratic local convergence in large samples, even without consistent initial estimates, in the context of MLE or GMM.
It then proposes a novel algorithm termed Sequential Linearly Constrained (SLC) algorithm. It is easily applicable to a broader class of structural models. A key advantage of the SLC is that it can be implemented without explicitly computing the Jacobian of the equilibrium constraints, and it can be several times faster than the NFXP approach. This study demonstrates its good performance through two numerical experiments: dynamic discrete game with time-varying unobserved heterogeneity and dynamic BLP model.
While the SLC algorithm performs well even without initial consistent estimates of $\theta$ and $Y$, incorporating such estimates is desirable as the algorithm can yield statistically efficient parameter estimates even after one iteration. Computationally efficient methods for obtaining consistent estimates exist for specific economic models---such as the CCP-based approach for dynamic discrete choice models with finite dependence (arcidiacono2011conditional)---, and such methods can be combined with the SLC to obtain statistically efficient parameter estimates. However, whether such computationally light methods exist for a broader class of models remains an open question. Developing these general methods with small computational costs for consistent parameter estimation remains a subject for future research.