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.
100,493 characters · 14 sections · 47 citation commands
Convexity Not Required: Estimation of Smooth Moment Condition Models
JEL Classification: C11, C12, C13, C32, C36.\newline Keywords: Non-linear estimation, over-identification, misspecification, nonlinear systems of equations, injectivity, local and global identification.
\baselineskip=18.0pt \thispagestyle{empty} \setcounter{page}{0}
The Generalized and Simulated Method of Moments (GMM, SMM) are commonly used to estimate structural Economic models. To find estimates, modern computer software provides researchers with a large set of free and non-free numerical optimizers, which, after inputting some tuning parameters, return a guess for the parameters of interest. While sampling properties of estimators are often derived, their practical implementation often receives a less detailed treatment. There is now a vast literature on statistical learning with a convex loss function, using stochastic gradient-descent. However, these results need not directly apply to GMM, as it often involves non-convex minimizations. A number of authors have pointed out the lack of robustness of off-the-shelf methods, and knittel2014 illustrate this in the context of demand estimation. This is perhaps not surprising since non-convex optimization is subject to a curse of dimensionality andrews1997 and becomes increasingly challenging when the number of parameters is moderate or large.
The main contribution of the paper is to show that convexity is not required for some methods to perform well in GMM estimation specifically: some algorithms are globally convergent under a global rank condition involving the Jacobian of the moments and the weighting matrix. This defines a class of non-convex problems that is as hard as convex problems for optimization. Since this is perhaps surprising, the following gives some intuition behind the result. Given sample moments $\overline{g}_n(\theta)$ with Jacobian $G_n(\theta)$, one can minimize the GMM objective function $Q_n(\theta) = 1/2 \overline{g}_n(\theta)^\prime W_n \overline{g}_n(\theta)$ iteratively, by minimizing successive quadratic approximations. To this end, convex optimizers rely on a quadratic expansion of $Q_n(\theta)$ using its gradient and Hessian. This quadratic approximation yields a proper minimization problem only if the Hessian is strictly positive definite, i.e. $Q_n$ is convex.
Another approach, is to linearly expand the sample moments using the Jacobian and plug the linearized moments into the GMM objective. Since the approximate moments are linear, this yields a proper minimization problem as long as $G_n$ has full rank. Gauss-Newton (gn) relies on this approach. Gradient-descent can be motivated by either the quadratic of linear approximation. In the just-identified case, it is well known that gn is locally convergent when $G_n$ has full rank around the solution. This paper goes further by showing that gn and gradient-descent are globally convergent when the product of $G_n$, $W_n$, and an average of $G_n$ has full rank everywhere. Unlike existing results, this applies to just and over-identified moments. The condition can be relaxed for the product to only be non-singular in a specific direction, towards the global minimizer. Under this weaker condition, gn with a Levenberg-Marquardt regularization and gradient descent are globally convergent. Importantly, for correctly specified models, the conditions imply that there are no local optima, besides the global minimizer; a necessary condition for global convergence of gradient-based optimizers. It is shown that these convergence results are robust to 1) moderate misspecification, and 2) moderately non-linear reparameterizations. However, the results may or may not hold depending on the choice of weighting matrix. In particular, when $W_n$ is ill-conditioned, convergence can be significantly slower.
Several conditions found in the convex and non-convex optimization literatures imply the weaker condition introduced of this paper. These include strong, star, and quasar convexity of the objective function. It also relates to the Polyak-Łojasiewicz condition, an important inequality which has gathered much interest in machine learning to prove convergence of gradient-descent. Strong monotonicity of the moments, a condition for solving just-determined system of non-linear equations, and strong injectivity, introduced here for just and over-identified models, also imply the weaker condition. Hence, the condition introduced in this paper is a common denominator of several existing conditions. In terms of econometrics properties, the conditions are sufficient for the parameters to be both locally and globally identified, when the model is correct or moderately misspecified.
A simple MA(1) estimation from gourieroux1996 illustrates the results analytically and numerically. The problem is non-convex: the scalar Hessian can be positive, negative, or zero; yet the conditions hold. As predicted, the recommended gn algorithm converges. Newton-Raphson provably diverges, and off-the-shelf optimizers can be unstable. When the model is moderately misspecified, gn remains globally convergent. In line with theory, significant misspecification can produce non-global optima which hinder the global convergence of gradient-descent and Gauss-Newton.
Two empirical applications further illustrate the results. The first application revisits the numerical results of knittel2014 for estimating random coefficient demand models on Nevo's generated cereal data. The same gn algorithm systematically converges from a wide range of starting values. In contrast, R's more sophisticated built-in optimizers can be inaccurate and often crash without additional error-handling. The second application estimates a small New Keynesian model with endogenous total factor productivity by impulse response matching. Matlab's built-in optimizers have better error-handling so that crashes are less problematic. Nonetheless, these optimizers' performance can be mixed whereas gn performs well for nearly all starting values.
Numerically, in all three applications, the GMM objective is non-convex at most values. The strong injectivity condition holds at most values, an indication that gn and gradient-descent are appropriate. The later converges very slowly, however. These findings explain the good performance of gn relative to more commonly used methods. The main takeaway is that non-convexity need not be a deterrent to structural estimation: simple algorithms can converge quickly and globally under alternative conditions.
\paragraph{Structure of the paper.} Section (ref) contains the main assumptions and results. Section (ref) reviews existing conditions found in the literature and relates the main assumptions with these conditions. Section (ref) suggests a numerical procedure to check whether the main assumption holds or not and a way to set the tuning parameter. Section (ref) illustrates the results with one numerical and two empirical applications. Appendices (ref) and (ref) give the proofs to the main results and additional results. The Supplemental Material consists of: Appendices (ref)-(ref). Appendix (ref) provides additional local convergence results, which complement the main global convergence results in the paper. Appendix (ref) provides of survey of empirical practice in the American Economic Review between 2016 and 2018. Appendix (ref) gives R code to replicate the numerical MA(1) example. Appendix (ref) provides additional simulation and empirical results. Appendix (ref) gives additional details about the methods found in the survey of Appendix (ref).
\paragraph{Notation:} In the following $\lambda_{\min}$, $\lambda_{\max}$ return the smallest and largest eigenvalues of a square positive semidefinite matrix. For an arbitrary rectangular matrix $A$ of size $d_g \times d_\theta$ with $d_g \geq d_\theta$, $\sigma_{\min}$, $\sigma_{\max}$ are the smallest and largest singular values of $A$ defined as $\sigma_{\min}(A) = \sqrt{ \lambda_{\min}(A^\prime A) }$ and $\sigma_{\max}(A) = \sqrt{ \lambda_{\max}(A^\prime A) }$; $A$ has full rank if, and only if, $\sigma_{\min}(A) > 0$.
Let $\overline{g}_n(\theta) = 1/n \sum_{i=1}^n g(\theta;x_i)$ be the sample moments and $G_n(\theta) = \partial_\theta \overline{g}_n(\theta)$ their Jacobian. Their population counterparts are $g(\theta) = \mathbb{E}[g(\theta;x_i)]$ and $G(\theta) = \partial_\theta g(\theta)$. $W_n$ is a weighting matrix which, for simplicity, does not depend on $\theta$ -- this excludes continuously-updated estimations. The sample GMM objective function is: \[ Q_n(\theta) = \frac{1}{2}\overline{g}_n(\theta)^\prime W_n \overline{g}_n(\theta), \] and the goal is to find the global minimizer $\hat{\theta}_n$ of $Q_n$ in $\mathbb{R}^{d_\theta}$. The population objective $Q(\theta) = \frac{1}{2} g(\theta)^\prime W g(\theta)$, defined similarly using the limit $W$ of $W_n$, has a global minimizer $\theta^\dagger$. Throughout, it will be assumed that the sample $Q_n$ is continuously differentiable. More specifically, this paper considers derivative-based optimizers of the form:
for $k = 0,1,\dots$, some staring value $\theta_0 \in \mathbb{R}^{d_\theta}$ and a matrix $P_{k,n}$, called conditioning matrix, assumed to be symmetric. The tuning parameter $\gamma \in (0,1]$ is called the learning rate. There are several ways to motivate ((ref)) as a minimization algorithm in the context of GMM estimation. They are conceptually similar but implicitly rely on a different set of assumptions. The first is to consider a quadratic approximation of the GMM objective function $Q_n$: \[ Q_n(\theta) \simeq Q_n(\theta_k) + \partial_\theta Q_n(\theta_k)(\theta - \theta_k) + \frac{1}{2 \gamma} (\theta - \theta_k)^\prime \partial^2_{\theta,\theta^\prime} Q_n(\theta_k) (\theta - \theta_k), \] here $\gamma$ penalizes the quality of the quadratic approximation. For linear models, such as OLS and IV regressions, $Q_n$ is quadratic so that $\gamma = 1$ is feasible. For non-linear models, the approximation is inexact, and $\gamma < 1$ is generally required. Minimizing the right-hand-side with respect to $\theta$ yields a Newton-Raphson (nr) iteration: $\theta_{k+1} = \theta_k - \gamma [\partial^2_{\theta,\theta^\prime} Q_n(\theta_k)]^{-1}\partial_\theta Q_n(\theta_k)$ with $ \partial_\theta Q_n(\theta_k) = G_n(\theta_k)^\prime W_n \overline{g}_n(\theta_k)$ and $P_{k,n} = [\partial^2_{\theta,\theta^\prime} Q_n(\theta_k)]^{-1}$. A quasi-Newton (qn) iterations replaces the Hessian matrix $\partial^2_{\theta,\theta^\prime} Q_n(\theta_k)$ with an approximation computed sequentially over $k$. The most popular qn software implementation is called bfgs. Importantly, the quadratic approximation implicitly requires that is $H_n(\theta) = \partial^2_{\theta,\theta^\prime} Q_n(\theta_k)$ strictly positive definite around $\theta_k$ so that ((ref)) yields a minimizer of the quadratic approximation.
Another way to motivate ((ref)) is to consider a linear approximation of the moments and plug it into the GMM objective function:
where now $\gamma$ penalizes the quality of the linear approximation. Take the first order condition in the last display to find ((ref)) with $P_{k,n} = (G_n(\theta_k)^\prime W_n G_n(\theta_k))^{-1}$, a Gauss-Newton (gn) iteration. The quadratic approximation requires the Hessian $H_n$ of $Q_n$ to be strictly positive definite at $\theta_k$. A gn iteration minimizes the linear approximation as long as the Jacobian $G_n$ of $\overline{g}_n$ has full rank at $\theta_k$ so that $G_n(\theta_k)^\prime W_n G_n(\theta_k)$ is strictly positive definite. Standard regularity condition imply local convexity around $\hat{\theta}_n$. Still, convexity is more challenging to satisfy away from the solution since $\|\overline{g}_n(\theta_k)\| \gg 0$ can result in a non-definite Hessian $H_n(\theta_k) = G_n(\theta_k)^\prime W_n G_n(\theta_k) + (\overline{g}_n(\theta_k)^\prime W_n \otimes I_d) \partial_\theta \text{vec}[ G_n(\theta_k)^\prime ]$, depending on the last term. This suggests that quadratic-based methods (nr, bfgs) and linear-based methods (gn) can behave differently when $Q_n$ is globally non-convex. Gradient-Descent (gd) can be motivated by either a linear or a quadratic approximation. The following summarizes the choice of $P_{k,n}$ for each algorithm:
The following gives the main assumptions on the population moments used to describe the large sample properties of the estimator $\hat{\theta}_n$ and optimization algorithms.
Assumption (ref) consists mainly of standard conditions to derive asymptotic properties for $\hat{\theta}_n$. The iid assumption can be relaxed to allow for time-series dependence. The parameter space is unbounded to accommodate the unconstrained optimization. The technical condition (v) and the next Assumption imply that $Q_n$ has a strictly quadratic lower bound. This ensures consistency without assuming compactness or uniform consistency of the sample moments. The quantity $\sigma_{\min}[G(\theta)]$ in the local identification condition refers to the smallest singular value of $G(\theta)$. The main Assumption (ref) below will rely on the following quantities: \[ \overline{G}(\theta) = \int_0^1 G(\omega \theta + (1-\omega)\theta^\dagger) d\omega, \quad \overline{G}(\theta_1,\theta_2) = \int_0^1 G(\omega \theta_1 + (1-\omega)\theta_2) d\omega. \] The matrix $\overline{G}(\theta)$ is an average derivative over the path from $\theta$ to the solution $\theta^\dagger$. The matrix plays a role in the mean-value identity: $g(\theta_1)-g(\theta_2) = \overline{G}(\theta_1,\theta_2)(\theta_1-\theta_2)$ (see Lemma (ref)).
Assumption (ref) gives the main conditions used in this paper for global GMM estimation of just and over-identified models.\footnote{The factor $\rho$ is assumed to be set, without loss of generality, such that $\sigma_{\min}[\overline{G}(\theta)] > \underline{\sigma}$ under (a) and $\|\overline{G}(\theta)(\theta-\theta^\dagger)\| > \underline{\sigma} \|\theta-\theta^\dagger\|$ under (b) for $\underline{\sigma}$ found in Assumption (ref).} Assumption (ref) (a) replaces the convexity condition $0 < \underline{\lambda}_H \leq \lambda_{\min}[H_n(\theta)] \leq \lambda_{\max}[H_n(\theta)] < \overline{\lambda}_H < \infty$ used to derive convergence results for gd, nr and qn.\footnote{See Nesterov2018, especially equations (1.2.25), (1.2.27) and Theorem 1.2.4 for gd.}, which may not hold for GMM. A sufficient, but restrictive, condition for Assumption (ref) (a) is that $g$ is the derivative of a convex function, for instance a Probit log-likelihood function. Further sufficient conditions are listed in Section (ref). Assumption (ref) (a) implies Assumption (ref) (b); the latter is the weaker condition. Assumption (ref) (a) implies that $G(\theta)$ has full rank for all $\theta$, Assumption (ref) (b) only requires $G(\theta)^\prime W \overline{G}(\theta)$ to be non-singular in the relevant direction $(\theta-\theta^\dagger)$. For over-identified models, both conditions (a) and (b) depend on the choice of weighting matrix $W$. Indeed, unlike square matrices, the product of full rank rectangular matrices does not automatically have full rank,\footnote{Take $G(\theta_1)^\prime = (1, 0)$ and $G(\theta_2)^\prime = (0, 1)$, both have full rank and yet $G(\theta_1)^\prime G(\theta_2)=0$ is singular.} and the weighting matrix changes the way $G$ and $\overline{G}$ are multiplied. It is possible for the product to be singular even when $G$ and $\overline{G}$ have full rank. Importantly, Assumption (ref) may or may not hold depending on the choice of weighting matrix $W$. If Assumption (ref) is not satisfied using the preferred weighting matrix, the algorithm remains locally convergent. A two-step estimation, with a weighting matrix for which Assumption (ref) holds in the first step, would provide a valid estimation strategy in that case. Assumption (ref) is invariant to some one-to-one reparameterizations, this is shown in the next section.
Under Assumption (ref), the parameters are both locally and globally identified (i.e. Assumption (ref) (iv) and (i)). Conversely, Assumption (ref) (iv) implies that Assumption (ref) (a) holds locally around $\theta^\dagger$. The condition requires that it holds globally rather than locally.\footnote{See Lemmas (ref), (ref) and Propositions (ref), (ref).} Under Assumptions (ref) and (ref), a sample analog of Assumption (ref) holds for the following quantities: \[ \overline{G}_n(\theta) = \int_0^1 G_n(\omega \theta + (1-\omega)\hat{\theta}_n) d\omega, \quad \overline{G}_n(\theta_1,\theta_2) = \int_0^1 G_n(\omega \theta_1 + (1-\omega)\theta_2) d\omega, \] with probability approaching 1, this is shown in Lemma (ref). When Assumption (ref) cannot be verified analytically, a related condition which does not involve the minimizer can be checked numerically on the sample moments and their Jacobian. This is considered in Section (ref).
Lemma (ref) shows that, although the parameter space is unbounded and $Q_n$ is non-convex, $\hat{\theta}_n$ is a consistent estimator under Assumptions (ref) and (ref).
Assumption (ref) requires $P_{k,n}$ to be finite and strictly positive definite. This is always the case for gd since $P_{k,n} = I_d$, and holds for gn under Assumption (ref) (a). If the moments only satisfy Assumption (ref) (b), Assumption (ref) does not necessarily hold for gn since the Jacobian $G_n(\theta_k)$ can be singular, but it remains valid for gd. When Assumption (ref) fails, one approach is to regularize the inverse using the so-called Levenberg-Marquardt (LM) algorithm to gn by setting $P_{k,n} = (G_n(\theta_k)^\prime W_n G_n(\theta_k) + \lambda I_d )^{-1}$ so that $\overline{\lambda}_P < \lambda^{-1} < \infty$ and $P_{k,n}$ is finite. Note that Assumption (ref) does hold for gn under strong injectivity conditions introduced in the next Section. nocedal-wright:06 list several additional approaches to enforce Assumption (ref), mainly for convex optimizers.
The following provides the main results: the global convergence properties of gradient-based algorithms. In the following, the initial value $\theta_0$ is taken from $\Theta$, a compact subset of $\mathbb{R}^{d_\theta}$. This is a technical assumption; although the optimization is unconstrained, the sample moments are not uniformly consistent on $\mathbb{R}^{d_\theta}$ which complicates the analysis. The following shows global convergence, uniformly over $\theta_0 \in \Theta$. The main idea is to show that, with probability approaching $1$, the optimization path $(\theta_k)_{k\geq 0}$ is restricted to a compact set, determined by $\theta_0$, where the sample moments are uniformly consistent. Without loss of generality, $\Theta$ is assumed convex and large enough that $\theta^\dagger \in \text{interior}(\Theta)$. In addition, local convergence results can be found in Appendix (ref), those results are new in the case of overidentified and misspecified models as they allow for $\|\overline{g}_n(\hat{\theta}_n)\|_{W_n}\neq 0$.
Theorem (ref) provides global convergence results that are comparable to the convex case. Because the factor $(1-\overline{\gamma})$ is less than $1$, the distance to the solution $\|\theta_{k+1}-\hat{\theta}_n\|$ decreases exponentially fast with $k$, as in the convex case. Several factors affect convergence. The constants $\underline{\lambda}$, $\overline{\lambda}$ coincide with $C_2 = 1/2 \rho^2 \underline{\sigma}^2/[\overline{\sigma}^2 \overline{\lambda}_W]$, $C_3 = 1/2 \overline{\sigma}^2 \overline{\lambda}_W$ in Proposition (ref) below. The convergence rate $1-\overline{\gamma}$ depends on $C_1 = 1/2 \rho^2 \underline{\sigma}^2 / [\overline{\sigma}^2 \overline{\lambda}_W]$, from the same Proposition.
Through these constants, it appears that identification strength - here measured by $\rho \underline{\sigma}$ - and the choice of weighting matrix $W_n$ affect the convergence properties. In particular, a weighting matrix that is ill-conditioned can lead to slower convergence. This can make optimization challenging. When the sample moments are highly correlated, the optimal weighting matrix can be ill-conditioned. Using equal weighting, a diagonal weighting matrix, or regularizing the optimal weighting matrix with $W_n = ( \hat{V}_n + \lambda I_d )^{-1}$, where $\hat{V}_n$ estimates the variance of $\sqrt{n}\overline{g}_n(\theta^\dagger)$, could improve numerical stability.
The size of $\|\overline{g}_n(\hat{\theta}_n)\|_{W_n}$ further affects convergence. The constant $C$ coincides with $C_4 = \overline{\lambda}_W^{1/2} L$ in Proposition (ref) below. The constant $L$ measures the non-linearity of the sample moments, $L=0$ corresponds to linear models. For linear models, $C=0$ implies that $\|\overline{g}_n(\hat{\theta}_n)\|_{W_n}$ does not affect convergence. Non-linear models have $L>0$ which makes optimization more sensitive to $\|\overline{g}_n(\hat{\theta}_n)\|_{W_n}$ for overidentified models.
In applications, $\|\overline{g}_n(\hat\theta_n)\|_{W_n}$ can be relatively large so that misspecification becomes a concern. Understanding the robustness of Theorem (ref) to non-negligible deviations from $Q(\theta^\dagger)=0$ is then empirically relevant. The following considers models where the quantity: \[ Q_n(\hat{\theta}_n) \overset{p}{\to} Q(\theta^\dagger) := \varphi/2 >0 \] does not vanish asymptotically which implies that $\|\overline{g}_n(\hat{\theta}_n)\|_{W_n}$ matters for convergence, even in large samples. Since $G_n$ cannot be full rank at $\theta=\hat\theta_n$ when the model is both just-identified and misspecified, the results presented here solely consider over-identified models.\footnote{The solution $\hat\theta_n$ is s.t. $G_n(\hat\theta_n)^\prime W_n \overline{g}_n(\hat\theta_n)=0$, misspecification implies $\overline{g}_n(\hat\theta_n) \neq 0$, and since $W_n$ has full rank, it must be that $G_n(\hat\theta_n)$ is singular for just-identified models. For over-identified models, $\overline{g}_n(\hat\theta_n)$ is in the null space of $G_n(\hat\theta_n)^\prime W_n$, which allows $G_n(\hat\theta_n)$ to be full rank.}
Theorem (ref) shows that convergence is robust to `moderate' amounts of misspecification. For linear models, $L=0$ implies that ((ref)) reads $\varphi < +\infty$, which is not restrictive. In ((ref)), the choice of $W_n$, nonlinearity, and identification strength restrict the amount of misspecification allowed in ((ref)). The restrictions ((ref)) are discussed further with Proposition (ref) below. The convergence rate $1-\overline{\gamma}$ also depends on $\varphi$, which slows convergence. In the limit, its expression is given by $(1-\overline{\gamma})^2 = 1 - \gamma \underline{\lambda}_P C_1 / 2$, where $C_1 = (\rho\underline{\sigma} - \overline{\lambda}_W^{1/2}L \sqrt{\varphi})^2/[C_3 + C_4 \sqrt{\varphi}]$. The constants $C_3$, $C_4$ appear in Proposition (ref) below. The first of the two terms in the upper bound in ((ref)) ensures that $\overline{\gamma}>0$ is feasible. Having $\varphi \neq 0$ makes convergence slower and estimation more challenging. When $\varphi$ is arbitrarily large, global convergence can fail. This is explained in the next Section, and illustrated with an MA(1) example. Since the magnitude of $\varphi$ depends on the choice of moments $\overline{g}_n$ and weighting matrix $W_n$, a careful selection of these two might mitigate this issue.
\paragraph{Convexity, monotonicity and the Polyak-Łojasiewicz condition.}
The following briefly reviews some convexity conditions found in the literature and an important relaxation called the Polyak-Łojasiewicz (PL) condition. The latter has gathered much attention in the machine learning literature in recent years. Because Assumption (ref) is stated on population quantities, the following discussion will focus on $Q$.
For general minimization of an objective $Q$, gd, nr and qn are globally convergent for $\theta^\dagger$ if $Q$ is $\mu$-strongly convex, i.e. if for some $\mu >0$: \[ Q(\theta_2) \geq Q(\theta_1) + \partial_\theta Q(\theta_1)(\theta_2-\theta_1) + \frac{\mu}{2}\|\theta_1-\theta_2\|^2, \] for all $\theta_1,\theta_2 \in \mathbb{R}^{d_\theta}$. When $Q$ is twice continuously differentiable it is strongly convex if its Hessian $H(\theta) = \partial^2_{\theta,\theta^\prime} Q(\theta)$ is strictly positive definite everywhere with $0 < \underline{\lambda}_H < \lambda_{\min}[H(\theta)] \leq \lambda_{\max}[H(\theta)] < \overline{\lambda}_H < \infty$. Under strong convexity, for $\gamma > 0$ sufficiently small and any $\theta_0$: \[ Q(\theta_{k+1}) - Q(\theta^\dagger) \leq (1-\eta) \left( Q(\theta_{k}) - Q(\theta^\dagger) \right), \] for some $\eta \in (0,1)$ which depends on $\gamma$, the choice of algorithm, i.e. $P_{k,n}$, and the eigenvalues of $H$. Iterating on this inequality indicates that the fit improves rapidly from any starting value $\theta_0$: $Q(\theta_{k}) - Q(\theta^\dagger) \leq (1-\eta)^{k} \left( Q(\theta_{0}) - Q(\theta^\dagger) \right)$. Under strong convexity, $Q$ has a unique global minimizer and no local optima. The literature has considered a number of relaxations of strong convexity under which gd is globally convergent. This includes the so-called star convexity condition introduced by nesterov2006: \[ Q(\theta^\dagger) \geq Q(\theta) + \lambda \partial_\theta Q(\theta)(\theta^\dagger-\theta) + \frac{\mu}{2}\|\theta-\theta^\dagger\|^2\] for some $\mu \geq 0$ and $\lambda = 1$. Fast convergence results for $\theta$ require $\mu > 0$. This is similar-looking to strong convexity but only involves the pairs $(\theta_1,\theta_2) = (\theta,\theta^\dagger)$. For these functions, the convexity property only holds on line segments toward $\theta^\dagger$. Star convexity implies that $\theta^\dagger$ is the unique global minimizer of $Q$. This condition can be further weakened to \textit{quasar convexity}, which allows for $\lambda > 1$ in the inequality above. hinder2020, Figure 1, plot several functions that satisfy these conditions.
karimi2016, guminov2017 showed that a number of relaxations of strong convexity imply the so-called Polyak-{\L}ojasiewicz (PL) inequality, named after polyak1963 and lojasiewicz1963, which requires that:
for all $\theta \in \mathbb{R}^{d_\theta}$ and some $\mu >0$. When $Q$ satisfies the PL inequality, $\partial_\theta Q(\theta) = 0$ implies $\theta$ is globally optimal, i.e. $Q(\theta) = Q(\theta^\dagger)$. The arg-minimizer may not be unique, however, unlike strong convexity. If the PL inequality holds and $\partial_\theta Q$ is Lipschitz continuous, it can be shown that for $\gamma >0$ small enough: $Q(\theta_{k+1}) - Q(\theta^\dagger) \leq (1-\eta) \left( Q(\theta_{k}) - Q(\theta^\dagger) \right)$ for gd karimi2016. This does not imply that $\theta_{k+1}$ converges to $\theta^\dagger$, however, unless the arg-minimizer is unique. Because strong convexity implies the PL inequality, karimi2016 argue that the latter holds locally over a larger area than strong convexity, predicting better optimization performance. They also note that it is difficult to characterize which functions satisfy the PL inequality. They show that $Q(\theta) = h(A\theta)$, with $h$ strongly convex and $A$ a non-zero matrix, satisfies the PL inequality.
Closely related to the GMM setting, a smaller literature has considered conditions for solving non-linear systems of equations of the form: $g(\theta) = 0$, typically with $g$ and $\theta$ of the same dimension. An important reference is dennis1996, who cast the problem as minimizing $Q(\theta) = \|g(\theta)\|^2$, similar to GMM, and derive global convergence results to a local minimum under convexity conditions (Theorems 6.3.3-6.3.4). deuflhard2005 studies global convergence under alternative conditions. For just and under-determined systems, several authors considered a strong monotonicity condition: \[ ( g(\theta_1) - g(\theta_2) )^\prime (\theta_1-\theta_2) \geq \mu \|\theta_1-\theta_2\|^2,\] with $\mu > 0$, e.g. solodov2000, polyak2020. Note that when $g = \partial_\theta F$, then $g$ is strongly monotone if, and only if, $F$ is strongly convex. Hence, global convergence under strong monotonicity is related to global convergence under strong convexity of $F$. In that case, $g$ is said to be cyclically monotone Rockafellar. These results do not consider $g(\theta^\dagger) \neq 0$ which is particularly relevant here.
A companion paper, Forneron2022b, considers correctly specified GMM estimation with non-smooth sample moments that may not satisfy Assumption (ref). There are two important differences in that setting: 1) the Jacobian $G_n$ is not defined, and 2) $Q_n$ can have local optima. The methods considered here are not sufficient to find a global optimum, and there is a curse of dimensionality for global convergence. The two papers are complementary. \paragraph{Relation between the different conditions.} Narrowing to the GMM setting specifically, the following shows that the PL inequality holds in the population for correctly specified models under Assumption (ref). A related result is derived under misspecification.
As discussed above, Assumption (ref) (a) implies Assumption (ref) (b). The latter confers most of the properties required for minimizing $Q$. It can be useful to re-write the condition in terms of $g$: Assumption (ref) (b) $\|G(\theta)^\prime W [g(\theta)-g(\theta^\dagger)]\| > \rho \underline{\sigma}\|\theta-\theta^\dagger\|$. For correctly specified models, $g(\theta^\dagger) = 0$ implies $G(\theta)^\prime W g(\theta) = \partial_\theta Q(\theta)$. The only critical point is $\theta = \theta^\dagger$. Hence, Assumption (ref) excludes local optima and saddle points when the model is correctly specified.\footnote{A critical point is a $\theta$ such that $\partial_\theta Q(\theta) = 0$. Assuming $Q$ is twice differentiable, it is a local minimum if $\partial^2_{\theta,\theta^\prime} Q(\theta)$ is positive semidefinite, maximum if $\partial^2_{\theta,\theta^\prime} Q(\theta)$ is negative semidefinite, and a saddle point if $\partial^2_{\theta,\theta^\prime} Q(\theta)$ is indefinite, i.e. has both positive and negative eigenvalues.}
Proposition (ref) shows that Assumption (ref) (b), together with bounds on $W$ and Lipschitz continuity of $G$ imply the PL inequality (1) for $Q$. In addition, (2) implies global identification and is needed to derive the convergence rate of $\theta_k$. Strong convexity also implies (1) and (2).
Proposition (ref) gives a condition under which quasar-convexity and the PL inequality imply Assumption (ref) (b). On compact sets, Assumption (ref) (i), (iii), (iv) together imply a $C_2 >0$ exists for correctly specified models. Assumption (ref) (b) does not imply quasar-convexity.\footnote{Quasar-convexity implies $(\theta-\theta^\dagger)^\prime G(\theta)^\prime W \overline{G}(\theta)(\theta-\theta^\dagger) \geq \frac{\mu}{2\lambda} \|\theta-\theta^\dagger\|^2$ for correctly specified models. This is more restrictive than Assumption (ref) (b) when $\theta$ is not scalar.} The following considers strong monotonicity and introduces a strong injectivity condition:
It can be shown that the strong injectivity property holds on compact convex sets under the Gale-Nikaidô-Fisher-Rothenberg global identification conditions: $\text{det}(G(\theta)) > 0$ and $G(\theta)$ positive quasi-definite, for all $\theta \in \mathbb{R}^{d_\theta}$, where $\text{det}$ is the determinant.\footnote{$G$ is positive quasi-definite if, and only if, $G+G^\prime$ is positive definite. See fisher1966, rothenberg1971; and komunjer2012 for a discussion and alternative conditions.}
For over-identified models, ((ref)) is not sufficient. As discussed above, the weighting matrix $W$ plays a role in the convergence properties. The following extends ((ref)) appropriately:
Relative to ((ref)), the additional term ensures that $g(\cdot)$ is one-to-one in the row space of $G(\cdot)^\prime W$. Taking $(\theta_1,\theta_2) = (\theta,\theta^\dagger)$ yields Assumption (ref) (b). Note that ((ref)) implies that Assumption (ref) holds for gn (Lemma (ref)). ((ref)) is more challenging to verify than ((ref)) as it involves the weighting matrix $W$ and the Jacobian $G$. If ((ref)) holds for a just-identified subset of moments, then it is possible to regularize $W$ so that ((ref)) holds (Lemma (ref)).
Figure (ref) summarizes the results of Propositions (ref), (ref), (ref). Since $Q_n(\hat{\theta}_n) = 0$ for just-identified models that are correctly specified, the relationship also applies in the finite samples problems where these conditions are met. When $g$ and $\theta$ are scalar, Assumption (ref) implies strict monotonicity, $g$ is either increasing or decreasing, but does not imply convexity of $Q$, however, as the MA example below will illustrate.
It remains to determine if Assumption (ref) (b) is minimal for global convergence, or if can be weakened further. The following condition is necessary for gd and other gradient-based optimizers of the form ((ref)) to be globally convergent:
The following shows, under regularity conditions, that (N) implies Assumption (ref) (b).
The case of misspecified models is more complicated, as the following shows that the equivalence between the PL inequality and Assumption (ref) (b) is not automatic.
Proposition (ref) (1) is only informative when the amount of misspecification is moderate, i.e. $\varphi < \rho^2 \underline{\sigma}^2/[\overline{\lambda}_W L^2]$. When this holds, there are no local optima besides $\theta^\dagger$. It also implies the PL inequality (1') holds. To recover convergence for $\theta$, the lower bound in (2) should be informative which further requires $\sqrt{\varphi} < C_2 / C_4$.\footnote{The derivations give the following bounds $C_2 = 1/2 \frac{\rho^2\underline{\sigma}^2}{\overline{\sigma}^2 \overline{\lambda}_W}$ and $C_4 = \overline{\lambda}_W^{1/2} L$ so that the condition reads $\sqrt{\varphi} < 1/2 \rho^2\underline{\sigma}^2 [\overline{\sigma}^2\overline{\lambda}_W^{3/2} L]^{-1}$. It is possible to relax this condition at the cost of more complicated derivations using a combination of global and local convergence arguments.} The degree of non-linearity - measured by $L$ - and the choice of weighting matrix - measured by $\overline{\lambda}_W,\underline{\lambda}_W$ and $\varphi$ - constrain the amount of misspecification permitted to get informative bounds. For correctly specified models, $Q_n(\hat{\theta}_n) = o_p(1)$ implies that both (1') and (2) hold asymptotically.
\paragraph{Further characterization of Assumption (ref) (Just-Identified).}
Like star-convexity, Assumption (ref) is stated relative to the unknown $\theta^\dagger$. The following Proposition gives several conditions under which Assumption (ref) (a) holds and properties implied by these conditions.
Condition (a) does not require knowledge of $\theta^\dagger$ and implies that $g(\cdot)$ is one-to-one. The latter is often assumed for indirect inference.\footnote{See e.g. gourieroux1993, Assumption (A4).} Condition (a) also implies ((ref)) with $\mu = \underline{\sigma}$. When the Jacobian can be linearly rearranged into a symmetric positive definite matrix $S(\theta) = U^{-1} G(\theta) V^{-1}$, then condition (a) holds. These problems can be thought of as implicitly convex in the special case where where $S$ is the second derivative of a convex function. For a given $\theta \in \mathbb{R}^{d_\theta}$, the decomposition (b) always exists: the singular value decomposition gives $G(\theta) = U(\theta) S(\theta) V(\theta)$ where $U(\theta),V(\theta)$ are unitary and $S(\theta)$ is diagonal with positive entries. A lesser known result, due to frobenius1910 shows that any square matrix can be written as the product of two real symmetric matrices; here $G(\theta) = S_1(\theta) S_2(\theta)$. The Jordan normal form of $G(\theta)$ can be used to compute this factorization bosch1986. If $G(\theta)$ is invertible, for all $\theta \in \mathbb{R}^{d_\theta}$, and $U,V$ or one of $S_1,S_2$ do not vary with $\theta$, in the singular value or Frobenius decomposition, then (b) holds. Under condition (c), $g$ is cyclically , and thus strongly, monotone.
Strong convexity is preserved by affine transformations and reparameterization that satisfy particular component-wise monotonicity constraints on the reparameterization boyd2004. Proposition (ref) shows that Assumption (ref) (a) is also preserved by affine transformations and moderately non-linear reparameterizations $h$. Hence, under Assumption (ref) (a), optimization should be locally robust to the choice of parameterization. Assumption (ref) (b) is preserved if $h$ is strongly injective, a mild requirement. In particular, invertible affine transformations preserve Assumption (ref) (b). Similar statements for overidentified models can be found in Propositions (ref), (ref). Propositions (ref) and (ref) together imply that if $Q$ is strongly convex for a particular parameterization, e.g. reduced-form coefficients, then Assumption (ref) (b) holds for $Q \circ h$ where $h$ satisfies the conditions above, where $h$ is the mapping from reduced form to structural coefficients.
The global convergence results hinge on Assumption (ref) (b) as it confers the objective several key properties for optimization. In some cases it may be feasible to verify analytically that one of the conditions in Figure (ref) or Proposition (ref) hold. For some models, it is possible to construct moments that identify the parameters, typically using injectivity arguments. In that case, ((ref)) holds which implies Assumption (ref) (b) holds, under regularity conditions.
For more complex models, it may only be possible to evaluate numerically over a representative set of points, whether one of these conditions is likely to holds, or not. Since Assumption (ref), and its sample counterpart Assumption (ref), depend on the unknown minimizer or $Q$, resp. $Q_n$, it is not possible to check the conditions numerically before the performing the estimation. It is possible to check a stronger condition which does not take an estimate $\hat{\theta}_n$ as input, however.
In the main results, the constant $\rho \underline{\sigma}$ can be arbitrarily small. In practice, however, when $\rho\underline{\sigma} \to 0$ the convergence rate $(1-\overline{\gamma}) \to 1$ is arbitrarily slow. The following approximates an upper bound for $(1-\overline{\gamma})$, assuming correct specification, and the corresponding number of iterations $\underline{k}$ required to achieve $Q_n(\theta_{\underline{k}}) - Q_n(\hat{\theta}_n) \leq \varepsilon [Q_n(\theta_0) - Q_n(\hat{\theta}_n) ]$ for a user-chosen $\varepsilon \in (0,1)$. In practice, these bounds can be very conservative. The value $\underline{k}$ mainly indicates whether global convergence is practically feasible (e.g. $\underline{k} \leq 10^3$) or not (e.g. $\underline{k} \geq 10^{12}$).
When $\hat{\theta}_n$ is unknown, before the estimation is performed, it is only possible to verify a stronger condition. The following considers a sample analog of ((ref)), introduced above:
for some $\mu_n > 0$. In the following, the finite grid of pairs $\Theta_K = \{ (\theta^1_1,\theta^2_1),\dots,(\theta^1_K,\theta^2_K) \}$ will be used for that purpose. It construction is discussed in more detail below. Suppose $\theta^1_k \neq \theta^2_{k}$ for each $k$, compute:
where $P_{k,n}$ is computed using $\theta^1_k$ and the algorithm of choice.\footnote{Note that $L_{Q,P,k}$ involves a Hessian-vector product, which can be computed using only gradients: $H_n(\theta^1_k)(\theta^1_k-\theta^2_{k}) \simeq [\partial_\theta Q_n(\theta^1_k + \epsilon[\theta^2_{k}-\theta^1_k])-\partial_\theta Q_n(\theta^1_k)]/\epsilon$, for $\epsilon$ small.} Then compute: \[ (1-\overline{\gamma})^2 = \max\left(0,1 - [\hat{\mu}_n \hat{C}_{3,n}]^2 /[4 \hat{L}_{Q,P,n}]\right), \quad \underline{k} \geq \frac{\log(\varepsilon)}{\log(1-\overline{\gamma})}, \] where $\hat{\mu}_n = \min_k \mu_k$, $\hat{C}_{3,n} = \min_k C_{3,k}$, and $\hat{L}_{Q,P,n} = \max_k L_{Q,P,k}$. The normalization using $P_{k,n}$ ensures that these values are invariant to linear reparameterizations of the parameters and/or moments for gn or nr. As a reference, with the normalization linear models have $\hat{\mu}_n=1$ under the standard rank condition. To ensure the product $P_{k,n} G_n(\theta_k^1)$ is well behaved, it is recommended to compute a pseudo-inverse of $G_n(\theta_k^1)^\prime W_n G_n(\theta_k^1)$ in the case of gn. This yields $\mu_k = 0$ and $\overline{\gamma}=0$ when $G_n(\theta_k^1)$ is numerically close to singular in the relevant direction. If the conditions fail or the bounds indicate that convergence is not practically feasible, typically when $\hat{\mu}_n < 10^{-2}$ for gn,\footnote{Dividing $\mu_n$ by 10 approximately reduces $\overline{\gamma}$ by a factor of 100, by a local expansion argument. Convergence becomes significantly slower when $\hat{\mu}_n$ approaches $0$.} gradient-based methods need to be modified to ensure global convergence, using multiple starting values or a hybrid approach with theoretical guarantees, see Forneron2022b for an explicit algorithm in that setting.
\paragraph{Constructing $\Theta_K$.} Take a set $\Theta$ large enough that $\hat{\theta}_n \in \Theta$ is likely. The grid $\Theta_K$ should be dense in $\Theta$ so that, as the number of points $K$ increases, any $\theta \in \Theta$ is arbitrarily close to some value in the grid. For $\Theta = [0,1]^{d_\theta}$, the Sobol and Halton sequences have this property, and are readily available in statistical software (R, Matlab, Python, Julia). In general, when $\Theta = [\underline{\theta}_1,\overline{\theta}_1] \times \dots \times [\underline{\theta}_{d_\theta},\overline{\theta}_{d_\theta}]$, where $\underline{\theta}_1,\overline{\theta}_1$ denote lower (resp. upper) bounds on each coefficient, a sequence can be constructed from the Sobol or Halton sequence, denoted $(\vartheta_{i,k})$, $i=1,\dots,d_\theta$, $k=1,\dots,K$, by setting $\theta_{i,k} = \underline{\theta}_i + (\overline{\theta}_i-\underline{\theta}_i)\vartheta_{i,k}$.
The results are stated for a fixed learning rate. In practice, adaptive choices of $\gamma_k$ are common, using a line search for instance. If the adaptive algorithm is tuned to satisfy the requirements for global convergence, then it is also globally convergent. To preserve convergence properties, additional tuning parameters are typically involved nocedal-wright:06. A backtracking line search, a simple and popular way to set the learning rate nocedal-wright:06, is used as a benchmark comparison for the fixed learning rate used in the applications.
By construction, $J_k^\prime p_k \geq 0$ so that the final $\gamma_k$ decreases the value of the objective function. The while loop terminates once the so-called Armijo condition is met:\footnote{See nocedal-wright:06, Nesterov2018 for discussions.} $Q_n(\theta_{k+1}) \leq Q_n(\theta_{k}) - c \gamma_k J_k^\prime p_k$. For just-identified models, the termination criterion is feasible if $c$ is sufficiently small.\footnote{A sample analog of Proposition (ref) implies that $Q_n(\theta_{k+1}) \leq (1 - \overline{\gamma})^2 Q_n(\theta_k)$, for any $\theta_k \in \Theta$, when $\gamma \in (0,1)$ small enough for some $\overline{\gamma} \in (0,\gamma)$. Proposition (ref) (1)-(2) further imply for just-identified models that $Q_n(\theta_k) - Q_n(\hat{\theta}_n) \leq c_n J_k^\prime p_k$ for some $c_n>0$. The Armijo condition is feasible if $c$ is small enough.} Having $\theta_k = \hat{\theta}_n$ implies $p_k = 0$; the condition holds for any $\gamma_k \in (0,1]$. A common choice is $c = 10^{-4}$, $\gamma_{\text{init}} = 1$, $\rho = 0.8$. These were used in all examples.\footnote{When there are bounds for parameters values, one can set $Q_n(\theta_{k+1}) = +\infty$ if $\theta_{k+1}$ is outside the bounds. Another approach is to project $\theta_{k+1}$ inside the bounds when $\gamma_k$ is too large.}
The first example illustrates the main results using a simple MA(1) process: \[ y_t = e_t - \theta^\dagger e_{t-1}, \quad e_{t} \overset{iid}{\sim} \mathcal{N}(0,1), \quad \theta^\dagger \in (-1,1), \] for $t=1,\dots,n$. $\theta^\dagger$ is the parameter of interest. Set $p \geq 1$, following gourieroux1996, $\theta^\dagger$ is estimated by matching coefficients from an auxiliary AR(p) model: $y_t = \beta_1 y_{t-1} + \dots + \beta_p y_{t-p} + u_t$. For $p = 1$, $\hat\beta_1 \overset{p}{\to} -\theta^\dagger/(1+\theta^{^\dagger 2})$ defines the moment condition: \[ \overline{g}_n(\theta) = \hat\beta_1 + \frac{\theta}{1+\theta^2},\] with Jacobian $G_n(\theta) =(1-\theta^2)/(1+\theta^2)^2 > 0$ for any $\theta \in (-1,1)$ and $G_n(\theta) = 0$ for $\theta \in \{-1,1\}$. It has full rank on any interval of the form $[-1+\varepsilon,1-\varepsilon]$, $\varepsilon \in (0,1)$. However, Figure (ref) shows that the Hessian $\partial^2_{\theta,\theta} Q_n(\theta)$ can be positive, negative, or equal to zero depending on the value of $\theta$ -- $Q_n$ is non-convex, especially when $\overline{g}_n(\theta)$ is large. Now notice that: $\overline{g}_n(\theta) = \partial_\theta F_n(\theta) \text{ where } F_n(\theta) = \hat\beta_1 \theta + \frac{1}{2}\log(1+\theta^2),$ which not a GMM objective but is nevertheless convex on $[-1,1]$, strongly convex on any $[-1+\varepsilon,1-\varepsilon]$, $\varepsilon \in (0,1)$. Hence, $\overline{g}_n$ is cyclically monotone and statisfies Assumption (ref) (a). Note that implicitly, gn minimizes the convex $F_n$ -- whereas nr explicitly minimizes the non-convex $Q_n$. This is specific to the just-identified case ($p=1$), since an $F_n$ cannot be defined in the over-identified case ($p>1$).
Table (ref) shows the search paths for nr and gn with a fixed $\gamma=0.1$ as well as R's built-in optim's bfgs implementation and the bound-constrained l-bfgs-b. nr diverges, because the objective is locally concave at $\theta_0 = -0.6$. This is surprising given how close $\theta_0$ is to $\theta^\dagger$. Although $Q_n$ is locally convex around $\hat\theta_n$, which is useful for local optimization, the corresponding neighborhood can be fairly small from a practical standpoint.
gn converges steadily from the same $\theta_0$. bfgs is more erratic, especially when $\theta_k \simeq -0.5$, i.e. $k = 1$, leading to a search outside the unit circle ($k=2$), before reaching an area where the iterations are better behaved ($k=3$ onwards). While here this is not too problematic, the objective function is well defined outside the bounds, this is more concerning in applications where the model cannot be solved outside the bounds -- this is illustrated in Section (ref). A natural solution is to introduce bounds using l-bfgs-b. The search, however, remains somewhat erratic as seen in the Table. Compare these to bfgs$^\star$ and l-bfgs-b$^\star$ which minimize $F_n$, instead of $Q_n$, using the same optim. Like \textsc{gn}, they steadily converge to $\hat\theta_n$.
For $p=12$, the model is over-identified and the conditions are more challenging to check analytically. Figure (ref) indicates that the strong injectivity condition ((ref)) appear to hold, the value is bounded away from zero, except at the boundary. The choice of weighing matrix affects the constant $\mu$ in ((ref)), as it appears to be smaller with optimal weighting.
Table (ref) shows that nr, bfgs and l-bfgs-b all fail to converge from $\theta_0 = 0.95$, a starting value with negative curvature, with identity weighting.\footnote{l-bfgs-b relies on projection descent which maps search directions outside the unit circle to $-1$ or $1$ where $\partial_\theta Q_n(-1) = \partial_\theta Q_n(1) = 0$, a stationary point for ((ref)).} Compare with gn, which steadily converges to $\hat\theta_n$. Starting closer to the solution, bfgs and \textsc{l-bfgs-b} also fail to converge using $\theta_0 = 0.6$; \textsc{gn} remains accurate (not reported). R codes can be found in Appendix (ref).
\paragraph{The impact of misspecification on estimation.} As discussed in Section (ref), when the degree of misspecification $\varphi$ becomes large, local optima may appear and making gradient-based optimizers non-globally convergent. To illustrate this issue, consider the true DGP: \[ y_t = e_t - \theta_1 e_{t-1} - \theta_2 e_{t-2}, \quad e_{t} \overset{iid}{\sim} \mathcal{N}(0,1), \quad \theta^\dagger \in (-1,1),\] where $\theta_2 \in \{0,0.4,0.8\}$ determines the degree of misspecification. The following sets $\theta_1 = - 0.1$ to ensure invertibility as $\theta_2$ varies. The moments are the same as above with $p=12$.
Figure (ref) compares $Q_n$ and its gradient $\partial_\theta Q_n$, with identity and optimal weighting, at different degrees of misspecification. On intervals $[-1+\varepsilon,1-\varepsilon]$, $Q_n$ has no local optima for $\theta_2 \in \{0,0.4\}$, in line with Proposition (ref). For the larger $\theta_2=0.8$, there are local optima: equal weighting has two (one maximum and minimum), optimal weighting has four (two maxima and minima). Also, the objective becomes flatter as $\theta_2$ increases, making optimization more challenging. Table (ref) shows how this translates into estimation properties. As predicted, gn is robust to moderate misspecification but only converges to a local minimum under heavier misspecification. Other methods (l-bfgs-b, bfgs) systematically fail to converge.
The following revisits the results for random coefficient demand estimation in knittel2014 with the `fake' cereal data generated by nevo2001.\footnote{It available in the R package BLPestimatoR brunner2017. The data consists of 2,256 observations for 24 products (brands) in 47 cities over two quarters in 94 markets. The specification is identical to Nevo's, with cereal brand dummies, price, sugar content (sugar), a mushy dummy indicating whether the cereal gets soggy in milk (mushy), and 20 IV variables. } This is a non-linear instrumental variable regression with sample moment conditions: $\overline{g}_n(\theta,\beta) = \frac{1}{n} \sum_{j,t} z_{jt}[ \delta_{jt}(\theta) - x_{tj}^\prime \beta ]$, where $z_{jt}$ are the instruments, $x_{jt}$ the linear regressors in market $j$ at time period $t$. The $8$ parameters of interest are the random coefficients $\theta$,\footnote{8 parameters are the unobserved standard deviation and the income coefficient on the constant term, price, sugar, and mushy.} which enter $\delta_{jt}$, recovered from market shares $s_{jt}$ using the fixed point algorithm of Berry1995. The $25$ linear coefficients $\beta$ are nuisance parameters concentrated out by two-stage least squares for each $\theta$. The replication sets the maximum number of iterations for the contraction mapping to $20000$ and the tolerance level for convergence to $10^{-12}$. This is important for the optimization to be well-behaved; see e.g. brunner2017, conlon2020. The range of starting values used here is much wider than in these papers,\footnote{conlon2020 draw “starting values from a uniform distribution with support 50% above and below the true parameter value.”} which explains why optimizers are more prone to crashing here than in their replications. Initial values are constructed as follows: the Sobol sequence generates values in $[0,1]^8$, the coefficients for standard deviations are adjusted to lie in $[0,10]$, those for income in $[-10,10]$. Values for which the contraction mapping produces an error are discarded until $50$ valid starting values are available.
Table (ref) and Figure (ref) compare the performance of quasi-Newton (bfgs), Nelder-Mead (nm), Simulated-Annealing (sa), and Nelder-Mead after Simulated-Annealing (sa+nm), using R's default optimizer optim, with Gauss-Newton (\textsc{gn}) and Gradient-Descent (\textsc{gd}) for 50 different starting values.\footnote{The solution of the contraction mapping is not well defined for all values in $\Theta$, so we use the first $50$ values produced by the Sobol sequence such that $\delta_{jt}$ is finite for all $j,t$. } As reported in knittel2014, optimization can crash often.\footnote{The optimizers will crash when the fixed point algorithms fail to return finite values. This is typically the case when the search direction was poorly chosen at the previous iteration.} Crashes could be avoided using error handling (try-catch statements). However, this may not be enough to produce accurate estimates as the next application will illustrate.\footnote{conlon2020 illustrate that modifications to the fixed-point algorithm and specific optimizer implementations to handle near-singularity of the Hessian can also improve performance for \textsc{bfgs}.} Only \textsc{gn} systematically produces accurate estimates; \textsc{bfgs} crashes 46% of the time and has one highly inaccurate estimate. Derivative-free optimizers (\textsc{nm}, \textsc{sa}, \textsc{sa+nm}) can produce inaccurate estimates. \textsc{gd} can be very slow to converge. Using a backtracking line search, \textsc{gn} converges in 11 iterations on average, compared to 8816 for \textsc{gd} -- which has a higher maximum number of iterations set at 10000, compared to 150 for \textsc{gn}. Increasing the maximum number of iterations for \textsc{gd} would improve the estimates at the expense of further computation time.
Figure (ref), illustrates the convergence of gn for the first 5 starting values. In line with the predictions of Theorem (ref), though $Q_n$ is non-convex, gn iterations steadily converge to the solution. This type of “Gauss-Newton regression” is related to salanie2022 who compute two-stage least-squares for linearized BLP.
The second application revisits moran2018's estimation of a model with endogenous total factor productivity (TFP) growth moran2018. They estimate parameters related to Research and Development (R&D) by matching the impulse response function (IRF) of an identified R&D shock to R&D and TFP in a small-scale Vector Auto-Regression (VAR) estimated on U.S. data.
The parameters of interest are $\theta = (\eta,\nu,\rho_s,\sigma_s)$ which measure, respectively, the elasticity of technology creation to R&D, R&D spillover to adoption, the persistence coefficient and size of impulse to the R&D wedge. The sample moments are $\overline{g}_n(\theta) = \hat \psi_n - \psi(\theta)$, $\hat \psi_n$ and $\psi(\theta)$ are the sample and predicted IRFs, respectively. The latter is computed using Dynare in Matlab. To minimize $Q_n$, the authors use Sims's csminwel (sims in the Table, Figures)\footnote{Details about csminwel and code can be found at: \url{http://sims.princeton.edu/yftp/optimize/}.} algorithm with a reparameterization which bounds the coefficients.\footnote{The replication uses the mapping $\theta_j = \underline{\theta}_j + \frac{\overline{\theta}_j-\underline{\theta}_j}{1 + \exp(-\vartheta_j)}$, where each $\vartheta_j$ is unconstrained. The original study relied on $\theta_j = 1/2(\overline{\theta}_j + \underline{\theta}_j) + 1/2(\overline{\theta}_j - \underline{\theta}_j) \frac{\vartheta_j}{\sqrt{1 + \vartheta_j^2}}$, which we found to make optimizers very unstable.} Although this type of reparameterization is commonly used, the Jacobian is singular at the boundary; this matters for both local and global convergence, according to the results. As in the demand estimation, initial values are constructed using the Sobol sequence and adjusted to match the bounds used in the original study, reported in the last two rows of Table (ref).
In the original paper, the authors initialize the estimation at $\theta_0 = (\eta_0,\nu_0,\rho_{s0},\sigma_{s0}) = (0.20,0.20,0.30,0.10)$, very close to $\hat\theta_n$. Here, 50 starting values are generated within the bounds in Table (ref). The model is estimated using csminwel and the same set of optimizers used in the previous replication. Table (ref) reports the results with and without the non-linear reparameterization. Similar to the MA(1) model with $p=12$, without the reparameterization, several optimizers return values outside the parameter bounds, which motivates the constraints in these cases. gn correctly estimates the parameters for all starting values but crashes twice for starting values for which both $\eta$ and $\nu$ are close to their lower bounds where the Jacobian is nearly singular. With the reparameterization, gn crashed more often, nines times in total, but is otherwise accurate. With backtracking line search, crashes are fewer for gn, and converges in 14 iterations, on average, with or without reparameterization, compared to 832 for gd without reparameterization and 2912 with reparameterization (both with cap of 10000). The crashes might also occur at values strictly within the parameter bounds for which Dynare cannot solve the model and returns an error. There is no obvious way to modify gn or \textsc{gd} to avoid this problem.
The other two gradient-based optimizers, bfgs and sims(csminwel), never crash because of better error handling in Matlab. They produce valid estimates less often than gn. Figure (ref) illustrates that csminwel is sensitive to reparameterization. Likewise, derivative-free methods can be inaccurate, as illustrated in Table (ref) and Figure (ref); some crashes occur despite Matlab's error handling. Finally, Figure (ref) shows 5 optimization paths for which gn does not crash with and without the reparameterization. Appendix (ref) gives additional results for larger values of $\gamma \in (0,1]$ and error handling.
Table (ref) illustrates the strong injectivity conditions, Assumption (ref) (b), and convexity for the MA(1) model and the two empirical applications. A grid of $100$ Sobol points was used to construct values within parameter bounds, respectively sb and lb, as described in Section (ref); the first grid value is enforced to take only values from the bounds. The objective $Q_n$ is locally convex at $\theta$ if the Hessian is positive definite, i.e. $H_n(\theta)>0$. The Table reports an estimate for $\mu_n$, $\rho\underline{\sigma}$ and the proportion of grid values where $Q_n$ is locally convex. To evaluate $\rho\underline{\sigma}$, the same step from Section (ref) were used setting $\theta_k^2 = \hat{\theta}_n$ for all $k$. The $P_{k,n}$ for gn was used so that the $\mu_n$, $\rho\underline{\sigma}$ reported here are invariant to linear reparameterizations of both the parameters and the moments.
For the MA(1) model, strong injectivity and Assumption (ref) (b) fail at the boundary where $\theta = \pm 1$. This is visible in the results for lb. For sb, both conditions hold as illustrated in the Table. Optimal weighting has some effect on the conditions and the predicted convergence properties. Convexity fails more often with a single moment condition ($p=1$).
For BLP, the conditions appear to hold and predict fast convergence for sb, used in conlon2020, where $Q_n$ is almost everywhere locally convex. With wider bounds (lb), convexity almost always fails, but the estimates for $\mu_n$, $\rho\underline{\sigma}$ are very close to sb. This confirm the good optimization properties for gn reported above.
For the DSGE model, $\mu_n$, $\rho\underline{\sigma}$ are of the same order of magnitude as the other applications without reparameterization. With reparameterization, the conditions can fail at the boundary which is visible in the Table under lb. With and without reparameterization, $Q_n$ is rarely locally convex, which confirms the challenges bfgs and csminwell can have.
In both empirical applications, the estimates for $\overline{\gamma}$ and $\underline{k}$ tend to be very small and large, respectively, despite $\mu_n$, $\rho\underline{\sigma}$ being away from zero. This reflect the large amount of non-linearity, measured by $C_{3,K}$ and $L_{Q,P,K}$. As discussed in Section (ref), these estimates can be fairly conservative which is clearly the case here. Also, because the moments are evaluated numerically, using a fixed-point algorithm for BLP, and the derivatives are computed by finite differences, the second-order derivatives can be fairly inaccurate. This issue is explained in Appendix (ref). Innacurate second-order derivatives can make optimizers like bfgs and csminwell numerically unstable, but will also affect the value of $L_{Q,P}$, which tend to be very large in the empirical applications resulting in a very conservative bound for $\overline{\gamma}$.\footnote{The estimate is $L_{Q,P,K} = 2 \cdot 10^4$ for BLP with large bounds, and $L_{Q,P,K} = 9$ with small bounds.}
Non-convexity of the GMM objective function is considered to be an important challenge for structural estimation. This paper considers alternative conditions under which there are globally convergent algorithms. The results are robust to non-convexity, moderately non-linear one-to-one reparameterizations, and moderate misspecification. Though off-the-shelf methods might fail to converge due to the non-convexity of the optimization problem, the paper has shown that this does not necessarily imply that it will be difficult in practice. Econometric theory emphasizes the role of the weighting matrix $W_n$ on the statistical efficiency of the estimator $\hat\theta_n$. Here, Assumption (ref) may or may not hold, depending on the choice of weighting matrix $W_n$. Its condition number $\kappa_W$ also affects local convergence which highlights an important role for the weighting matrix: it may facilitate or hinder the estimation itself.