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.
97,173 characters · 18 sections · 75 citation commands
Noisy, Non-Smooth, Non-Convex Estimation of Moment Condition Models
JEL Classification: C11, C12, C13, C32, C36.\newline Keywords: Generalized and Simulated Method of Moments, Non-Asymptotic bounds.
\baselineskip=18.0pt \thispagestyle{empty} \setcounter{page}{0}
Generalized Method of Moments (GMM) and Simulated Method of Moment (SMM) estimations are routinely used to evaluate structural economic models on empirical data. Textbook examples include estimating models of dynamic discrete choice hotz1993, market equilibrium Berry1995, general equilibrium christiano2005, and export participation eaton2011. In practice, estimates are computed using numerical optimization methods. The focus of this paper is on a common problem faced by practitioners: optimization for GMM, SMM and related methods tend to be very difficult, and while econometric theory often takes optimization for granted, practitioners do need to successfully minimize the GMM objective function to get valid empirical results.
Three issues are common to structural estimations. Oftentimes i. the objective function is non-convex: fast convex optimizers can return local optima which are not consistent estimates, ii. because of sampling uncertainty and simulation noise: optimization can either succeed or fail with some probability (even if the population problem is well posed), and iii. the sample objective function can be non-smooth or discontinuous: fast derivative-based methods become unstable. Issue i. depends on both the model and choice of moments. Also, convexity is not invariant to non-linear reparameterizations of the coefficients. Issues ii. and iii. are finite sample by nature; iii. is common in SMM estimations of discrete choice models.
This paper tackles these estimation challenges in two steps. First, it proposes a new Gauss-Newton algorithm that combines non-smooth moments with smoothed Jacobian estimates. The algorithm is shown to converge quickly, in an area that is nearly as large as in the smooth population problem. The resulting estimates are asymptotically unbiased under mild conditions on the bandwidth. Existing methods that rely on smoothing produce estimates that can exhibit asymptotic bias unless the bandwidth vanishes quickly enough; faster than the optimal rate for optimization. Such undersmoothing is not required here, making the estimates more robust to the choice of bandwidth. The local convergence analysis relies on standard local identification conditions.
Second, the paper shows that augmenting each Gauss-Newton iteration with a grid-search step results in a globally convergent algorithm. After a finite number of iterations, the combined steps preserve the fast local convergence rate, albeit with a slightly slower rate. This contrasts with commonly used global optimizers that converge slowly. The global convergence results leverage standard global identification conditions.
The idea here is that the local step alone converges quickly to the target estimates, but only from a local initial value. The global step alone is globally convergent from any initial value, but slow. When using them concurrently, from an arbitrary initial value, the global step can improve upon the local one for a certain number of iterations -- determined by the global identification condition and the choice of sequence used in the global step. Afterwards, convergence is driven by the local step. This makes convergence more robust without sacrificing too much performance. Here, the algorithm adaptively transitions from global to local optimization without user input. The fast convergence occurs without smoothing the sample moments in either the local or global step. This is achieved by exploiting the specific formulation of the GMM objective function.
This paper adopts a non-asymptotic (finite sample) approach to study the optimization and statistical properties of the algorithm simultaneously. This allows to better understand the effects that tuning parameters can have on both the estimation and the estimates. For instance, smoothing should facilitate optimization but could also degrade the properties of the estimator. Here, the derivations show how both are affected by the choice of smoothing parameters. The finite-sample approach is also conceptually attractive because it gives results that are valid across repeated samples. Textbook optimization results often consider fixed objective functions nocedal-wright:06,bertsekas2016. However, in empirical work, each dataset is associated with its own, unique, objective function. In that sense, guarantees over repeated samples provide a sense of robustness against sampling uncertainty.
The first extension of the main results considers momentum, also known as the Polyak heavy-ball, which is often used to accelerate the convergence of stochastic gradient descent algorithms. Here, it can be used to reduce the effect of sampling noise on optimization while maintaining the rate of convergence. In practice, it smoothes the optimization path, optimal values for the tuning parameters are tabulated. The second extension proposes a computationally attractive Monte Carlo quasi-Newton implementation of the algorithm which makes it applicable to models where smoothing is not analytically tractable.
Simulated and empirical examples illustrate the results. A dynamic discrete choice model illustrates the algorithm's advantages over a benchmark which minimizes a smoothed objective function. From a distant starting value, the proposed algorithm performs well whilst the benchmark systematically fails to find consistent estimates. When initialized at the true value, the benchmark is consistent but sensitive to the choice of smoothing parameter. In contrast, the algorithm is fairly robust to the smoothing parameter. This is in line with theory. A thorough comparison with several global optimizers highlights the good performance of the proposed algorithm - both in terms of numerical accuracy and computation time. Similar to this paper, several of these optimizers are hybrid: they combine local and global iterations. However, they often rely on heuristic rules and lack theoretical guarantees for convergence, unlike the proposed method. This is reflected in the relative performance of these methods. Then, a simple Aiyagari model illustrates the results on a small but fairly challenging estimation problem where commonly used optimizers fail to converge. The algorithm performs well for a range of tuning parameters. Finally, the empirical application considers a model of joint retirement decision from honore2018. Their estimation of $30$ coefficients takes more than 5 hours whereas the algorithm introduced in this paper finds accurate estimates within 11 minutes.
Lastly, note that the results cover smooth GMM estimations as a special case. The idea to combine a local and global search step to design a fast converging optimizer naturally extends to non-GMM settings such as smooth maximum likelihood and non-linear least-squares estimations. A key difference is that the local optimization results here rely on approximate local linearity of moments whereas in these settings derivations would require local strong convexity and a convex optimizer. Convexity requires further assumptions beyond the first derivatives and can be more restrictive than the assumptions used here. This is illustrated for a simple quantile estimation below. Because these estimations would involve a different set of assumptions and derivations, they are not covered in this paper.
\paragraph{Outline of the paper.} After a brief literature review, Section (ref) briefly describes the setting and introduces the algorithm with a simple illustration. Section (ref) provides the main assumptions, then derives local and global convergence results in finite samples. Section (ref) provides two extensions of the algorithm and derives their properties. Section (ref) provides simulated and empirical applications to illustrate the properties and performance of the method compared to a number of popular global optimizers. The Appendices provide proofs as well as additional simulation and empirical results.
Optimization is a defining feature of non-linear models evaluated by M-estimation:
where $Q_n$ is the sample objective function, $\theta$ are parameters of interest belonging to a compact, finite-dimensional set $\Theta$. While some estimators have closed-form solution, this is the case of OLS regressions, non-linear models generally require numerical methods to find $\hat\theta_n$. Statistical software provide empirical researchers with a wide array of optimizers.
When choosing among these methods, two properties are important to consider: convergence guarantees, typically over a class of functions, and convergence rates within that class. The former matters since a failed optimization yields inconsistent estimates, the latter is critical when $Q_n$ is costly to evaluate. The following expands on these two properties.\footnote{Since most optimization results consider a fixed objective $Q_n$, sampling uncertainty is not discussed yet.} The goal here is to find, after $b \geq 1$ iterations, a $\theta_b$ such that:
at a desired tolerance level $\text{err} >0$. Suppose $Q_n$ has a certain set of features, one could ask: how large does $b$ need to be to guarantee the error $\|\theta_b-\hat\theta_n\|$ is at most $\text{err}$? This is the complexity of the optimization problem Nemirovsky1983. If $Q_n$ is smooth and strongly convex,\footnote{$Q_n$ is strongly convex if its Hessian $H_n$ is continuous and $0 < \underline{\lambda} \leq \lambda_{\min}[H_n(\theta)] \leq \lambda_{\max}[H_n(\theta)] \leq \overline{\lambda} < \infty$.} then $b \geq O( |\log(\text{err})| )$ so that to divide $\text{err}$ by $10$, an additional $O(\log[10])$ iterations are required. Gradient descent, and (quasi)-Newton methods can achieve this rate of convergence. However, when $Q_n$ is non-convex, they may only converge to a local optimum. If $Q_n$ is $r$-times continuously differentiable but otherwise arbitrary, then $b \geq O( \text{err}^{-d_\theta/r} )$ iterations are required, where $d_\theta = \text{dim}(\theta)$. For continuous problems, $r=1$, this rate is achieved by a grid search. Here, the curse of dimensionality is apparent: to divide $\text{err}$ by $10$, the number of iterations needs to be multiplied by $10^{d_\theta} = 10,000$ for $d_\theta=4$.\footnote{To put this into context, the average GMM estimation published between 2016 and 2018 in the American Economic Review involved $d_\theta = 10$, the median $d_\theta = 6$ jjlz22.} Nesterov2018 and andrews1997 illustrate that as $d_\theta$ increases, optimization quickly is practically infeasible under this rate of convergence. Computer software provides a number of global optimizers, however, as griewank1981 points out, those that do not cover the set $\Theta$ are “heuristic,” in the sense that global convergence is not guaranteed.\footnote{A sequence $(\theta_k)_{k \geq 0}$ will be said to cover $\Theta$ if it is dense, i.e. $\sup_{\theta \in \Theta} (\inf_{ 0 \leq \ell \leq k}\|\theta-\theta_\ell\|) \to 0$, as $k \to \infty$. These types of methods are referred to as space-covering in griewank1981.} Several popular methods are used as benchmark for the method introduced in this paper in Section (ref), and their known theoretical properties are reviewed there.
Convergence guarantees are a very desirable property since they ensure that a given method will perform well in a variety of settings. A common denominator of GMM estimations are the conditions required for consistency and asymptotic normality. This paper solely relies on these conditions to derive convergence results so that the algorithm proposed below performs well in a range of empirical settings and datasets.
Since many GMM problems are non-convex, a number of authors have proposed strategies to speed-up convergence. Building on robinson1988, andrews1997 proposed a stopping-rule to produce a consistent first-step estimate, with $\text{err} = o_p(1)$, which is followed by a few Newton-Raphson iterations. The main challenge is in computing the consistent estimate which still requires a global search. chernozhukov2003 introduce a quasi-Bayesian framework where MCMC sampling replaces optimization. The random-walk Metropolis-Hastings algorithm converges under weak conditions. However, convergence rates, which measure performance, are mostly derived for log-concave distributions.\footnote{See e.g. mengersen1996, brooks1998, belloni2009.}
A common strategy to handle non-smooth sample moments is to use smoothing. The complexity bounds above then apply to the smoothed optimization problem. mcfadden1989 and bruins2018 suggest to replace an indicator function by a smooth approximation in simulation-based estimation of discrete choice models. Several papers consider smoothing in quantile regressions kaplan2017,fernandes2021,he2021. A poor choice of smoothing parameter can have a significant impact on the estimates and their statistical properties kaplan2017. Another approach is to apply subgradient methods to the unsmoothed quantile loss function lee2022. Other applications of smoothing for econometric estimation include fermanian2004, altissimo2009, kristensen2012. For these, it is possible to apply numerical bias correction to reduce smoothing bias as in kristensen2017; this is asymptotically valid under further regularity conditions bruins2018. For indirect inference, frazier2019 consider an alternative simulation method for discrete outcomes to get smooth simulated moments. This applies to a specific class of models and moments. ackerberg2009, sauer2021 propose to use importance sampling to smooth the objective function; however, a poor choice of sampling distribution inflates simulation noise ackerberg2009. hong2015 study the asymptotic properties of M-estimates defined by a first-order condition computed by finite-differences; optimization properties are not derived. Finite-differences yield non-smooth derivative estimates when the objective is itself non-smooth, this makes the properties of the associated derivative-based optimizer unclear. Here, the local and global search step involve unsmoothed moments making the estimates more robust to the choice of bandwidth. This is illustrated in the next Section. Despite using unsmoothed moments, convergence results are comparable to the smooth population problem, in finite samples. This contrasts with generic non-smooth optimizers, including subgradient methods, which are generally slower than smooth ones.
Let $\overline{g}_n(\theta) = \frac{1}{n} \sum_{i=1}^n g(\theta;x_i)$ be the sample moments, where $x_i$ are iid random variables and $\theta$ are the parameters of interest. The goal is to find an estimate of $\theta^\dagger \in \Theta$ solving the moment condition $g(\theta^\dagger):=\mathbb{E}[g(\theta^\dagger;x_i)]=0$. In practice, this entails finding an approximate minimizer $\hat\theta_n \in \Theta$ of the sample objective function $Q_n(\theta) = \|\overline{g}_n(\theta)\|_{W_n}^2 = \overline{g}_n(\theta)^\prime W_n \overline{g}_n(\theta)$ satisfying:
where $(W_n)_{n \geq 1}$ is a sequence of symmetric and strictly positive definite weighting matrices with limit $W >0$. To simplify the analysis, $W_n$ is assumed not to depend on $\theta$.\footnote{This excludes continuously-updated and allows for two-step estimations.} Algorithm (ref) is designed to find such an approximate minimizer. It relies on the Jacobian $G_{n,\varepsilon}$ of smoothed moments $\overline{g}_{n,\varepsilon}$ using a bandwidth $\varepsilon$, computed by convolution smoothing. Let $\phi$ be the standard Gaussian density, the smoothed moments and Jacobian are:
The computation of $G_{n,\varepsilon}$ is discussed in the next subsection and Section (ref), Algorithm (ref). Algorithm (ref) alternates between a local and a global search step. The main idea is that, from the current guess $\theta_{b}$, the local step proposes a new candidate $\theta_{b+1}$. When the moments are non-smooth, the smoothed Jacobian $G_{n,\varepsilon}$ is used to approximate the unfeasible Gauss-Newton local iteration. The global step proposes a competing candidate $\theta^{b+1}$, if the fit is strictly better, $\theta_{b+1}$ is discarded and $\theta^{b+1}$ becomes the current iterate. Otherwise, the value from the local step is kept. Note that, away from the solution, the local step might produce values $\theta_{b+1}$ outside of $\Theta$. One possibility is to set $\|\overline{g}_n(\theta_{b+1})\|_{W_n} = +\infty$, then the global step always binds. Another standard approach is to project $\theta_{b+1}$ inside $\Theta$, that is to set each coefficient that is outside its bounds slightly below/above its respective upper/lower bound. Appendix (ref) provides R code for Algorithm (ref) with the extensions from Section (ref).
Algorithm (ref) returns the best iterate $\tilde{\theta}_n$, which is a valid estimator if at least one of the $\theta_j$ is close to minimizing $Q_n$. For inference, standard errors for the $j$-th coefficient of $\theta$ can be readily computed as $\text{se}(\tilde{\theta}_{n,j}) = \sqrt{\tilde{V}_{n,jj}/n}$; $\tilde{V}_{n,jj}$ denotes the $j$-th diagonal entry of $\tilde{V}_{n}$. The variance $\Sigma$ can be estimated in a standard fashion, e.g. $\Sigma_n = 1/n \sum_{i=1}^n g(\tilde{\theta}_n;x_i) g(\tilde{\theta}_n;x_i)^\prime - \overline{g}_n(\tilde{\theta}_n)\overline{g}_n(\tilde{\theta}_n)^\prime$ for an iid sample. Note that Gauss-Newton iterations alone can be globally convergent without convexity under rank conditions jjlz22. Here, the global step ensures that convergence will hold over a broader class of problems. A common approach to finding a global solution is the multi-start strategy which runs a local optimizer many times from different starting values. The number of starting values required for global convergence depends on 1) the choice of algorithm and 2) the objective function; so it is unknown in practice. This paper simply adds a grid search step to get global convergence in a single run of the Algorithm. As discussed earlier, covering the set $\Theta$ is required to guarantee convergence. Running both steps at each iteration yields faster convergence than using them in succession as in the two-step approach of robinson1988, andrews1997. Differences with these methods are discussed in more detail in Appendix (ref).
Potential applications of Algorithm (ref) include quantile regressions and SMM estimations of static or dynamic discrete choice models. Some smooth problems are made non-smooth by numerical approximations. For instance, if an income process is discretized to solve for the optimal policy, the approximated model is discrete and the moments are potentially non-smooth. Continuous models with discrete changes can also exhibit non-smoothness. In finance, this is the case of jumps-diffusion models estimated by SMM carrasco2000.
The algorithm has three tuning parameters: the learning rate $\gamma$, the bandwidth $\varepsilon$, and the covering sequence $(\theta^b)_{b \geq 0}$. The following provides specific recommendations for each of them.
\paragraph{Learning rate.} The Gauss-Newton normalization $[ G_{n,\varepsilon}(\theta_b)^\prime W_n G_{n,\varepsilon}(\theta_b) ]^{-1}$ ensures local convergence for any choice of $\gamma \in (0,1)$. In contrast, gradient descent methods require $\gamma > 0$ to be sufficiently small. The main tradeoff in choosing $\gamma$ is that a larger learning rate leads to faster local convergence but over a smaller local area around the solution. In the applications, a smaller value $\gamma=0.1$ balances speed and the area of convergence quite well.
\paragraph{Bandwidth.} For optimization, setting $\varepsilon = c \cdot n^{-1/4}$ is rate optimal. The constant $c$ can be chosen to minimize a mean squared prediction error as follows. Construct a grid $\theta^1,\dots,\theta^B$, this can be the same as the covering sequence described below. Then for each $j = 1,\dots,B$ find $\theta^j_1,\dots,\theta^j_K$ the K-nearest neighbors to $\theta^j$ from the same grid and compute: \[ \text{MSE}_{j}^\ell(\varepsilon) = \| \overline{g}_n(\theta^j) - \overline{g}_n(\theta^j_\ell) - G_{n,\varepsilon}(\theta^j_\ell)(\theta^j - \theta^j_\ell) \|^2. \] Compute the average $\text{MSE}(\varepsilon)$ of $\text{MSE}_{j}^\ell(\varepsilon)$ over $j$ and $\ell$. This measures how well $G_{n,\varepsilon}$ can predict $\overline{g}_n(\theta^j)$ from $K$ nearby directions $\theta^j_\ell$. The moments and Jacobians can be computed once, ideally in parallel, before evaluating the $\text{MSE}$s. Take $\varepsilon_1,\dots,\varepsilon_m = c_1 \cdot n^{-1/4},\dots,c_m \cdot n^{-1/4}$ and pick the $\varepsilon$ among these which minimizes the $\text{MSE}$ above. For regressions, it can be useful to normalize variables to have mean $0$ and variance $1$ so that the bandwidth is scale-free. The applications use $B = 250$, $K=10$, set $c_1,\dots,c_m$ to $8 = 2^{3},2^2,\dots,2^{-4},2^{-5} = 0.003$.
The choice of a MSE criterion is related to the role of $G_{n,\varepsilon}$ in the local step. It is used to approximate $\overline{g}_n(\theta)$ around $\theta = \theta_{b-1}$, the local step minimizes this linear surrogate. The MSE criterion measures the accuracy of the linear approximation.
\paragraph{Covering Sequence.} the choice of sequence $(\theta^b)_{b \geq 0}$ matters much for the performance of Algorithm (ref). The main property required is that the so-called dispersion of the sequence: \[ D(\theta^1,\dots,\theta^b) = \sup_{\theta \in \Theta} \left( \inf_{j=1,\dots,b} \|\theta-\theta^j\| \right) \] be as small as possible for each $b$. The dispersion measures how well the sequence $(\theta^1,\theta^2,\dots)$ covers the parameter space $\Theta$, i.e. how quickly, as $b$ increases, the distance between any $\theta$ and the closest $\theta^j$ decreases. Intuitively, the goal is to have one of $\{\theta^1,\dots,\theta^b\}$ that is as close as possible to $\hat{\theta}_n$, wherever it may be within $\Theta$, at each iteration $b$. This is only possible if the parameter space is bounded, so it is necessary to set bounds for the parameters. Bounds are always required for global non-convex optimization. Many parameters are naturally bounded: the discount rate is between $0$ and $1$, the risk-aversion coefficient is typically assumed to be between $0$ and $40$. For less obvious cases, it can be useful to estimate a simpler or reduced-form model to get a gross approximation of the parameter's magnitude. A good example is salanie2019 who approximate non-linear demand coefficients using a Taylor expansion and then use two-stage least-squares to approximate the coefficients.
Given a set of bounds for each coefficient, write the parameters space as $\Theta = [\underline{\theta}_1,\overline{\theta}_1] \times \dots \times [\underline{\theta}_{d_\theta},\overline{\theta}_{d_\theta}]$. In the special case where $\Theta = [0,1] \times \dots \times [0,1]$, there are sequences for which the dispersion $D(\theta^1,\dots,\theta^b)$ is close to its theoretical lower bound, these are called low-dispersion sequences. This is discussed in more detail in the next Section. The Sobol and Halton sequences are two popular low dispersion sequences judd1998,Lemieux2009. Both are readily available in R, Julia, Matlab, and Python. Finally, to construct the covering sequence for $\Theta$, take one of the Sobol or Halton sequence, denoted $\vartheta^b$, and adjust each of its coordinates $\theta^b_j = \underline{\theta}_j + [\overline{\theta}_j - \underline{\theta}_j] \vartheta_j^b$ to match the desired bounds.
Several factors can indicate that the user-specified $\Theta$ is too small to contain the solution $\hat{\theta}_n$. 1) If the local step exits $\Theta$ in some direction, and the global step cannot improve the fit inside $\Theta$, then the bounds could be widened in that direction. 2) If, after many iterations, $Q_n(\tilde{\theta}_n)$ is significantly non-zero, the bounds may be too narrow and could be enlarged. It may also be that the model is misspecified, a case that is not covered in the result.
\paragraph{Computing the Jacobian.} First, the following Lemma applies a change of variable argument similar to powell1989 or hazan2016 to derive some useful identities.
While in theory the smoothing can be applied to any estimation problem, $\overline{g}_{n,\varepsilon}$ and $G_{n,\varepsilon}$ typically do not have closed-form. This is the case for the applications in Section (ref). Nevertheless, an unbiased Monte Carlo estimator of $G_{n,\varepsilon}$ can be computed as follows:
using random draws $Z_\ell \overset{iid}{\sim} \mathcal{N}(0,I_{d_\theta})$ with $L \geq 1$. The mean-zero adjustment $\overline{g}_{n,\varepsilon}(\theta)$ yields better finite-$L$ properties. Other choices of densities $\phi$ could be used as long as they are: 1) continuously differentiable, and 2) satisfy $\int Z \phi(Z)dZ = 0$, $\int \|Z\|^2 \phi(Z) dZ < \infty$, $\int \|\partial_Z \phi(Z)\| dZ < \infty$. The Monte Carlo approximation $\hat{G}_{n,\varepsilon}$ does require $\phi$ to be a proper distribution, which excludes higher-order kernels.
A similar estimator was proposed in polyak1987, polyak1990 for minimizing smooth and globally convex objectives by stochastic gradient descent (sgd), and was later studied in nesterov2017. Convergence for gradient descent is slower than Gauss-Newton. chen2015 proposed a similar Jacobian for sieves. In an extension of the main results, a more computationally attractive quasi-Newton implementation of $\hat{G}_{n,\varepsilon}$ is introduced. Its non-asymptotic local convergence properties are derived. All simulated and empirical examples in the paper are based on this implementation. When the sample Jacobian $G_n$ exists, and satisfies appropriate regularity conditions, it can be used in the Local Step of Algorithm (ref) instead of $G_{n,\varepsilon}$.
\paragraph{Regularization.} Away from the solution, the Jacobian $G_{n,\varepsilon}$ could be nearly singular. This is a concern since inverting an ill-conditioned matrix can be inaccurate and unstable. In this case, it is recommend to invert $G_{n,\varepsilon}(\theta_b)^\prime W_n G_{n,\varepsilon}(\theta_b) + \lambda I_d$ instead, with $\lambda \geq 0$, in the local step to stabilize it. This is know as the Levenberg–Marquardt (LM) algorithm. It can be particularly useful when $G_{n,\varepsilon}$ is evaluated by Monte Carlo and/or quasi-Newton to reduce sensitivity of the inverse to the approximation error. Section (ref) gives specific recommendations for setting the penalty $\lambda_b$ at each iteration it that setting.
\paragraph{Stopping Criteria.} Standard convergence criteria for global optimization include that: a maximum number of iterations has been reached, or $Q_n$ has not declined more than some tolerance level (e.g. $\text{tol} = 10^{-6}$) for a certain number of iterations (e.g. $100$). When these occur, the optimization terminates. In the applications, Algorithm (ref) was terminated when a maximal number of iterations was reached, to keep the computation time fixed.
andrews1997 proposes an asymptotically valid data-dependent stopping criterion for GMM specifically. The following sketches the idea in the context of Algorithm (ref). For simplicity, suppose $W_n$ is the optimal weighting matrix and standard regularity conditions hold. Take $c_{p,n}$ to be the $1-\alpha$ quantile of a $\chi^2_{p-d_\theta}$ distribution, where $p$ is the number of estimating equations, $p-d_\theta$ is the number of overidentifying restrictions. andrews1997 suggests $\alpha = 5\%$ and tabulates values. One can also use the quantiles of a $\chi^2_{p}$ distribution, which is less conservative. Suppose at iteration $k \geq 1$ that $Q_n(\theta_k) \leq c_{p,n}/n$. Then, under regularity conditions, $\theta_k$ is a consistent estimator with $\| \theta_k - \theta^\dagger \| \leq O_p(n^{-1/2})$ and $\| \theta_k - \hat{\theta}_n \| \leq O_p(n^{-1/2})$. From there, run the algorithm for an additional $j \geq 1$ iterations. The algorithm's properties roughly imply that $\|\theta_{b} - \hat{\theta}_n\| \leq (1-\overline{\gamma})^j \| \theta_k - \hat{\theta}_n \|$ for $b=k+j$ and a $\overline{\gamma}$ positive that is smaller than $\gamma$. Then applying one of the two stopping criteria above gives a valid estimator, using $j_{\max} \gg \log(n)$ additional iterations.
The following illustrates several aspects of Algorithm (ref), including differences between the local step and smoothing the objective function as well as convergence issues under non-convexity. Consider estimating the $t$-th quantile of $x_i \overset{iid}{\sim} F$, $t\in (0,1)$. The population and sample moments are $F(\theta)-t=0$, and: \[\overline{g}_n(\theta) = \frac{1}{n} \sum_{i=1}^n \left[\mathbbm{1}_{x_i \leq \theta} -t\right] = 0, \] or more compactly $F_n(\theta)-t=0$; where $F_n$ is the empirical CDF of $x_i$. For any choice of bandwidth $\varepsilon$, the smoothed moment condition and its Jacobian are:\footnote{This is derived using $\mathbb{E}_{Z}(\mathbbm{1}_{x_i \leq \theta + \varepsilon Z}) = \Phi([\theta-x_i]/\varepsilon)$.} \[\overline{g}_{n,\varepsilon}(\theta) = \frac{1}{n} \left[\sum_{i=1}^n \Phi\left(\frac{\theta-x_i}{\varepsilon}\right)-t\right],\quad G_{n,\varepsilon}(\theta) = \frac{1}{n\varepsilon} \sum_{i=1}^n \phi\left( \frac{\theta-x_i}{\varepsilon} \right), \] where $\Phi$ and $\phi$ are the standard Gaussian CDF and density. Notice that $G_{n,\varepsilon}(\theta)=f_{n,\varepsilon}$ is the kernel density estimator of the density $f$. In this example, the local step in Algorithm (ref) involves the non-smooth CDF $F_n$ and the smoothed density $f_{n,\varepsilon}$: \[ \theta_{b+1} = \theta_b - \gamma \frac{ F_n(\theta_b)-t}{ f_{n,\varepsilon}(\theta_b) }. \] The main appeal is that an exact solution $\hat\theta_n$ with $F_n(\hat\theta_n)-t=0$ is a fixed-point of these iterations, regardless of the bandwidth $\varepsilon$. This implies that if the Algorithm is convergent, then the solution is robust to the choice of bandwidth.
Applying Gauss-Newton to the smoothed moments further substitutes $\overline{g}_n$ for $\overline{g}_{n,\varepsilon}(\theta) = F_{n,\varepsilon}(\theta) - t$, $F_{n,\varepsilon}(\theta) = \int_{-\infty}^\theta f_{n,\varepsilon}(u)du$ is the smoothed CDF. This results in the iterations: \[ \theta_{b+1} = \theta_b - \gamma \frac{ F_{n,\varepsilon}(\theta_b)-t}{ f_{n,\varepsilon}(\theta_b) }. \] Here, a fixed-point $\hat\theta_{n,\varepsilon}$ must satisfy $F_{n,\varepsilon}(\hat\theta_{n,\varepsilon})-t=0$. This is the $t$-th quantile of $F_{n,\varepsilon}$ which satisfies $\hat\theta_{n,\varepsilon} = \hat\theta_n + O_p(\varepsilon^2)$. The smoothing bias, captured by $\varepsilon^2$, is only negligible when $\varepsilon = o(n^{-1/4})$ which is more restrictive.
Although estimating a quantile is simple, the GMM formulation can be non-convex. For simplicity take $x_i \sim \mathcal{N}(0,1)$ and, to focus on non-convexity specifically, consider the smooth population loss $Q(\theta) = [ \Phi(\theta) - t ]^2$ with $t = 0.7$. Numerical calculations show that the Hessian $\partial^2_{\theta} Q(\theta)$ is positive for $\theta \in [-0.69,1.15]$, there $Q$ is locally convex, but is negative for $\theta \in [-3,-0.7]$ and $\theta \in [1.16,3]$, there $Q$ is locally concave. As a result, an off-the-shelf BFGS optimizer fails to minimize the population $Q$ from $\theta = -3$, it returns $\theta = 5.02$ instead of $\theta^\dagger = 0.52$. The sample problem is also non-convex.
Even in this basic example, local optimization does not necessarily converge from an arbitrary value. To resolve this: the global step will produce, if the covering sequence has low dispersion, a value $\theta^b$ that is sufficiently close to $\theta^\dagger$ from which the local step is convergent. From that iteration onwards, the local step is convergent. Hence, the global step ensures the optimizer converges from an arbitrary starting value.
Figure (ref) illustrates the difference between Algorithm (ref) and smoothing the GMM objective. It compares estimates based on R's quantile function\footnote{The sample quantile function in R returns an estimate obtained by a linear interpolation between order statistics, see e.g. hyndman1996 for details. As such, it does not exactly solve the moment condition $\overline{g}_n(\theta)=0$, and its finite sample properties can differ slightly from the GMM estimator.} with both methods. Estimates from Algorithm (ref) and the sample quantiles are nearly identical, for all $\varepsilon$. The blue curve corresponds to $\hat\theta_{n,\varepsilon}$, its bias is visible when $\varepsilon = 0.5,0.2$. Since $n^{-1/4} = 0.25$, the choice $\varepsilon = 0.1$ represents the undersmoothing regime. The bandwidth minimizing the MSE criterion is also reported under Auto, it generally close to the smaller value, i.e. $\varepsilon \simeq 0.13$. Table (ref) further illustrates the bias of the estimates and its implications for inference. Inferences using $\hat\theta_{n,\varepsilon}$ are size distorted compared to sample quantiles and Algorithm (ref). The size of the test also depends on the choice of bandwidth. The table also reports the size when standard errors are computed as in Algorithm (ref) (size$_1$) or, as a benchmark, using the analytical expression and Silverman's rule of thumb to evaluate the unknown density term (size$_2$).
In this example, smoothing replaces the indicator function with the same $\Phi(\cdot/\varepsilon)$ for all $i$. This is generally not the case. Consider a quantile regression with moments $\overline{g}_n(\theta) = \frac{1}{n} \sum_{i=1}^n [\mathbbm{1}( y_i - x_i^\prime \theta )-t]w_i$, $w_i$ are the instruments. The smoothed moments are $\overline{g}_{n,\varepsilon}(\theta) = \frac{1}{n} \sum_{i=1}^n [\Phi( \frac{y_i - x_i^\prime \theta}{\varepsilon\|x_i\|_2} )-t]w_i$. Here, the bandwidth is effectively $\varepsilon \|x_i\|_2$ which is individual specific. Similarly, the smoothed Jacobian is $G_{n,\varepsilon}(\theta) = \frac{-1}{n \varepsilon} \sum_{i=1}^n \frac{x_i^\prime}{\|x_i\|_2}\phi( \frac{y_i - x_i^\prime \theta}{\varepsilon\|x_i\|_2} )w_i$ which, compared to smoothing the indicator function directly, has the additional $\|x_i\|_2$ term in the denominator. Notice that all moments of $x_i/\|x_i\|_2$ are finite even when $x_i$ has infinite moments.
This section gives the main assumptions and derives the main results. The finite sample convergence results are derived. The role of the tuning parameters is discussed and the role of covering sequences $(\theta^b)_{b \geq 0}$ on the global convergence properties is explored in more details.
Assumption (ref) are the conditions required for the population quantities. Conditions ii, iii, vi imply $Q(\theta) = \|g(\theta)\|^2_W$ behaves like $\|\theta-\theta^\dagger\|^2$ around the solution $\theta^\dagger$. Conditions ii, iv hold in the particular case where $g(\cdot)$ is twice continuously differentiable with bounded second derivative. Condition iii and v are textbook local and global identification assumptions. These two conditions are particularly important for the results. In practice, it is important to check if the two conditions hold since they determine whether the optimization is well-posed and affect the speed of convergence. Often, this is difficult to verify analytically; forneron2024 provides a numerical procedure to detect whether these conditions fail.
Assumption (ref) gives additional conditions for the sample moments. The assumptions allow for a variety of estimation problems, including simulated method of moments estimations with non-smooth moments. As an example, consider a simulated method of moments estimation of a static discrete choice model using: \[ g(\theta;x_i) = \frac{1}{S} \sum_{s=1}^S z_i ( y_i - \mathbbm{1} \{ w_i^\prime \theta + u_i^s \geq 0 \} ), \] where $x_i = (y_i,z_i,w_i,u_i^1,\dots,u_i^S)$, $y_i$ is a binary outcome, $z_i$ are instruments, $w_i$ regressors, $u_i^s$ are simulation draws and $S \geq 1$ the number of simulated samples. For a given value of $\theta$, $\mathbbm{1} \{ w_i^\prime \theta + u_i^s \geq 0 \}$ is a simulated discrete outcome. For this model, Assumption (ref) i holds if $z_i$ has finite variance. Condition ii holds with $\psi = 1/2$ under additional moment conditions for $w_i,z_i$ and if $u_i^s$ has a continuous density. The GMM quantile regression discussed above also falls within this class. Other examples include moments that satisfy what is known as a stochastic Lipschitz condition: \[ \| g(\theta_1;x_i) - g(\theta_2;x_i) \| \leq L(x_i)\|\theta_1 - \theta_2\|, \text{ for any } \theta_1,\theta_2 \in \Theta, \] where $L(x_i)$ has finite variance. This holds in the special case where $g(x_i;\theta)$ is smooth: i.e. continuously differentiable and its derivative is bounded by $L$, i.e. $\|\partial_\theta g(x_i;\theta)\| \leq L(x_i)$ for all $\theta$. For these models, condition i holds since $\Theta$ is bounded, and ii holds with $\psi = 1$. Assumption (ref) i, ii is related to the so-called stochastic equicontinuity property. pakes1989, andrews1994, newey-mcfadden-handbook list additional examples, including models which feature censoring, sample selection, simulated multinomial and dynamic discrete choice, and further discuss the so-called $L_2$-smoothness condition ii.
Under Assumptions (ref)-(ref), it can be shown that the approximate minimizer $\hat\theta_n$ is $\sqrt{n}$-consistent and asymptotically Gaussian: \[ \sqrt{n}(\hat\theta_n - \theta^\dagger) \overset{d}{\to} \mathcal{N}(0,V), \] where $V$ is the asymptotic variance of $\hat\theta_n$, see e.g. newey-mcfadden-handbook.
Given Assumption (ref), (ref) above, the following first provides large sample convergence properties of Algorithm (ref). Then, several finite sample results leading to large sample ones are stated; they indicate which rate is optimal for the smoothing parameter $\varepsilon$. Results for the smooth population objective are given in Appendix (ref), as a baseline without sampling uncertainty and non-smoothness. The results below are stated in terms of an approximate minimizer:
where $G = G(\theta^\dagger)$, it can be shown to satisfy ((ref)).
Theorem (ref) ((ref)) implies that after $b=k+j$ iterations, the iterate $\theta_b$ from Algorithm (ref) is an approximate minimizer of the GMM objective and ((ref)) implies that it is a $\sqrt{n}$-asymtotically normal estimator of the parameters of interest $\theta^\dagger$. As in the last step of Algorithm (ref), given an appropriate estimator $\Sigma_n$ of $\Sigma$, the smoothed Jacobian $G_{n,\varepsilon}(\tilde{\theta}_n)$ can be used to compute standard errors and make inferences on the true $\theta^\dagger$.
There are two phases in the convergence which contribute to the required number of iterations $b = k+j$. First, the Algorithm performs a global search. For a well-chosen covering sequence, $D(\theta^0,\dots,\theta^b) \to 0$ implies that there exists $k$ such that the global step necessarily produces a $\theta^\ell$, $0 \leq \ell \leq k$, that is close enough to $\hat{\theta}_n$ for 1) the local Gauss-Newton step to be convergent, and 2) there are no value away from $\hat{\theta}_n$ that can improve the fit -- i.e. if a $\theta^b$ lowers the objective value in the global step, with $b \geq k$, then it is also locally convergent. This implies that the Algorithm remains local after at most $k$ iterations. Importantly, 1) hinges on the local identification condition, and 2) on the global identification condition (respectively Assumptions (ref) iii. and v.).
After those $k$ iterations, Algorithm (ref) automatically transitions into a fast local convergence phase. Only an additional $j$ iteration, increasing slowly with $n$, are required to produce asymptotically valid estimates that satisfy ((ref)) and ((ref)). The main bottleneck of the Algorithm is typically the finite number of iterations $k$, which depends on the GMM objective and the choice of covering sequence. The latter is now discussed in more detail before presenting several finite sample local convergence results leading to the main Theorem above.
\paragraph{Choice of covering sequence.} What is particularly attractive in this result is the automatic transition from $k$ global steps to the $j$ local ones without user input. Because $k$ only depends on the GMM objective and $(\theta^\ell)_{\ell \geq 0}$, the Algorithm adapts to the complexity of the optimization problem in a single run. Notice that this property comes from using the local and global steps concurrently. Methods that use them successively, as in e.g. robinson1988, andrews1997 do not have this property as explained in Appendix (ref).
To get good convergence results, the covering sequence should reach a neighborhood of $\hat{\theta}_n$ using as few iterations $k$ as possible. This is measured by dispersion $D$ considered in the choice of covering sequence above. For each $k \geq 1$, covering number arguments imply that there exists $(\theta^0,\dots,\theta^k)$ which minimize the dispersion with $D(\theta^0,\dots,\theta^k) = O(k^{-1/d_\theta})$. The problem is that, for $k+1$, the optimal set $(\theta^0,\dots,\theta^k,\theta^{k+1})$ does not use the previous $(\theta^0,\dots,\theta^k)$. The optimal cover for a given $k$ does not necessarily define a proper sequence over $k$. Importantly, the ordering of the points in the sequence matters to ensure the discrepancy $D(\theta^0,\dots,\theta^k)$ is as small as possible for each $k \geq 1$.
In the special case where $\Theta = [0,1]^{d_\theta}$, niederreiter1983 showed that low-discrepancy sequences, such as the recommended Sobol or Halton sequence, satisfy: \[ D_{k-1} \leq O( \sqrt{d_\theta} \log(k)k^{-1/d_\theta}), \] which is close to the theoretical lower bound $O(k^{-1/d_\theta})$ at each $k$. In comparison, his result implies that uniform random draws $\theta^j \overset{iid}{\sim} \mathcal{U}_{\Theta}$ only achieve $D_{k-1} = O_p(k^{-1/[2d_\theta]})$, which is slower. Differences between random and low dispersion sequences are illustrated in Figure (ref), Appendix (ref). Finally, although $\Theta$ needs to be bounded to have $k < \infty$, it need not be finite-dimensional. For instance, gallant1987 prove compactness of the Hermite polynomial density space. This implies existence of a finite cover, a sufficient condition for $k < \infty$. Finite covers also exist for the class of smooth functions on bounded convex sets VanderVaart1996. Constructing a covering sequence with low-dispersion in these settings could be challenging, however. Additionally, chen2015 consider smoothing of the Jacobian in an non-parametric setting, as required for the local step. Extensions to the infinite-dimensional setting are not considered here.
\paragraph{Local Convergence.} The following considers the convergence properties of local step only:
where $G_b = G_{n,\varepsilon}(\theta_b)$. As noted earlier, the Gauss-Newton updates rely on a smoothed Jacobian to approximate an infeasible derivative when the problem is non-smooth. Local convergence requires $G_{n,\varepsilon}(\hat{\theta}_n)^\prime W_n \overline{g}_n(\hat{\theta}_n)$ to be as small as possible, in order to behave like a first-order condition for the minimization problem. To this end, Lemma (ref) below measures the stability of ((ref)) at $\theta=\hat\theta_n$.
The term $\Gamma_{n,\varepsilon}$ measures how close $\hat\theta_n$ is to being a fixed-point of ((ref)).\footnote{In just-identified models for which an exact solution $\overline{g}_n(\hat\theta_n)=0$ exists, it is possible to set $\Gamma_{n,\varepsilon}=0$ by using the exact minimizer.} Here $c_n \geq 1$ controls the probability level for which the bound holds over repeated samples, the bound is informative for $c_n \geq 1+C$. It will be useful for understanding the relation between the probability of a successful optimization and the choice of tuning parameters. The bandwidth $\varepsilon$ appears in $\Gamma_{n,\varepsilon}$ in both a numerator and a denominator. The optimal bandwidth for minimizing $\Gamma_{n,\varepsilon}$ is $\varepsilon \asymp (c_n n^{-1/2})^{(1+\psi)/2}$, i.e. $\varepsilon = O(n^{-1/2})$ for the continuous case ($\psi=1$), and $\varepsilon = O(n^{-3/8})$ for the discontinuous one ($\psi=1/2$). Using smoothed instead of sample moments in Algorithm (ref) would require to substitute $\overline{g}_n(\hat\theta_n)$ for $\overline{g}_{n,\varepsilon}(\hat\theta_n)$ in Lemma (ref) above. Then, two additional terms would appear in $\Gamma_{n,\varepsilon}$, respectively of order $\varepsilon^2$ and $\sqrt{\varepsilon}n^{-1/2}$. In particular, $\varepsilon^2$ would add smoothing bias to the estimates.
In the following, $\sigma_{\min}$ denotes the smallest singular value, i.e. $\sigma_{\min}[G(\theta)] = \sqrt{ \lambda_{\min}[G(\theta)^\prime G(\theta)] }$ where $\lambda_{\min}$ returns the smallest singular value of the square matrix $G(\theta)^\prime G(\theta)$.
Proposition (ref) gives the local convergence properties of a single iteration of the local step. Recall that the local step relies on $G_{n,\varepsilon}$ to linearly approximate the moments around $\theta_b$. Since the moments are non-smooth, the approximation does not become more accurate as $\theta$ approaches $\theta_b$. The approximation does become more accurate as the sample size $n$ increases and the bandwidth $\varepsilon$ decreases. With this in mind, the result implies that, with high probability, the local step produces a $\theta_{b+1}$ which reduces the distance to $\hat{\theta}_n$ by a factor $(1-\overline{\gamma})$ up to an approximation error, measured by $\Delta_{n,\varepsilon}$ which involves non-smoothness of the moments and smoothing bias in the Jacobian. Note that, unlike e.g. dennis1996 in the smooth case, the derivations do not rely on local convexity arguments which would define a different local neighborhood. For instance, the quantile problem $[\Phi(\theta)-t]^2$ is globally non-convex. Nevertheless, for the quantile example there always exists a $\underline{\sigma}$ such that $R_G = \Theta$ when it is bounded since the density $\phi$ is strictly positive everywhere.
The local convergence neighborhood is of size $\min(R_{n,G},R_{n,\varepsilon})$ which depends on several factors: features of the moments $R_G$, $\underline{\sigma}$, $L_G$, $\kappa_W$, and tuning parameters $\varepsilon,\gamma$.\footnote{The constants $C_a,C_{\sigma},M_{1,Z},C_2$ are made more explicit in Proposition (ref) in the Appendix.} $R_G$ measures the area where the singular values of $G$ are a least $\underline{\sigma}$, a measure of local identification. $L_G$ measures the non-linearity of the moments; $L_G=0$ implies $g(\cdot)$ is linear in $\theta$. $\kappa_W$ bounds the condition number of $W_n$; setting $W_n = I_d$ yields $\kappa_W=1$. A larger $L_G,\kappa_W$ and/or a smaller $R_G,\underline{\sigma}$ makes the area of convergence smaller. Very non-linear models (large $L_G$) or less sharply identified (small $R_G,\underline{\sigma}$) are harder to estimate since the area of local convergence is reduced. Using the optimal weighting matrix is more efficient from a statistical perspective. However, if it is poorly conditioned (large $\kappa_W$) then estimation is more difficult. This motivates less efficient choices such as $W = (\Sigma + \lambda)^{-1}$ where $\Sigma = \text{var}[g(x_t;\theta^\dagger)]$ and $\lambda > 0$ is a regularization parameter which brings $\kappa_W$ closer to $1$.
The bandwidth $\varepsilon$ also affects the area of convergence via $R_{n,\varepsilon}$. If $\varepsilon \to 0$ and $\sqrt{n}\varepsilon \to 0$, then $\min(R_{n,G},R_{n,\varepsilon}) = R - o(1)$, where $R = \min(R_G,[\gamma-\overline{\gamma}]\underline{\sigma}/[\gamma L_G \sqrt{\kappa_W}])$ is the area of local convergence in the smooth population problem (see Lemma (ref)). A larger learning rate $\gamma$ allows for faster convergence, measured by $\overline{\gamma} < \gamma$, but over a smaller neighborhood. For smooth problem, it is possible to have $\gamma \nearrow 1$ as $\theta_b$ approaches $\hat{\theta}_n$ to make convergence faster. Here, a larger $\gamma$ increases the effect of $\Delta_{n,\varepsilon}$ which can make local optimization less stable. Hence, the recommendation to use a smaller $\gamma$.
The optimal choice of bandwidth, for maximizing the area of local convergence $R_{n,\varepsilon}$, is $\varepsilon = O( n^{-1/4} )$. This choice of bandwidth is larger than for minimizing $\Gamma_{n,\varepsilon}$ above. Lastly, a smaller value $\underline{\sigma}$ makes Algorithm (ref) more sensitive to the choice of $\varepsilon$ given its role in $\underline{\sigma}_{n,\varepsilon}$. The optimal choice for maximizing $\underline{\sigma}_{n,\varepsilon}$ is also $\varepsilon = O( n^{-1/4} )$.
While Proposition (ref) derived properties for a given iteration of the local step, Theorem (ref) describes the full optimization path $(\theta_0,\theta_1,\dots)$ from a local starting value $\theta_0$. The final step of Algorithm (ref) returns $\tilde \theta_n$ such that $\|\overline{g}_n(\tilde\theta_n)\|_{W_n} = \min_{j=0,\dots,b} \|\overline{g}_n(\theta_j)\|_{W_n} \leq \|\overline{g}_n(\theta_b)\|_{W_n}$.
The Theorem states that only $b = O(|\log[\Gamma_{n,\varepsilon}]|) = O(\log[n])$ iterations are needed to achieve $\|\theta_b-\hat\theta_n\| = O_p(\Gamma_{n,\varepsilon}) = o_p(n^{-1/2})$ if $\varepsilon \to 0$ and $\sqrt{n}\varepsilon \to 0$. The term $\Gamma_{n,\varepsilon}$ captures, among other things, the smoothing bias in the estimates. Using an optimal rate $\varepsilon = O(n^{-1/4})$ yields $\Gamma_{n,\varepsilon} = O_p(n^{-3/4})$ for both the continuous ($\psi=1$) and discontinuous case ($\psi=1/2$). This is negligible compared to sampling uncertainty, as required to prove ((ref)) in Theorem (ref).
In general, the local step alone does not guarantee convergence. Theorem (ref), Appendix (ref), extends Theorem (ref) to a global convergence result when using both global and local steps.
This section provides two extensions for the local step in Algorithm (ref) and Proposition (ref). The first extension adds momentum to ((ref)), a.k.a. the Polyak heavy-ball polyak1964. It is commonly used in stochastic gradient descent to accelerate convergence. Here it can be used to either accelerate convergence or maintain the rate of convergence while reducing the effect of sampling noise on optimization. The second extension builds a quasi-Newton approximation of $G_{n,\varepsilon}$ which is useful when the smoothed Jacobian does not have closed-form.
The first extension considers a modification of the local search step:
where $\theta_{-1} = \theta_0$, $G_b = G_{n,\varepsilon}(\theta_b)$ and $\alpha \in [0,1)$ is the momentum parameter. Two new quantities affect convergence, the companion matrix $A(\gamma,\alpha)$ and the effective rate $\gamma(\alpha)$: \[ A(\gamma,\alpha) = \left(
\right),\quad 1-\gamma(\alpha) = \sigma_{\max}[A(\gamma,\alpha)],\] For any $\gamma \in (0,1)$, there exists $\alpha \in (0,1]$ such that $\gamma(\alpha) > \gamma$, and a $\alpha^\star$ which maximizes $\gamma(\alpha)$. Table (ref) provides a selection of combinations $(\gamma,\alpha^\star)$. The last row measures $\gamma/\overline{\gamma(\alpha)} < 1$ which implies that setting $\overline{\gamma}/\gamma > 1$ is now possible: momentum can affect the convergence rate / sensitivity to sampling noise tradeoff in Proposition (ref).
Notice that $\theta_{b+1} - \alpha(\theta_b-\theta_{b-1}) = \theta_b - \gamma (G_b^\prime W_n G_b)^{-1}G_b^\prime W_n \overline{g}_n(\theta_b)$ is the same as ((ref)) in Proposition (ref), so convergence of ((ref)') follows from the same derivations. The results are given for $\boldsymbol{\theta_b} = (\theta_b^\prime,\theta_{b-1}^\prime)^\prime$ to denote $\theta_b$ and its lagged value $\theta_{b-1}$; similarly $\boldsymbol{\hat\theta_n} = (\hat\theta_n^\prime,\hat\theta_n^\prime)^\prime$.
The convergence is now stated in terms of $\boldsymbol{\theta_{b+1}}-\boldsymbol{\hat\theta_n}$, similar to a first-order vector autoregression. The main advantage of using momentum is that it can allow to maintain the rate of convergence but using a smaller learning rate $\gamma$. This can potentially result is a larger area of local convergence, without sacrificing speed of local convergence.
In many empirical applications, the smoothed Jacobian $G_{n,\varepsilon}$ is not available in closed form. The following proposes a computationally attractive approximation.
Algorithm (ref) involves an additional tuning parameter: $L \geq d_\theta$. A larger value of $L$ reduces the Monte Carlo error but involves more lagged values $\theta_{b-\ell}$, $\ell = 0,\dots,L-1$, which slows convergence. The algorithm only adds one moment evaluation per iteration, rather than $L \gg 1$ for a direct Monte Carlo implementation. The theoretical results require $L \to \infty$, but the divergence can be slow. Updating several directions at each $b$ could speed-up convergence as it improves the approximation quality.\footnote{Updating several directions in parallel, e.g. $2$ or $4$, can improve the finite-sample optimization properties because, for the same choice of $L$, $\hat{G}_b$ depends on fewer lagged values, e.g. $L/2$ or $L/4$.} For simplicity, the results are only derived for the sample mean estimator: $\hat G_L(\theta_b) = \frac{1}{L} \sum_{\ell=0}^{L-1} \frac{1}{\varepsilon} [ \overline{g}_n(\theta_{b-\ell} + \varepsilon Z_{b-\ell}) - \overline{g}_n(\theta_{b-\ell}) ]Z_{b-\ell}^\prime$. The following extends Proposition (ref) when $\hat{G}_{L}(\theta_b)$ is used in ((ref)). It is assumed that the Algorithms runs for at most $1 \leq b_{\max} < \infty$ iterations (a stopping criterion).
Compared to Proposition (ref), the number of iterations $b_{\max}$ now plays a role in the results. This is because $\hat{G}_{L}(\theta_b)$ is computed by Monte Carlo and the worst-case approximation error can be made arbitrarily large over infinitely many iterations. The dependence is only logarithmic because $\hat G_L(\theta_b) \simeq \frac{1}{L} \sum_{\ell=0}^{L-1} \partial_\theta g(\theta_{b-\ell})Z_{b-\ell} Z_{b-\ell}^\prime,$ where $Z_{b-\ell} Z_{b-\ell}^\prime$ are iid Wishart distributed. Relying on properties of this distribution, the proof derives exponential tail bounds for the singular values of $\hat{G}_b$. This allows $b_{\max}$ to increase linearly or polynomially in $n$ so that the Monte Carlo error term $\delta_n$ only diverges logarithmically.
The probability bound is now $1-(5+C)/c_n$ to account for the Monte Carlo error. Comparing $\Delta_{n,\varepsilon}$ with $\hat{\Delta}_{n,\varepsilon}$ in Propositions (ref), (ref) highlights the importance of the lags, measured by $\mathcal{E}_b$, in the convergence. The term $\hat{\Gamma}_{n,\varepsilon}$ depends on $(nL)^{-1/2}$ so that $L \to \infty$ is required. This divergence can, in theory, be arbitrarily slow. In practice, setting $L \geq \max(25,1.5 \times d_\theta)$ typically yields good results.
The proof of the Proposition highlights several factors which can affect convergence when using $\hat{G}_L$. Algorithm (ref) relies on Monte Carlo draws $Z_{b-\ell}$ and past updates $Y_{b-\ell}$, two sources of approximation. If the matrix $G_{n,\varepsilon}^\prime W_n G_{n,\varepsilon}$ is ill-conditioned, because $\underline{\sigma}$ is small, the inverse $(\hat{G}_L(\theta_b)^\prime W_n \hat{G}_L(\theta_b))^{-1}$ can be sensitive to Monte Carlo noise. This is measured by $\hat{\sigma}_{n,\varepsilon}$ in the Proposition. Also, if $\theta_b$ changes rapidly, $\theta_{b-\ell} - \theta_b$ is large and past values $Y_{b-\ell}$ give a poor approximation of the current $G_{n,\varepsilon}(\theta_b)$ when the underlying moments $\overline{g}_{n}$ are very non-linear.
\paragraph{Regularization.} To counter these two issues, the Levenberg-Marquardt (LM) algorithm which inverts $(\hat{G}_L(\theta_b)^\prime W_n \hat{G}_L(\theta_b) + \lambda_b I_d)$ instead of $(\hat{G}_L(\theta_b)^\prime W_n \hat{G}_L(\theta_b))$ in the local step can be particularly useful. Here, $\lambda_b$ penalizes for the quality of the approximation $\hat{G}_L(\theta_b) - G_{n,\varepsilon}(\theta_b)$. The applications use $\lambda_b = \tilde{\lambda} \textsc{rmse}_b \times \|W \hat{G}_{L}(\theta_b)\|$ where $\textsc{rmse}_b$ is the root mean squared error in the regression of $Y_{b-\ell}$ on $\tilde{Z}_{b-\ell}$ and $\tilde{\lambda} = 10^{-3}$. This is implemented using R in Appendix (ref). When the moments are linear, $\textsc{rmse}_b = 0$ and $\lambda_b = 0$ for all $b$. When the moments are more non-linear, $\textsc{rmse}_b$ and $\lambda_b$ are larger. This shrinks the size of the update $\theta_{b} - \theta_{b-1}$. Increased stability comes at the cost of slower convergence. In simulated examples below, this also helps stabilize the local step when $\varepsilon$ is chosen to be very small.
The first simulated example adapts an example from bruins2018. The model is a simple panel dynamic discrete choice model. The data generating process (dgp) is given by:
where $x_{it}$ are iid, strictly exogenous regressors, $i=1,\dots,n$ and $t=1,\dots,T$. The parameter of interest is $\theta = (\beta,\rho)$. The model is estimated by matching OLS estimates from a linear regression of $y_{it}$ on $x_{it}$ and $y_{it-1}$ between sample and simulated datasets. Because of the indicator function $\mathbbm{1}$, the objective function is discontinuous in $\theta$. The simulations use $n=250$, $T=10$, $\text{dim}(\beta) = 14$, $\beta^\dagger_{1}=\dots=\beta^\dagger_{5} = 1/\sqrt{5}$, $\beta^\dagger_{6}=\dots=\beta^\dagger_{14}=0$, $\rho^\dagger = 0.7$. The distant starting value is $\theta_0 = (0,\dots,0)$.
Figure (ref) compares Algorithm (ref) with a local BFGS algorithm with smoothed moments, replacing $\mathbbm{1}\{\cdot > 0\}$ with $\Phi(\cdot / \varepsilon)$ in the dgp, where $\Phi$ is the CDF of the standard Gaussian. It also reports results for the unsmoothed problem using a battery of global optimizers: Particle Swarm Optimization clerc2010, Covariance matrix adaptation evolution strategy hansen2003, Dividing Rectangles jones1993, Multi-Level Single-Linkage with random numbers and low-dispersion sequence rinnooy1987, and Simulated Annealing followed by Nelder-Mead (sa-nm), the latter is commonly used in Economics. These methods apply to continuous and discontinuous problems. \textsc{direct}, Simulated Annealing and \textsc{mlsl} have global convergence guarantees for continuous problems, given appropriate tuning belisle1992,locatelli1998. Similar to Algorithm (ref), \textsc{mlsl} and \textsc{mlsl-lds} are hybrid methods. It is a type of multi-start algorithm which combines random draws or the Sobol sequence, respectively, with a local search algorithm. At each iteration, one or several local searches are run until convergence. See Appendix (ref) for convergence properties of multi-start Algorithms relative to Algorithm (ref). \textsc{cma-es} has convergence guarantees for a class of scale-invariant functions.\footnote{$Q$ is scale-invariant if $Q(\theta^\dagger + \Delta_1) \leq Q(\theta^\dagger + \Delta_2) \Leftrightarrow Q(\theta^\dagger + \rho \Delta_1) \leq Q(\theta^\dagger + \rho \Delta_2)$ for any $\rho >0$. See toure2023 for the results and their Assumptions F1, F2 for the conditions on $Q$.} \textsc{pso} is heuristic, i.e. does not have formal convergence guarantees.
The top panel compares the final value of the objective $Q_n$ among different methods and tuning parameters. BFGS only converges when initialized near the true value, the computational cost increases when the starting value is further away and when the bandwidth becomes smaller. Global optimizers do not perform as well despite a heavy computational cost (bottom panel). Algorithm (ref) is run using the distant starting value and the Sobol sequence. It is fairly robust to the choice of bandwidth. The MSE minimizing bandwidth (Auto) gives satisfying result, the time reported in the bottom panel includes the bandwidth selection step. Although mlsl-lds also relies on the Sobol sequence with a local search, its performance is much worse than Algorithm (ref).
Table (ref) compares estimates between BFGS with a smoothed indicator and Algorithm (ref). The table reports the estimates' $\tilde{\theta}_n$ bias as well as the mean absolute difference between $\tilde{\theta}_n$ and $\hat{\theta}_n$, found using Algorithm (ref) initialized at $\theta^\dagger$, run for $1250$ iterations with a small bandwidth. The mean absolute error is normalized by the standard deviation of the estimates to be on the same scale as a critical value for a t-test, e.g. $1.96$. Algorithm (ref) is generally close to $\hat{\theta}_n$ whereas smoothed GMM is more sensitive to the choice of bandwidth. For the smallest $\varepsilon = 0.005$, the Levenberg-Marquardt algorithm stabilizes the local step significantly but slows convergence so that $250$ iterations do not systematically produce a valid estimator.
The second simulated example compares the estimation of a simple heterogeneous agent model. These models are key for evaluating the distributional effects of macroeconomic policies. They are, however, very difficult to estimate. As a result, many report calibrated, rather than estimated, results. The following illustrates the properties of Algorithm (ref) for an SMM estimation of a textbook Aiyagari model. The parameters of interest are $\theta=(\beta,\gamma,\mu,\rho,\sigma)$ which are the discount factor, risk-aversion, average log-income, persistence of log-income, and the standard deviation of log-income shocks. This type of estimation is challenging because the model is discretized, approximately solved by value function iterations, and the $5$ parameters are estimated using $10$ non-smooth moments (quantiles). The model and solution method are described in more detail in Appendix (ref).
Table (ref) compares the minimized objective values, parameter estimates, and computation time for one simulated sample with a range of bandwidths $\varepsilon \in \{0.5,0.2,0.1,0.05,0.01\}$, $\varepsilon = 0.1$ is the MSE minimizing bandwidth.\footnote{ The LM algorithm is applied to improve stability of the local step over the range of tuning parameters. } and the global optimizers considered in the previous example.\footnote{The global optimizers have a maximum of $50000$ function evaluations, to keep runtime below $6$ hours.} The table suggests that $\gamma$ is the hardest parameter to estimate, as several methods find reasonable estimates for the others. pso finds accurate estimates, except for the coefficient risk-aversion $\gamma$. Figure (ref) and Table (ref), Appendix (ref) illustrate the optimization paths for Algorithm (ref) and considers additional values for the tuning paramaters.
The empirical example replicates the estimation of the joint duration model of employment by honore2018. They model the joint decision of optimal retirement age for married couples. In their setting, husbands and wives engage in a Nash bargaining game to determine when each will retire. Unlike single-agent duration models, their specification does not have a closed-form likelihood so they consider simulation-based estimation. Because the simulated outcomes are discrete, the resulting objective function is discontinuous and fairly difficult to minimize. They use a repeated succession of 5 optimizers.\footnote{Their procedure: “The following loop of procedures was used until a loop produced a change in the parameter estimate of less than $10^{-5}$. [\dots] (a) particle swarm [\dots], (b) Powell’s conjugate direction method, (c) downhill simplex using Matlab’s fminsearch routine (d) pattern search using Matlab’s built-in routine, (e) particle swarm focusing on the jump-parameters [\dots].” See honore2018.} The following compares their estimation from an accurate initial guess with Algorithm (ref) initialized using a randomized Sobol sequence.
Table (ref) replicates the first two columns of Tables 2 and 3 in honore2018. For brevity, only the first four coefficients are reported here. The estimates are similar for both specifications. However, the time required to compute the estimates is significantly reduced: the two specifications take 3.5 and 5.5 hours to estimate with their replication code. In comparison, Algorithm (ref) finds estimates in less than 15 minutes using the same code. For comparison, for the $12$ coefficient model, the random-walk Metropolis-Hastings (MH) algorithm has yet to converge after $100000$ iterations (i.e. 10 hours, see Figure (ref), Appendix (ref)) for the same distant starting value. MH converges more quickly from the same starting value as honore2018 (not reported). Table (ref) in Appendix (ref) provides additional results for $\varepsilon = 0.8 \cdot 10^{-2}, 2.11 \cdot 10^{-2}$. The latter is the MSE minimizing bandwidth. The range of bandwidths that provide accurate estimates is much narrower than in previous examples. Two features explain this issue: first, the optimal weighting matrix $W_n$ used for estimation has condition number $\kappa = 1.9 \cdot 10^{3}$, $4.4 \cdot 10^{3}$, for the smaller and larger specification, respectively. This is because several moments are strongly correlated. Second, the Jacobian estimated using Algorithm (ref) at the solution has condition number $2 \cdot 10^{3}$, $5 \cdot 10^3$, indicating a relatively small $\underline{\sigma}$. Note that the dummies corresponding to $\ge$ 62 yrs-old and $\ge$ 65 yrs-old, and the corresponding moments, are somewhat redundant -- using a more parsimonious specification could help improve on both of these issues.
This paper shows, using an explicit algorithm, that faster GMM estimation is possible under standard econometric assumptions. For local optimization, the main idea is that smoothing only the Jacobian produces less smoothing bias -- allowing for larger bandwidths than existing methods which benefit optimization properties. For global optimization, introducing a global search step allows to adapt to the complexity of the estimation problem, without additional user intervention or input. \baselineskip=14.0pt