Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
116,306 characters · 22 sections · 44 citation commands
SLIM: Stochastic Learning and Inference in Overidentified Models
\onehalfspacing
Modern machine learning methods have revolutionized data analysis across numerous disciplines. Stochastic approximation, most prominently stochastic gradient descent (SGD), has emerged as a foundational tool behind the success of the recent deep learning revolution. First introduced by Robbins-Monro-1951, SGD is celebrated for its computational simplicity and scalable online implementation. It has become indispensable for modern optimization, particularly in settings involving a large number of parameters, huge datasets, and complex nonlinear objective functions bottou2018optimization. By iteratively approximating gradients using random samples, SGD achieves substantial computational efficiency, making it highly effective for large-scale optimization problems that arise in contemporary data analysis.
SGD has been widely applied in the context of M-estimation, where the parameter of interest \(\theta_o\) is defined as the minimizer of a population loss function \(\mathbb{E}[\ell(z_i, \theta)]\), with \(\ell(z_i, \theta)\) a known real-valued function of the observation \(z_i\) and parameter \(\theta\). The canonical form of the algorithm is
where \(t\) indexes the iteration count. This approach has proven effective for solving large-scale optimization problems, particularly in machine learning and statistics.
However, the application of SGD in econometrics remains limited. A central reason is that many parameters of interest in economics are not defined as minimizers of objective functions, but rather as solutions to systems of moment conditions: \[ \mathbb{E}[g(z_i, \theta_o)] = 0, \] where \(g(z_i, \theta)\) is a vector-valued function with \(\dim(g) \geq \dim(\theta)\). This setting includes overidentified models, where \(\dim(g) > \dim(\theta)\), and is typically addressed using the generalized method of moments (GMM) introduced by GMM. Adapting SGD to such moment-based estimation problems poses significant challenges because the standard notion of a stochastic gradient does not directly apply. This reflects a fundamental methodological distinction between modern scalable optimization techniques and classical GMM estimation frameworks in econometrics.
To illustrate these differences, consider the standard GMM setup. Let \(\theta_o \in \Theta \subset \mathbb{R}^d\) denote the true parameter vector, where \(\Theta\) is the parameter space. Let \(z_{1:n} := (z_i)_{i=1}^n\) be a sample of size \(n\), and let \(g(z, \theta)\) be a vector of moment functions such that \(\theta \mapsto g(z, \theta)\) is differentiable for each \(z\). Define the sample averages of the moment function and its Jacobian as:
The standard GMM estimator is then given by \[ \hat{\theta}_{n,W_n} := \operatornamewithlimits{arg\hspace{0.1em} min}_{\theta \in \Theta} \bar{g}_n(\theta)' W_n \bar{g}_n(\theta), \] where \(W_n\) is a possibly data-dependent weighting matrix. When \(g(z, \theta)\) is nonlinear in \(\theta\), computing \(\hat{\theta}_{n,W_n}\) typically requires iterative numerical methods.
Establishing convergence guarantees for algorithms that compute nonlinear GMM estimators has remained a longstanding challenge, particularly in settings where the objective function is nonconvex or a consistent initial estimate is difficult to obtain without exhaustive search. robinson1988 analyzes the use of grid search to construct a consistent initial estimator. Given such an estimator \( \hat{\theta}_0 \), applying \( t \) iterations of a Newton–Raphson (NR) procedure yields a refined estimator \( \hat{\theta}_t \) that is consistent and asymptotically normal, provided the initial estimator is sufficiently accurate. In Robinson’s framework, this initial estimator is obtained through a global search over a regular grid, which becomes computationally infeasible in large-scale problems. Building on this idea, andrews1997 proposes a stopping rule to obtain a consistent first-step estimator, followed by a fixed number of NR updates. While this method avoids pre-specifying grid resolution, it still requires a global search to verify the stopping criterion, which limits its scalability. More recently, Forneron:Zhong:23 study optimization methods for nonlinear GMM without assuming convexity or requiring a consistent initial estimator, focusing on iterative procedures such as gradient descent and Gauss–Newton.
All of the above methods rely on full-sample evaluations of the moment function, its derivative, and related quantities at every update. However, in large-scale settings with very large \( n \), computing the full-sample estimator \( \hat{\theta}_{n,W_n} \) becomes computationally burdensome or even infeasible because it requires evaluating \( \bar{g}_n(\theta) \) and \( \bar{G}_n(\theta) \) over the entire dataset in each step. Moreover, the usual form of SGD in (ref) does not directly apply, since the gradient of the GMM objective function involves a product of two sample averages: \( 2 \bar{G}_n(\theta)' W_n \bar{g}_n(\theta) \). Consequently, neither a single observation nor a single mini-batch yields an unbiased estimate of this gradient.
To address these challenges, we propose a stochastic approximation method that updates the parameter vector using independent mini-batches of the moment function and its Jacobian. The construction is analogous to a U-statistic, yielding unbiased estimates of updating directions that form a martingale difference sequence. This structure permits the use of martingale convergence theory to establish almost sure convergence, while providing substantial scalability advantages for large-scale nonlinear GMM problems.
There are a couple of related papers in the literature. SGMM:2023 extend stochastic approximation to GMM, referred to as SGMM, and develop methods for online inference in the context of linear instrumental variable (IV) regression, where $g(z_i, \theta)$ is linear in $\theta$. SGMM employs a second-order optimization method, in the terminology of bottou2018optimization, and achieves asymptotic efficiency comparable to two-stage least squares (2SLS) and GMM. However, SGMM faces notable limitations when applied to nonlinear GMM: it requires a consistent initial estimator, and its second-order nature limits scalability relative to first-order methods. As a related approach, OGMM:2025 develop an online generalized method of moments (OGMM), designed for stationary and ergodic time series data that arrive in batches. They assume that the initial batch size grows to infinity and that a consistent GMM estimator can be constructed from this batch. Moment restrictions are subsequently updated using a first-order Taylor approximation. Although their method accommodates time series dependence, it still relies on the availability of a consistent initial estimator from a large batch.\footnote{The literature on stochastic approximation for overidentified moment restrictions is sparse. Beyond SGMM:2023 and OGMM:2025, there are a few papers on online linear IV regressions in computer science and engineering communities: see, for example, venkatraman2016online, Chen:Roy:Hu:Balasubramanian:2024, della2023online and the references therein. None of these works consider the asymptotic distributions of their online linear IV estimators and statistical inference procedures, however.}
A key unpleasant assumption in both SGMM and OGMM for nonlinear moment restrictions is the requirement of an initial consistent estimator for $\theta_o$. To address this issue, we develop SLIM (Stochastic Learning and Inference in overidentified Models), which is a modular stochastic approximation framework for nonlinear GMM that includes a required first-order algorithm and an optional (but recommended) second-order algorithm that can be applied subsequently to improve efficiency. Notably, our procedure does not require an initial consistent estimator. Inference can be conducted using either a random scaling approach, which is preferred when loading the full dataset is computationally expensive since the random scaling term is updated within either the first-order or second-order algorithm, or a plug-in approach, which is preferred when full-sample access is feasible.
Furthermore, we address the challenge of instability in single-pass stochastic approximation, where each observation is accessed only once, by formally establishing asymptotic theory for multi-pass stochastic approximation with mini-batches. This greatly enhances both the stability and accuracy of the method. While SGMM:2023 demonstrate the empirical benefits of multi-pass strategies in linear IV regression, they do not provide formal theoretical guarantees, which we develop in this paper.
The organization of the paper is as follows. Section (ref) introduces the setting and presents our proposed SLIM algorithm. It provides a high-level roadmap that outlines the procedure and offers a heuristic explanation of how and why stochastic approximation works in the context of GMM. It also includes a discussion of related methods in the literature. Section (ref) develops asymptotic theory under a random-sampling framework that incorporates both stochastic approximation and sampling uncertainty. Building on this theory, Section (ref) proposes random scaling methods for statistical inference. Section (ref) extends our approach to efficient estimation via second-order stochastic approximation and introduces plug-in inference methods. We compare random scaling and plug-in inference methods, showing that both are asymptotically valid but suited to different computational regimes. In addition, we extend the Sargan-Hansen $J$-test to plug-in, debiased plug-in, and online versions tailored to stochastic learning. Section (ref) presents the results of Monte Carlo experiments that illustrate our methodology using the Exact Affine Stone Index (EASI) demand model developed by EASI. In the EASI system, budget shares depend on implicit utility, which is a nonlinear function of unknown model parameters and is endogenous due to its dependence on the budget shares via the Stone index. We estimate these parameters using nonlinear GMM, designing a data-generating process based on the dataset constructed by EASI and scaling up the sample size in our experiments. To evaluate both computational and inferential performance, we compare standard full-sample nonlinear GMM with our proposed methods. When the sample size is 100{,}000, the number of moment conditions is 576, and the number of unknown parameters is 380, Stata’s conventional GMM routine succeeds in computing the estimator and its confidence set only after 18 hours on a high-performance laptop. In contrast, our efficient stochastic approximation methods complete the same task in under 1.4 hours, demonstrating substantial gains in scalability and practical feasibility for large-scale applications. We also demonstrate that the debiased plug-in $J$-test delivers satisfactory finite-sample inference, and SLIM scales efficiently as the sample size increases to one million. Finally, Section (ref) provides concluding remarks and outlines directions for future research. We present the proof of consistency immediately following the conclusion, as establishing the consistency of SLIM is of fundamental importance given that the algorithm does not require initialization with a consistent preliminary estimator. Appendix (ref) develops an alternative mode of asymptotic theory and inference under a fixed-sample framework that captures the error due to stochastic approximation. Appendix (ref) introduces a warm-start algorithm that can be used as an optional initialization step for SLIM. Appendix (ref) offers additional discussion that complements Section (ref). Appendix (ref) contains all proofs omitted from the main text.
We begin by providing a roadmap for stochastic approximation in nonlinear GMM. Figure (ref) outlines the recommended procedure. The roadmap emphasizes that an initial consistent estimator is not required, regardless of whether a warm start (described in Appendix (ref)) is used. The first main step, which is required, applies a first-order approach based on U-statistics (see the current section for details). The second step, which is optional but recommended, refines this estimator using a second-order update to improve efficiency (see Section (ref)). The procedure branches into two alternative inference methods: the random scaling approach, described in Section (ref), preferred when loading the full dataset is computationally expensive; and the plug-in approach, described in Section (ref), preferred when full-sample access is feasible.
In our setting, computational efficiency is as important as statistical efficiency. Thus, for the moment, we assume that the weighting matrix is fixed and positive-definite but not necessarily optimal. To reduce computational costs, we use stochastic approximation, where $\theta$ is iteratively updated using randomly selected subsamples of the data. Specifically, at each iteration, two independent mini-batches of observations are used to compute unbiased estimates of $\bar{g}_n(\theta)$ and $\bar{G}_n(\theta)$. We show that this approach significantly reduces the computational burden while remaining consistent and asymptotically normal under mild regularity conditions, even though it may sacrifice statistical efficiency compared to the fully efficient GMM estimator. We return to the issue of efficiency and optimal weighting in Section (ref).
Since the weighting matrix $W_n$ remains fixed throughout the iterations, we further simplify the optimization objective by setting $W_n$ to the identity matrix, without loss of generality.\footnote{It is implicitly assumed that it is computationally straightforward to compute $W_n^{1/2}$ and pre-multiply it to $g(z_i, \theta)$, so that we can re-define $g(z_i, \theta)$ to be $W_n^{1/2}g(z_i, \theta)$.} This leads to the following simplified formulation of the estimator:
The following subsection details our stochastic approximation.
We begin by specifying the input for our stochastic approximation algorithm. Let $\theta_0^* \in \Theta$ denote the initial value of the parameter vector, which is computed based on the information in $\mathcal{F}_n:= \sigma(z_{1:n})$ and remains bounded as $n$ increases. For example, $\theta_0^*$ can be a vector of predetermined constants (e.g., zeros) or a computationally feasible estimator derived from a random subsample of $z_{1:n}$ when $n$ is very large. A key advantage of our approach is that consistency and asymptotic normality of the final estimator do not depend on the initial estimator being consistent, although the quality of the starting value can have a non-negligible impact on finite-sample performance.
In addition, the algorithm takes as input the total number of iterations $N$, the learning rate schedule $\{\gamma_t : t=1, \ldots, N\}$, and the mini-batch sizes $B_G$ and $B_g$, where $B = B_G + B_g$ denotes the total mini-batch size. These quantities are predetermined before the algorithm begins. Let $\{\tilde{z}_j : j=1, \ldots, B N\}$ denote a sequence of i.i.d. draws sampled uniformly at random from $z_{1:n}$. That is, for $j = 1$ to $B N$, sample $i(j)$ uniformly from $\{1, \ldots, n\}$ and let $\tilde{z}_j = z_{i(j)}$. This sequence provides the mini-batches of data used for the stochastic updates. When $n$ is relatively large, we set $B N = n$; otherwise, we allow $B N \geq n$, which is often referred to as multi-pass algorithms. When $n$ is relatively small, multi-pass algorithms can improve the finite-sample performance as they allow the algorithm to iterate over the data multiple times.
At each iteration $t \geq 1$, the parameter vector is updated using a mini-batch gradient descent step, which can be expressed as:
where
In this update, conditional on $\theta_{t-1}^*$, the derivative matrix $\tilde{G}_{t} (\theta_{t-1}^*)$ serves as an unbiased estimate of the full-sample derivative matrix $\bar{G}_n(\theta_{t-1}^*)$, while the moment vector $\tilde{g}_{t} (\theta_{t-1}^*)$ is an unbiased estimate of the full-sample moment vector $\bar{g}_n(\theta_{t-1}^*)$. By construction, $\tilde{G}_{t} (\theta_{t-1}^*)$ and $\tilde{g}_{t} (\theta_{t-1}^*)$ are independent of each other conditional on $\theta_{t-1}^*$. On one hand, the use of independent mini-batches for these two components ensures that the stochastic updates remain unbiased conditional on $\theta_{t-1}^*$, while significantly reducing the computational burden relative to full-sample evaluations. On the other hand, mini-batches also provide variance reduction compared to stochastic approximation with $B_G = B_g = 1$. This variance reduction is crucial in finite samples, as each iteration involves a product of two noisy elements, despite both being unbiased and independent of each other.
At each iteration, the Polyak1990-ruppert1988efficient average is computed: \[ \bar{\theta}^*_t = \frac{1}{t} \theta^*_t + \frac{t-1}{t} \bar{\theta}^*_{t-1} \] and our proposed estimator is $\bar{\theta}^*_{N}$. The averaging step helps mitigate noise in the stochastic updates. Our proposed estimator is consistent and asymptotically normal under mild regularity conditions as both $n$ and $N$ tend to infinity. Furthermore, we show how to conduct inference using the stochastic path $\{\theta_t^*: t=1,\ldots, N \}$. In summary, the proposed first-order method is outlined in Algorithm (ref).
We begin by stating the basic assumptions that underpin our analysis.
The sequence $\{ \tilde{z}_j \}_{j=1}^{BN}$ is generated by uniform sampling with replacement from $z_{1:n}$, and is partitioned into $N$ mini-batches of size $B$. When $BN = c \cdot n$ for some $c > 1$, our algorithm, as described in Section (ref), samples each data point $c$ times on average. This framework is commonly referred to as a multi-pass algorithm. An alternative sampling strategy is discussed in Appendix (ref).
The learning-rate sequence $(\gamma_t)_{t=1}^N$ is assumed to decay at the rate $O(t^{-a})$ for some exponent $a\in(1/2,1)$. This ensures $\sum_{t=1}^\infty \gamma_t = \infty$ and $\sum_{t=1}^\infty \gamma_t^2 < \infty$, which are standard conditions guaranteeing convergence of the stochastic gradient-descent algorithm. We also assume the moment function $g(z,\theta)$ is twice differentiable in $\theta$, as is customary in the nonlinear GMM setting.
We develop asymptotic theory under the i.i.d. sampling. Extensions of this framework to settings with time series or cluster dependence are left for future research.\footnote{In Appendix (ref), we consider an alternative asymptotic framework in which the data are treated as fixed, and the objective is to approximate the full-sample GMM estimator $\hat{\theta}_n$, defined in (ref). This setting is particularly relevant when evaluating $\hat{\theta}_n$ is computationally intensive, such as when the sample size $n$ is very large.} To state the further regularity conditions, we introduce additional notation. Define the sample Lyapunov function
where $\bar{g}_n(\theta)$ and $\hat{\theta}_{n}$ are given in (ref) and (ref), respectively. The population moment function and its Jacobian are also denoted by
Let $I_{d}$ denote the $d$-dimensional identity matrix. For square matrices $A$ and $B$, we write $A \le B$, $A\ge 0$, and $A >0$ if $B - A$ is positive semidefinite, $A$ is positive semidefinite, and positive definite, respectively. We use $\|\cdot\|$ to denote the Euclidean norm on vectors and the spectral norm on matrices.
The bulk of Assumption (ref) reflects standard regularity requirements under a random-sampling framework. Condition (ref) assumes that the data are i.i.d. Condition (ref) imposes global point identification of the true parameter $\theta_o$, while (ref) ensures that $\theta_o$ is locally well-separated via a rank condition on the Jacobian $G(\theta_o)$ and the covariance matrix $\Omega_o$. Condition (ref) imposes integrability on $g(\cdot,\theta_o)$ and the envelope function for the Jacobian $G(\cdot, \theta)$. Condition (ref) requires Lipschitz continuity of $\theta \mapsto G(z, \theta)$, with an integrable Lipschitz constant $L(z)$. Condition (ref) assumes the existence of a constant $M$ such that the difference between $M I_{d}$ and the Hessian of the sample Lyapunov function $\bar{Q}_n(\theta)$ is uniformly positive semidefinite, and that the ratio of $n^{-1} \sum_{i=1}^n \|g(z_i, \theta)\|^{2p}$ to $\bar{Q}_n(\theta)^p + 1$ is uniformly bounded by $M$, both with probability approaching one. This condition is used to control the variability of the stochastic errors.
It remains to discuss Conditions (ref) and (ref), which impose requirements on the behavior of the sample Lyapunov function in the regions where $\|\theta - \hat{\theta}_n\| \le \delta$ and $\|\theta - \hat{\theta}_n\| \ge \delta$, respectively. The local lower bound in (ref) ensures that $\bar{Q}_n(\theta)$ is strongly convex in a neighborhood of $\hat{\theta}_n$. In turn, (ref) requires that both $\bar{Q}_n(\theta)$ and the squared norm of the gradient of the sample Lyapunov function, $\|\bar{G}_n(\theta)' \bar{g}_n(\theta)\|^2$, are bounded away from zero for all $\theta$ such that $\|\theta - \hat{\theta}_n\| \ge \delta$. Notably, we do not impose global convexity or compactness of the parameter space; instead, convergence is ensured by the two-region structure encoded in (ref) and (ref). These two conditions jointly serve a role analogous to global strong convexity in ensuring convergence of stochastic approximation algorithms. Condition (ref) provides local curvature near the solution, ensuring that stochastic iterates are attracted to the unique minimizer $\hat{\theta}_n$ eventually. Meanwhile, Condition (ref) prevents the algorithm from stagnating far from the solution by requiring that the gradient magnitude remains bounded away from zero outside the neighborhood. This two-region structure, consisting of local strong convexity near the optimum and a sharp separation condition away from it, allows us to establish stability and convergence without relying on global convexity or compactness. Such an approach is particularly well suited to overidentified nonlinear GMM settings.
To analyze the stochastic properties of our algorithm, it is useful to decompose $(\bar{\theta}^*_{N} - \theta_o)$ into two parts:
where the first term captures the discrepancy arising from stochastic approximation relative to the full-sample GMM estimator, and the second reflects the estimation error of the full-sample GMM estimator relative to the true parameter vector. In econometrics, the conventional analysis focuses only on the statistical error, implicitly assuming that the computational error is asymptotically negligible. Intuitively, this decomposition separates the error from using a stochastic algorithm (computational) from the inherent sampling uncertainty in estimation (statistical).
We now introduce some notation. Let $\mathbb{P}^*_n$ denote the distribution induced by the mini-batch sample $(\tilde{z}_j)_{j \ge 1}$, and let $\mathbb{E}^*_n$ denote the corresponding expectation. With this in place, we are ready to state our first theoretical result.
In other words, the probability of the event $\{\|\theta_N^* - \hat{\theta}_n\| \ge \varepsilon\}$ converges to zero uniformly in $n \in \mathbb{N}$ as the number of iterations increases, indicating that $\|\theta_N^* - \hat{\theta}_n\|$ decays at a rate depending on $N$ but not on the sample size $n$. Since $\mathbb{E}^*_n \bar{Q}_n(\theta_N^*)$ decays at rate $O(\gamma_N) = O(N^{-a})$ with probability approaching $1$, the convergence rate of $\theta_N^*$ toward $\hat{\theta}_n$ is $O(N^{-a/2})$; see Lemma (ref) in the appendix.
The standard parametric rate $O(N^{-1/2})$ can be achieved by averaging over the stochastic path, yielding the averaged estimator $\bar \theta_N^* = N^{-1} \sum_{t=1}^N \theta_t^*$ with reduced variance. This estimator also satisfies the (martingale) central limit theorem (CLT), as established in Theorem (ref). Its asymptotic variance closely matches that of the full-sample GMM estimator, with an additional term arising from the uncertainty in the estimation of $G(\theta)$. To describe it precisely, let us introduce the following matrix:
which represents $\operatornamewithlimits{plim}_{n \to \infty} \bar \Sigma_n$, the asymptotic variance of the stochastic error, where $\bar \Sigma_n$ is a sample analog of $\Sigma_o$ and formally defined in the appendix (see (ref)). The second term reflects the variability in the mini-batch approximation of $G(\theta)$, which decays at a rate inversely proportional to the corresponding mini-batch size. The next theorem presents the CLT applied to the averaged estimator.
Theorem (ref) establishes that \(\sqrt{n}(\hat{\theta}_n - \theta_o)\) and \(\sqrt{N} (\bar{\theta}_N^* - \hat{\theta}_n)\) are jointly asymptotically normal and asymptotically independent. While the convergence result for the full-sample estimator \(\hat{\theta}_n\) is standard, the result for the averaged stochastic estimator \(\bar{\theta}_N^*\) is new. It shows that \(\bar{\theta}_N^*\), centered at \(\hat{\theta}_n\), is asymptotically normal, achieves the \(O(N^{-1/2})\) convergence rate, and has an asymptotic covariance matrix whose “meat” component includes an additional term
which captures the noise introduced by approximating the Jacobian with finite mini-batches. This variance inflation diminishes as the mini-batch size \(B_G\) increases.
By leveraging global information from $(\bar \theta_t^*)_{t=1}^N$, it is possible to extract additional information about the asymptotic variance of the stochastic approximation estimator. This idea is formalized in the following theorem, which establishes the functional central limit theorem (FCLT) for the partial sum process of stochastic updates, generalizing Theorem (ref) under a slightly stronger moment condition.
Throughout, let $Z$ denote a standard normal vector and $W(r)$ a standard multivariate Wiener process on $[0,1]$, independent of $Z$. The Brownian Bridge is denoted by $\bar{W}(r) = W(r) -r W(1)$. When necessary to evade confusion, we specify their dimensions by adding a subscript, e.g., $Z_d$. Let $\rightsquigarrow$ denote the weak convergence in $\ell^\infty([0,1])$. Moreover, let \(\chi^2_\ell\) denote a chi-squared variable with \(\ell\) degrees of freedom.
Theorem (ref) strengthens Theorem (ref) by establishing weak convergence of the entire stochastic trajectory rather than only its terminal average. The limiting process is a multidimensional Wiener process, scaled by the same asymptotic variance as in the CLT. This result enables the development of inference procedures that utilize the full stochastic path without recomputing averages of the sample moment function or Jacobian over the entire dataset. As a result, inference can be conducted efficiently even in large-scale settings where loading the full dataset is computationally infeasible. In summary, Theorem (ref) provides the theoretical foundation for the random scaling inference methods introduced in the next section.
We conclude this section by commenting on the restrictions imposed on \( p \), \( a \), and \( N \) by our asymptotic framework. To illustrate, consider the case \( a = 0.501 \), which is the value used in our Monte Carlo experiments. The moment condition \( p > (1 - a)^{-1} \) then requires \( p > 2.004 \), which is relatively mild. This reflects the fact that a slower decay of the learning rate (i.e., smaller \( a \)) permits weaker moment conditions. At the same time, the requirement \( N^{1 - a} / n \to 0 \) imposes a restriction on how fast the number of stochastic updates \( N \) can grow relative to the sample size \( n \). For \( a = 0.501 \), this implies \( N = o(n^{2.004}) \), meaning that \( N \) cannot grow faster than approximately \(n^2\). This example highlights the trade-off in choosing \( a \): slower decay relaxes the required moment conditions but tightens the allowable growth rate of the number of iterations.
In this section, we develop methods for conducting inference on the true parameter vector $\theta_o$. We focus on testing $\ell \leq d$ linear restrictions of the form \[ H_0 : R \theta_o = c, \] where $R$ is a known $(\ell \times d)$ matrix of rank $\ell$, and $c$ is a known vector in $\mathbb{R}^\ell$.
We consider two approaches to inference: one based on random scaling and the other on plug-in estimation of the asymptotic variance matrix. Random scaling inference has been used in two largely disconnected branches of the literature. In time series econometrics, it has long been employed to estimate the long-run variance, which is the sum of all autocovariances, since its introduction by kiefer2000simple. More recently, lee2021fast demonstrated its usefulness for inference with SGD, after which it has been adopted in stochastic approximation methods across various estimation problems. For instance, lee2025fast applied it to quantile regression models, while SGMM:2023 used it for linear GMM estimation. Additional applications of random scaling can be found in li2021statistical,ChenLai2021,Du:Zhu:Wu:Na:2025, among others. The plug-in method is deferred to Section (ref), where it is naturally integrated with the efficient estimation procedure developed there.
Define the $(\ell \times \ell)$ random scaling matrix by
This matrix admits the following iterative updating rule for each $t \geq 1$:
with initial values $A_0 = 0_{\ell \times \ell}$ and $b_0 = 0_\ell$. The Wald statistic based on $V_N(R)$ is asymptotically pivotal, as established by the FCLT in Theorem (ref).
\noindentComputational advantage: The iterative structure in (ref)--(ref) enables inference based on random scaling to be implemented efficiently. Crucially, it avoids computing or storing full $d \times d$ matrices and instead requires updates and storage only at the reduced $\ell \times \ell$ level. This leads to substantial computational gains, particularly when the number of restrictions $\ell$ is much smaller than the parameter dimension $d$. For instance, in many applications we are interested in a single linear combination of the parameter vector, so that $\ell = 1$, making the gain especially pronounced.
Since we can decompose $(\bar \theta^*_{N} - \theta_o)$ as the sum of $(\bar \theta^*_{N} - \hat{\theta}_n)$ and $(\hat{\theta}_n - \theta_o)$, the relative rates at which $n$ and $N$ diverge affects asymptotics. We consider four regimes: two intermediate cases and two polar extremes. The leading case is discussed below, and three other cases are in Appendix (ref).
We begin with the scenario where $N$ and $n$ tend to infinity at the same rate. Define the Wald statistic as
Here, the scaling factor $(n^{-1} + (NB_g)^{-1})^{-1}$ accounts for the combined uncertainty from both stochastic approximation and random sampling. From the decomposition \[ R \bar\theta^*_N - c = R(\bar\theta^*_N - \hat \theta_n) + R(\hat \theta_n - \theta_0), \] under $H_0 : R \theta_o = c$, and Theorem (ref), the asymptotic variance of $R \bar\theta^*_N - c$ is given by
because $B_g \Sigma_o \ge G_o' \Omega_o G_o$. Hence, we obtain the following corollary as a direct implication of Theorem (ref).
In Corollary (ref), we have the following distributional equivalence:
which follows from the fact that \( V_N(R) \) is asymptotically a functional of the process $\bar{W}(r)= W(r) - r W(1)$, up to an unknown variance matrix. Since \( W(1) \) is projected out, \( \operatorname{Cov}(\bar{W}(r) , W(1)) = 0 \) and thus the process $\bar{W}(r)$ is independent of \( W(1) \) due to the Gaussianity .
Recall that the test statistic $T_{n,N}$ defined in (ref) takes the form
The analysis across all asymptotic regimes considered in Section (ref) and Appendix (ref) shows that asymptotically valid inference is possible under any divergence rates of $n$ and $N$, by employing the following decision rule:
where $\mathrm{cv}$ is a critical value that depends only on the rank of $R$ and the nominal significance level. This general validity stems from the scaling factor \[ \tau_n := \left(n^{-1} + (N B_g)^{-1} \right)^{-1}, \] which adapts automatically to different asymptotic regimes. When $N/n \to 0$, we have $\tau_n / (N B_g) \to 1$; when $n/N \to 0$, we have $\tau_n / n \to 1$. As a result, the procedure accommodates both sources of uncertainty—stochastic approximation and random sampling—without requiring explicit knowledge of their relative magnitudes. If inference precision is a priority, it is advisable to increase both $N B_g$ and $B_G$ as much as the available computational budget allows.
We conclude this section with a summary of the random scaling inference procedure and practical guidance for its implementation. The test statistic $T_{n,N}$ enables valid inference across a wide range of asymptotic regimes by adapting to the relative growth of $n$ and $N$. Its implementation is computationally efficient due to recursive formulas for $V_N(R)$, which require only $\ell \times \ell$ operations and avoid computing full $d \times d$ matrices. The method also bypasses the need to estimate the asymptotic variance explicitly, making it particularly attractive for large-scale problems.
For practical use, we recommend choosing mini-batch sizes $B_G$ and $B_g$ to balance statistical precision and computational cost. Increasing $B_G$ improves efficiency but raises computational burden. Run Algorithm (ref) for a sufficiently large number of iterations $N$, typically with $N \geq C \cdot n$ for some constant $C > 1$, unless $n$ is extremely large. Inference can be conducted by adding the random scaling inference component to Algorithm (ref). This procedure provides a unified and robust method for inference across different asymptotic regimes. The limiting distribution in (ref) is pivotal and depends only on the number of restrictions $\ell$. When $\ell = 1$, critical values $\mathrm{cv}_{1 - \alpha/2}$ are given in Abadir:Paruolo:97. For instance, $\mathrm{cv}_{1 - \alpha/2} = 6.747$ when $\alpha = 0.05$, and $\mathrm{cv}_{1 - \alpha/2} = 5.323$ when $\alpha = 0.10$. The corresponding $(1 - \alpha)$ confidence interval for $R \theta_o$ is given by \[ \left[ R \bar{\theta}_N^* \pm \mathrm{cv}_{1 - \alpha/2} \cdot \sqrt{n^{-1} + (N B_g)^{-1}} \cdot (B_g V_N(R))^{1/2} \right]. \] For general $\ell > 1$, critical values are available from kiefer2000simple.
Although the stochastic approximation steps in (ref) yield asymptotically normal estimators, they are generally not efficient. The potential gains from improving efficiency can be substantial, especially when the sample size \(n\) is not exceedingly large. To address this limitation, we propose optional (but recommended) refinement steps that incorporate second-order information.
Let \(\bar{\theta}^*_N\) denote the averaged estimator after the \(N\)-th iteration. Compute:
where $W_{\mathrm{MB}}$ is based on $M_{\mathrm{MB}}$ randomly sampled mini-batches, \({\dagger}\) denotes the generalized inverse, and $M_{\mathrm{MB}}$ is a tuning parameter that diverges as \(n \to \infty\).
Let \(B_{G,t}\) be an increasing and diverging sequence after the \(N\)-th iteration; for example, \(B_{G,t} = B_{G,0} + \log (t-N) \cdot 1( t > N)\), where $B_{G,0}$ is the mini-batch size used in the first step. Define \(B_t := \sum_{s=1}^{t} (B_{G,s} + B_g)\) as the total number of mini-batch draws up to iteration \(t\). Choose a predetermined value \(T > 0\) (e.g., \(T = C \cdot n\) for some constant \(C > 0\)), and sample \(\{\tilde{z}_j\}_{j=B(N+1)}^{B_T}\) uniformly at random from the original dataset \(\{z_1,\ldots,z_n\}\). For \(t = N+1\) to \(t = T\), sequentially update:
where \(\tilde{G}_{t} (\theta)\) and \(\tilde{g}_{t} (\theta)\) are re-defined as
replacing the definitions in (ref). This procedure replaces lines 7 and 8 in Algorithm (ref) with (ref) and (ref), respectively.
It is important to emphasize that \(\Phi_n\) and \(W_{\mathrm{MB}}\) are computed only once and reused throughout the refinement step. Since \(\bar{\theta}_N^*\) is already a consistent estimator, there is no need to recompute \(\left(\Phi_n' W_{\mathrm{MB}} \Phi_n\right)^{\dagger}\) or \(W_{\mathrm{MB}}\). However, updating \(\tilde{G}_{t}(\theta)\) and \(\tilde{g}_{t}(\theta)\) at each iteration remains essential, as it leverages the Lyapunov-function-based convergence guarantees of the algorithm.
\paragraph{Choice of \(\Phi_n\) and \(W_{\mathrm{MB}}\)} Alternative choices for \(\Phi_n\) and \(W_{\mathrm{MB}}\) are possible, as long as they satisfy the consistency conditions: \[ (\Phi_n' W_{\mathrm{MB}} \Phi_n)^{-1} \overset{\mathbb{P}}{\to} (G_o' \Omega_o^{-1} G_o)^{-1}, \quad \text{and} \quad W_{\mathrm{MB}} \overset{\mathbb{P}}{\to} \Omega_o^{-1}. \] When the sample size \(n\) is relatively small, a natural alternative to $W_{\mathrm{MB}}$ is the full-sample analog: \[ W_n = \left( \frac{1}{n} \sum_{i=1}^{n} g(z_i, \bar \theta_N^*) g(z_i, \bar \theta_N^*)' \right)^{\dagger}. \] By contrast, \(\Phi_n\) is based on the full sample by default, since it does not involve matrix inversion and its mini-batch counterpart tends to yield similar results due to averaging. However, when \(n\) is too large for full-sample computation, one can approximate \(\Phi_n\) using either of the following alternatives: (i) a simple average over a randomly selected subsample of the data, or (ii) an average over multiple randomly drawn mini-batches. Both approaches offer computational savings and yield similar results due to averaging.
Define the modified Lyapunov function,
and the corresponding GMM estimator,
We are now ready to state the additional assumptions required for this section.
Conditions (ref)–(ref) mirror Conditions (ref)–(ref), but are formulated for the weighted objective \(\bar{Q}_{n,W}(\theta)\) in place of \(\bar{Q}_n(\theta)\). Condition (ref) allows for a sequence \(B_{G,t}\) that either converges to a finite constant or diverges. Condition (ref) governs the joint growth of \(N\) and \(T\) relative to the sample size \(n\). To interpret this condition, suppose \(a = 0.501\), as used in our Monte Carlo experiments. The requirement \((T - N)^{1 - a}/n \to 0\) implies that \(T - N\) must grow slower than \(n^{1/(1 - a)}\), or approximately \(n^{2.004}\). Likewise, the conditions \(N^{1 - a}/n \to 0\) and \(N^a/n = O(1)\) are satisfied when \(N = O(n^{1/a})\), or approximately \(O(n^{1.996})\). Therefore, both conditions are satisfied as long as \(N\) and \(T - N\) grow at a polynomial rate in \(n\), with exponents slightly below 2, ensuring the validity of the refinement procedure in large samples.
Denote the additional variance component by \[ V_{{G}} = (G_o' \Omega_o^{-1} G_o)^{-1} \mathbb{E}\left[ \left(G(z_i, \theta_o) - G(\theta_o)\right)' \Omega_o^{-1} \left(G(z_i, \theta_o) - G(\theta_o)\right) \right] (G_o' \Omega_o^{-1} G_o)^{-1}, \] which reflects the variance of the mini-batch estimator $\tilde{G}_t(\theta_{t-1}^*)$ as $\theta_{t-1}^* \to \theta_o$.\footnote{Note that \(V_G/B_G\) differs from (ref) because we now employ the optimal weighting matrix.}
Theorem (ref) establishes the CLT and FCLT for the efficient second-order algorithm, extending the results in Theorems (ref) and (ref). The limiting distributions capture the efficiency gains achieved through the second-order update rule in (ref). In particular, the variance of the asymptotic distribution for \(\sqrt{(T - N) B_g} (\bar{\theta}^*_T - \hat{\theta}_{n,W})\) closely approximates that of the efficient GMM estimator, up to an additional variance term \(V_G / B_G\), which arises from the stochastic approximation \(\tilde{G}_t(\theta)\) to $\bar{G}_n(\theta)$. This additional variance vanishes in the limit as \(B_{G,t} \to B_G = \infty\), in which case \(\bar{\theta}_T^*\) attains the same asymptotic efficiency as \(\hat{\theta}_{n,W}\), up to a known scalar factor.
\paragraph{Random Scaling Inference}
Just as Theorem (ref) provides a foundation for random scaling inference, Theorem (ref)(ii) implies that the algorithm described in Section (ref) remains applicable with minor modifications. Specifically, we begin updating the random scaling matrix \(V_t\) at the start of the refinement stage and modify the test statistic to \[ (n^{-1} + ((T - N) B_g)^{-1})^{-1} (R \bar{\theta}_T^* - c)' (B_g V_T(R))^{-1} (R \bar{\theta}_T^* - c), \] while continuing to use the same critical values as before.
When \(R \theta_o\) is scalar (i.e., \(\ell = 1\)), the \((1 - \alpha)\)-level random scaling confidence interval for \(R \theta_o\) becomes
where \(\mathrm{cv}_{1 - \alpha/2}\) is the critical value defined in Section (ref).
\paragraph{Plug-In Inference}
We now turn to an important special case in which \(B_{G,t} \to B_G = \infty\). In this setting, the refinement procedure achieves near-efficiency and allows for plug-in inference based on consistently estimated quantities. This leads to the following corollary.
The corollary provides a Wald-type test statistic based on \(\bar{\theta}_T^*\), using a plug-in estimate of the asymptotic variance implied by Theorem (ref)(i). Since plug-in inference is most appropriate when the sample size \(n\) is relatively modest, we estimate \(G_o\) and \(\Omega_o^{-1}\) using full-sample quantities evaluated at \(\bar{\theta}_T^*\). Specifically, define \[ \Phi_n(\bar \theta_T^*) := \frac{1}{n} \sum_{i=1}^{n} G\bigl(z_i, \bar \theta_T^* \bigr), \quad W_n(\bar \theta_T^*) := \left( \frac{1}{n} \sum_{i=1}^{n} g(z_i, \bar \theta_T^*) g(z_i, \bar \theta_T^*)' \right)^{\dagger}. \] Then, a feasible version of the Wald statistic in (ref) is
and the null \(H_0 : R \theta_o = c\) is rejected if this statistic exceeds the \(\chi^2_\ell\) critical value.
As a counterpart to (ref), a \((1 - \alpha)\) plug-in confidence interval for the scalar case of \(R \theta_o\) is given by
where \(z_{1 - \alpha/2}\) denotes the \(1 - \alpha/2\) quantile of the standard normal distribution.
Both random scaling and plug-in approaches provide asymptotically valid inference for the refined GMM estimator $\bar{\theta}_T^*$, but each offers distinct practical advantages. Random scaling avoids computing full-sample quantities such as the sample moment function and Jacobian, making it especially useful when the dataset is too large to fit in memory. By contrast, the plug-in method relies on full-sample estimates of the Jacobian and the optimal GMM weighting matrix, which can deliver more accurate and tighter confidence intervals when $n$ is moderate and full-sample access is feasible. The choice between the two methods should depend on computational constraints and the precision required.
To sharpen the discussion, it is useful to distinguish three scenarios. In the first, a strictly online setting possibly involving real-time decision-making, each data batch is observed only once and cannot be revisited. This also covers cases in which privacy concerns prevent storage of the full dataset. In such situations, only a single-pass stochastic algorithm is feasible since neither a consistent initial estimator nor multi-pass processing is available. For inference, random scaling remains viable, as it permits continuous updating without revisiting past data, while plug-in methods are infeasible because they require recomputation of full-sample statistics after estimation.
The second scenario involves sequential streaming of very large datasets. Here, the data cannot be loaded into memory all at once and must be processed in parts. Estimation again requires a stochastic algorithm due to memory constraints, but multi-pass processing is now possible. For inference, both random scaling and plug-in methods can be applied, although the plug-in approach entails more demanding memory usage.
In the third scenario, the dataset fits into memory but batch algorithms that process the entire dataset in each iteration are computationally burdensome. In this case, stochastic algorithms are preferred for estimation because of their computational efficiency. Both random scaling and plug-in inference methods are feasible, and the computational difference between them may be negligible since the plug-in method involves only sample averages and matrix inversions rather than repeated iterations over the data. Our numerical experiments in Section (ref) correspond to this third setting.
In summary, random scaling is indispensable in strictly online environments, advantageous for sequential streaming, and competitive in memory-intensive settings, while the plug-in approach is attractive when full-sample access is feasible and tighter confidence intervals are desired. Table (ref) provides a compact overview.
We propose Sargan--Hansen specification tests based on the refined GMM estimator $\bar{\theta}_T^*$. The first is the plug-in $J$-statistic, \[ J := n \, \bar{g}_n(\bar{\theta}_T^*)' W_{\mathrm{MB}} \bar{g}_n(\bar{\theta}_T^*). \]
Because $\bar{\theta}_T^*$ does not solve the first-order condition $\Phi_{n}(\theta)' W_{n} \bar{g}_{n}(\theta)=0$ of the full-sample GMM criterion $\bar{g}_{n}(\theta)' W_{n} \bar{g}_{n}(\theta)$, we debias the plug-in moment condition before constructing the test. Let $\bar{\Phi}=\Phi_{n}(\bar{\theta}_T^*)$, $\bar{W}=W_{n}(\bar{\theta}_T^*)$, and $\bar{g}=\bar{g}_{n}(\bar{\theta}_T^*)$. Define \[ \bar{g}_{n}^{D}=\Big(I-\bar{W}^{1/2}\bar{\Phi} (\bar{\Phi}' \bar{W} \bar{\Phi})^{-1}\bar{\Phi}'\bar{W}^{1/2}\Big)\bar{W}^{1/2}\bar{g}, \] which reduces to $\bar{W}^{1/2}\bar{g}$ if and only if $\bar{\theta}_T^*$ coincides with the full-sample GMM estimate. The associated debiased $J$-statistic is \[ J_{\mathrm{D}} := n \, \bar{g}_{n}^{D}{}' \bar{g}_{n}^{D} = n \, \bar{g}' \!\left(\bar{W}-\bar{W}\bar{\Phi}(\bar{\Phi}' \bar{W}\bar{\Phi})^{-1}\bar{\Phi}'\bar{W}\right)\! \bar{g}. \]
Finally, we consider an online version. Starting at iteration $t=N+1$, update the averaged moment function recursively as \[ \bar{g}_t^* = \frac{t-N-1}{t-N}\,\bar{g}_{t-1}^* + \frac{1}{t-N}\,\tilde{g}_t(\theta_{t-1}^*), \qquad t=N+1,\ldots,T, \] with $\bar{g}_N^*=0$. The corresponding online $J$-statistic is \[ J^* := \Big(n^{-1} + ((T-N)B_g)^{-1}\Big)^{-1} \, \bar{g}_T^{*'} W_{\mathrm{MB}} \bar{g}_T^*. \]
The next theorem establishes the asymptotic distributions of $J$, $J_{\mathrm{D}}$, and $J^*$.
Theorem (ref) shows that the plug-in statistic converges to a mixture of two independent chi-squared distributions, with weights determined by the limit of $n / \{(T - N) B_g\}$. The debiased version eliminates this extra variability and converges to the standard $\chi^2_{d_g - d}$. It also shows that the online statistic also converges to $\chi^2_{d_g - d}$. Equivalently, $ n \, \bar{g}_T^{*'} W_{\mathrm{MB}} \bar{g}_T^* \overset{d}{\to} (1 + \tau) \chi^2_{d_g - d}.$
Theorem (ref) highlights the key differences between these statistics. Under the null, the classical full-sample Sargan--Hansen test has limiting distribution $\chi^2_{d_g - d}$. By contrast, the plug-in statistic $J$ includes an additional $\tau \chi^2_{d}$ term, reflecting the variability introduced by using $\bar{\theta}_T^*$ instead of the efficient full-sample GMM estimator $\hat{\theta}_{n,W}$. The debiased version corrects for this source of error, while the online statistic $J^*$ absorbs it through stochastic averaging and converges to a scaled chi-squared distribution.
When $\tau$ is small (e.g.\ $n \ll (T - N) B_g$), all three tests are well approximated by $\chi^2_{d_g - d}$. When $\tau$ is non-negligible, the additional variability takes the form $\tau \chi^2_{d}$ for the plug-in statistic and $\tau \chi^2_{d_g - d}$ for the online statistic, which may lead to differences in size and power. In practice, the online version is preferable when full-sample computation is difficult, while the debiased plug-in test is preferred otherwise.
In this section, we report the results of Monte Carlo experiments. First, we discuss the rule-of-thumb methods for selecting the learning rates and mini-batch sizes.
We describe hyperparameter choices for both the warm-start and main stages. In the warm-start stage, limited prior information makes learning rate selection difficult. For Algorithm (ref) in the appendix, we set $ \gamma_e = \gamma_{0,\mathrm{ws}} e^{-a}, $ with \(a = 0.501\). We explore several candidate pairs \((\gamma_{0,\mathrm{ws}}, B_{\mathrm{ws}})\), adopting a linear scaling rule where \(\gamma_{0,\mathrm{ws}}\) increases linearly with \(B_{\mathrm{ws}}\).
For the main stage, we recommend \(\gamma_t = \gamma_0 t^{-a}\) with the same exponent. To determine \(\gamma_0\), we compute \(\tilde{G}_{t}(\bar{\theta}^*_{\mathrm{ws}})\) using mini-batches of size \(B_{\mathrm{ws}}\) and the warm-start estimate \(\bar{\theta}^*_{\mathrm{ws}}\). If \(n\) is moderate, we use all \(\lfloor n / B_{\mathrm{ws}} \rfloor\) batches; otherwise, we sample a large random subset. Let \(\mathbb{T}\) denote the selected batch indices.
Next, in the main algorithm, suppose we set the mini-batch size to \(B_{\mathrm{main}} = B_g = B_G\), which may differ from \(B_{\mathrm{ws}}\). To specify \(\gamma_0\), define \[ \Psi_0 \;:=\; \text{median}\Bigl\{ \bigl\| \tilde{G}_{t}(\bar{\theta}^*_{\mathrm{ws}})' \tilde{G}_{t}(\bar{\theta}^*_{\mathrm{ws}})\bigr\|_2 : t \in \mathbb{T}\Bigr\}, \] where \(\|\cdot\|_2\) denotes the spectral norm. We then propose the following for \(\gamma_0\):
where \(s_0 > 0\) is a tuning constant. This choice controls early-stage instability and scales \(\gamma_0\) appropriately with the batch size. Since \(\gamma_0\) is one of the most critical hyperparameters, we strongly recommend tuning \(s_0\) based on performance, as shown in our Monte Carlo results.\footnote{In a related context, bengio2012practical emphasized the importance of tuning the initial learning rate, noting: “If there is only time to optimize one hyper-parameter and one uses stochastic gradient descent, then this is the hyper-parameter that is worth tuning.”}
EASI introduced the Exact Affine Stone Index (EASI) demand model. Specifically, the implicit Marshallian budget shares ($\mathbf{w} \in \mathbb{R}^J$) for $J$ goods are given by:
where $y \in \mathbb{R}$ is the implicit utility defined as:
$\bm{z} \equiv (z_1,\ldots,z_L)' \in \mathbb{R}^L$ is a vector of demographic characteristics with $z_0 = 1$, $\mathbf{p} \in \mathbb{R}^J$ denotes a vector of log prices, $x \in \mathbb{R}$ is the logarithm of nominal total expenditures, and $\bm{\varepsilon} \in \mathbb{R}^J$ is a vector of unobserved preference shocks. Unknown parameter vectors and matrices in (ref) are: $\{ \mathbf{b}_r: r=0,\ldots,5 \}$, $\mathbf{C}$, $\mathbf{D}$, $\{ \mathbf{A}_l: l=0,\ldots,L \}$, and $\mathbf{B}$.
The EASI demand system is nonlinear in its parameters because $y$ depends on $\{ \mathbf{A}_l : l = 0, \ldots, L \}$ and $\mathbf{B}$. Furthermore, $y$ can act as an endogenous right-hand-side variable since the budget shares influence $y$ through the Stone index, $\mathbf{p}^{\prime} \mathbf{w}$. Unknown parameters in the EASI demand functions can be estimated using nonlinear GMM, provided that a sufficient number of valid instruments are available.
Their empirical analysis examines annual expenditures across $J = 9$ categories: food at home, food out, rent, household operation, household furnishing, clothing, transportation operation, recreation, and personal care. The dataset incorporates price information spanning 12 years across 4 Canadian regions, resulting in 48 distinct price vectors. The model includes $L = 5$ demographic characteristics: age minus 40, binary indicators for male, non-ownership of a car, and receipt of social assistance, along with a linear time trend (calendar year minus 1986). Their dataset comprises $n=4,847$ observations of single-person households renting their homes, aged 25 to 64, with positive expenditures on rent, recreation, and transportation. For further details on the estimation sample, refer to their Section II.B and Table 1.
As is standard practice, EASI excluded the last demand equation (personal care) to satisfy the adding-up constraint for the budget shares and estimated eight share equations. Consequently, only $J-1 = 8$ log prices, expressed relative to the omitted category, need to be included. This results in a total of $[6 + 2L + (L+2)(J-1)](J-1) = 576$ unknown parameters, which will be further reduced to 380 if Slutsky symmetry is imposed by assuming the symmetry of $\{ \mathbf{A}_l: l=0,\ldots,L \}$, and $\mathbf{B}$.
For computational implementation of nonlinear GMM, we regard the Stata code developed by Pendakur:2015 as the benchmark, which utilizes Stata's GMM command. We use the same instruments as in Pendakur:2015:
where $\bm{p}$ (slightly abusing notation) represents the $J-1$ vector of log prices relative to the omitted category. The rationale for these instruments is that $x$ is exogenous but correlated with the endogenous variable $y$. As noted by Pendakur:2015, the original instruments proposed in EASI are more complex and less transparent, as they require preliminary estimation steps. The model is just-identified when symmetry is not imposed and overidentified under the symmetry restrictions. In the remainder of the paper, we focus on the overidentified case.
As in EASI, we begin by describing the estimated Engel curves. They focused on the estimated budget shares evaluated at $\mathbf{p} = \bm{0}_J$, $\bm{z} = \bm{0}_L$, and $\bm{\varepsilon} = \bm{0}_J$, where $\bm{0}_M $ denotes the zero vector of dimension $M$: $$ x \mapsto \widehat{\mathbf{w}}(x) =\sum_{r=0}^5 \widehat{\mathbf{b}}_r x^r, $$ where $\{\widehat{\mathbf{b}}_r\}_{r=0}^5$ are the GMM estimators of $\{\mathbf{b}_r\}_{r=0}^5$. As a function of $x$, it corresponds to the Engel curves for the reference individual facing the log prices of $\mathbf{p} = \bm{0}_J$. Here, the reference individual is defined as a 40-year-old, car-owning female in 1986 who did not receive social assistance (i.e., $\bm{z} = \bm{0}_L$) with $\bm{\varepsilon} = \bm{0}_J$.
Figure (ref) shows the estimated Engel curves along with 90% pointwise confidence intervals. The domain of the curves is restricted to [-0.7, 0.9], which is well within the support of $x$. The estimated Engel curves presented in Figure (ref) differ slightly from those in EASI, as the latter are based on nonlinear 3SLS and use a different set of instruments as explained in Section (ref). For instance, the slope of the house operation Engel curve becomes negative when the nominal expenditure, $x$, is large, whereas their estimate remains nearly a positive linear segment in that range. We also find that the confidence intervals based on the GMM tend to be wider than those of the 3SLS particularly when $x$ is large.
We now turn to the estimated price effects. Instead of going over all the empirical results reported in EASI, we focus on one noteworthy case: EASI remarked, for the reference individual who faces $\mathbf{p} = \bm{0}_J$ at median expenditure where $x=0$, that “a rent price increase of 10 percent would be associated with a budget share 0.63 percentage points higher when expenditure is raised to equate utility with that in the initial situation.” This remark is based on the estimate of the relevant component in $\bm{A}_0$. In other words, the own-price compensated semi-elasticity for the rent budget share is 0.063. Using the same instruments as in Pendakur:2015, we find that now the point estimate is 0.053 and its 90% confidence interval is [0.010, 0.096].
In terms of computation, it took 27 minutes to obtain all the empirical results reported above using a laptop computer with an Apple M1 Pro chip (8-core CPU, 14-core GPU), 16GB of RAM, and macOS Ventura 13.5. As it can be computed in reasonable time with $n=4,847$, the sample size used in EASI, we carry out experiments based on this empirical example, while increasing the scale from $n=20,000$ to $n=10^6$.
Because of the simultaneity of $\mathbf{w}$ and $y$, it is challenging to simulate $\mathbf{w}$ while allowing for non-zero higher-order polynomial terms of $y$. As the main goal of the experiment is to compare the standard GMM with our stochastic approximation approach in terms of inference precision and computational cost in a well-identified setting, we simplify the data generating process in the experiments as follows.
We first set the true values of parameters at the GMM estimates reported in the previous subsection. Obtain the sample variance of the regression residuals, denoted by $\widehat{\sigma}_{j}^2$, where $j=1,\ldots,J-1$. Compute the sample average of observed budget shares, denoted by $\overline{\bm{w}}_{\textrm{sample}}$. Normalize the scale of the demographic characteristics, $\bm{z}_i$, by dividing each component by its maximum absolute value.
Simulating observations, we first randomly draw $\left\{\left(x_i, \mathbf{p}_i, \bm{z}_i\right), i=1, \ldots, n\right\}$ with replacement from the original dataset, with $\bm{z}_i$ normalized as described above, where $n$ can be larger than the original sample size. Then, generate $y_i$ by replacing $\mathbf{w}$ in (ref) with $\overline{\bm{w}}_{\textrm{sample}}$:
where $\widehat{\mathbf{A}}_l$ and $\widehat{\mathbf{B}}$ are the GMM estimates obtained in the preliminary step. Finally, for each $i=1,\ldots,n$, generate $\mathbf{w}_i$ by
where the $j$th component of $\bm{\varepsilon}_i$ is randomly drawn from $N\left(0, \hat{\sigma}_{j}^2\right)$ and the components of $\bm{\varepsilon}_i$ are mutually independent. By design, $y$ is exogenous; however, GMM using the instruments given in (ref) is still valid. One caveat is that the simulation design does not guarantee that each element of $\mathbf{w}_i$ is bounded between 0 and 1; however, it is of little concern because the moment conditions we consider do not impose the boundedness requirement either.
We employ two different types of sample sizes: We first consider $n \in \{ $20,000, 50,000, 100,000$\}$, which reflect scales that can be encountered in contemporary demand estimation scenarios. We then use $n = $ 1,000,000, an extremely large sample size, for which full-sample GMM is practically infeasible due to computational limitations. The number of simulation replications is set to 500. Our simulations are conducted on the High-Performance Computing (HPC) cluster at the University of Connecticut. The cluster features a range of CPU architectures, including Intel Skylake, AMD EPYC 64-core, and AMD EPYC 128-core processors, with compute nodes offering between 34 and 126 cores and between 187 GB and 503 GB of RAM. Computation time is measured separately using 5 replications on a laptop computer equipped with an Apple M1 Pro chip (8-core CPU, 14-core GPU), 16GB of RAM, and macOS Ventura 13.5.
We begin by conducting simulations with $n \in \{ $20,000, 50,000, 100,000$\}$. In this setting, we estimate the model using our second-order stochastic approximation approach. For the comparative purpose, we also consider the Stata GMM by Pendakur:2015. Since simulations for our method are conducted in MATLAB, to minimize simulation discrepancies between two platforms, simulated datasets are first generated and then same datasets are used in both.
Let $\theta$ denote the vectorization of model parameters $\{ \mathbf{b}_r: r=0,\ldots,5 \}$, $\mathbf{C}$, $\mathbf{D}$, $\{ \mathbf{A}_l: l=0,\ldots,L \}$, and $\mathbf{B}$ under the Slutsky symmetry restrictions. Our procedure follows the steps outlined below. First, we apply the warm-start algorithm described in Algorithm (ref) with $\theta_{0,\mathrm{ws}}^\ast=\bm{0}_{380}$. We set $\gamma_{0,\mathrm{ws}} = 0.1$ and $\alpha=0.501$. Second, we use our stochastic approximation, based on Algorithm (ref), to obtain $\theta_N^\ast$ and $\bar{\theta}_N^\ast$. In both steps, we employ the weighting matrix, \[ W_{2SLS} = I_{J-1} \otimes \left(\frac{1}{n} \sum_{i=1}^{n} \bm{q}_i \bm{q}_i^\prime \right)^{-1}, \] assuming the independence of unobserved preference shocks across equations. The initial step size $\gamma_0$ is selected based on the rule-of-thumb method proposed in Section (ref). We use $s_0=5$ and $8$ for $n=20,000$, which yield average values of $\gamma_0 = 0.082$ and $0.051$, respectively. For $n=50,000$ and $100,000$, we use $s_0=3$ and $8$, resulting in average values of $\gamma_0 = 0.136$ and $0.051$, respectively. Finally, we apply the refinement step for efficient estimation, as proposed in Section (ref). We use the mini-batch based optimal weight in ((ref)), with $M_{\mathrm{MB}}=50,000$, under the same independence restriction.
The batch sizes are set to $B_G=B_g=512$ in both the warm-start and mini-batch stochastic approximation steps. In the refinement step, we set $B_{g}=512$ and $B_{G,t} = B_{G,0} + \log(t-N)$ with $B_{G,0} = 512$, so that $B_{G,t}$ diverges as $T$ increases, ensuring asymptotic efficiency. We begin updating the efficient average estimator at the start of the refinement step: $ \bar{\theta}_t^\ast = \frac{t-N-1}{t-N} \bar{\theta}_{t-1}^\ast + \frac{1}{t-N} \theta_t^\ast $ for $t=N+1,...,T$. For inference based on stochastic approximation, we consider the two methods described in Section (ref): random scaling (RS) inference and plug-in (PI) inference (see (ref) and (ref) for the corresponding confidence interval formulas).
We first focus on estimating the Engel curves. We use the root integrated mean squared error (RIMSE) as a precision measure, that is defined as $$ \sqrt{\sum_{j=1}^{J-1} \int_{\underline{x}}^{\overline{x}} \mathbb{E} \left( \sum_{r=0}^5 \left( \widehat{b}_{j,r} - b_{j,r}\right) x^r \right)^2 dx}, $$ where $b_{j,r}$ is the $j$th element of $\mathbf{b}_r$, $[\underline{x}, \overline{x}] = [-0.7, 0.9]$ and the integral is approximated by partitioning the domain into an interval width of 0.1.
Table (ref) reports the results. In column (1), we present simulation results from Stata (full-sample) GMM. Columns (2)-(5) report results from our SLIM method, using varying numbers of iterations and step sizes. The computational cost represents the average computation time in hours. We fix the number of iterations at ($N,T$)=(20,000, 40,000) and (20,000, 60,000) for all sample sizes. The number of epochs in the warm-start step is set to $E_{\mathrm{ws}}=$ 25, 4, and 1 for $n=$ 20,000, 50,000, and 100,000, respectively, which yield approximately 37,000 to 38,000 iterations ($=\lfloor n / B_{\mathrm{ws}}\rfloor \times ( \lfloor n / B_{\mathrm{ws}} \rfloor )-1 ) \times E_{\mathrm{ws}}$). Under this design, the computation time for our method remains within 1 to 1.5 hours, regardless of the sample size.
The table shows that our SLIM demonstrates a substantial gain in computational efficiency compared to the full-sample GMM, while maintaining comparable estimation precision. Our approach tends to yield accurate estimates within 1.05 to 1.35 hours across all sample sizes considered. In contrast, the computation time for the full-sample GMM increases substantially with sample size. For instance, when $n=$ 50,000, the full-sample GMM takes 7.84 hours, while our method yields slightly larger RIMSEs in only 1.05 to 1.33 hours. The computational advantage becomes even more pronounced when $n=$ 100,000. Our method produces accurate estimates within 1.06 to 1.35 hours, whereas the full-sample GMM exceeds 18 hours, surpassing our computational budget. These results highlight that the computational gain of our procedure becomes increasingly significant when working with large datasets.
We next examine the finite sample performance based on the price effect on the rent budget share for the reference individual with $\mathbf{p} = \bm{0}_J$ and $x=0$, which corresponds to the third diagonal element of $\bm{A}_{0}$, namely $\bm{A}_{0,33}=0.0533374$.
Table (ref) summarizes the simulation results. It first reports the root mean square error (RMSE), bias and standard deviation (SD) to evaluate estimation accuracy. We observe similar patterns to those in the Engel curves estimation in Table (ref). Our SLIM produces accurate estimates, comparable to the full-sample GMM estimates in terms of RMSE, but with significantly shorter computation time. As in the Engel curve estimation, the full-sample GMM exceeds the time limit when $n=$ 100,000, while our method delivers accurate estimates in slightly over one hour.
In this simulation, we also assess the accuracy of statistical inference. The table reports the empirical coverage probabilities at the 0.95 nominal level and the average lengths of the associated confidence intervals. For our approach, we consider the RS inference as well as the PI inference. Overall, both procedures tend to yield accurate empirical coverage rates, which are comparable to the full-sample GMM inference, particularly when $s_0=3$ and $5$. Finally, the table shows that the average confidence interval lengths of our PI method are comparable to those of the full-sample GMM.
We then examine how the SLIM estimator and its associated random scaling confidence intervals behave as the number of iterations \(T\) increases. Figure (ref) shows the average SLIM path and the corresponding 95% random scaling confidence intervals over \(T\), computed by averaging across simulation replications. We set \((s_0, E_{\mathrm{ws}}) = (5, 25), (3, 4)\), and \((3, 1)\) for \(n = 20{,}000\), \(50{,}000\), and \(100{,}000\), respectively. The confidence intervals begin to stabilize around \(T = 40{,}000\) for all sample sizes.
We next evaluate the performance of the SLIM based Sargan--Hansen overidentification tests developed in Section (ref). For this experiment, we generate instruments based on \[ x - c \, \mathbf{p}^{\prime} \mathbf{w} \] in place of $x$ in ((ref)), where $c \in \{0, 1, 2\}$ controls the strength of endogeneity. When $c = 0$, the instruments are exogenous, while larger values of $c$ correspond to stronger endogeneity. We consider the proposed debiased overidentification test, as well as the plug-in and online overidentification tests. For comparison, we also employ the standard overidentification test based on the full-sample GMM.
Table (ref) reports the empirical rejection probabilities at the 0.05 nominal level. We observe that our debiased overidentification test achieves accurate empirical size, while the plug-in and online overidentification tests rarely reject the null hypothesis when $c=0$. The power of the debiased test increases with both the degree of endogeneity and the sample size. The plug-in and online tests fail to detect the endogeneity of instruments for $c=1,2$ when \(n = 20{,}000\), \(50{,}000\), while they show good power properties when \(n = 100{,}000\) and $c=2$. The standard overidentification test performs comparable to our debiased test, but requires much longer computation time as the sample size grows. As in the previous simulation experiments, the computation time exceeds our budget when $n=$100,000, preventing us from simulating the empirical rejection probabilities.
Lastly, we examine the performance of our approach using an extremely large sample size of $n=$ 1,000,000. In the setting, the full-sample GMM is infeasible due to computational limitations. We evaluate the performance of various versions of our stochastic approximation approach. We first employ the efficient second-order SLIM, as in the previous simulation. We also consider two versions of first-order SLIM methods: one based on Algorithm (ref) and one with the refined weight (that is, omitting $(\Phi_n' W_{\mathrm{MB}} \Phi_n)^{\dagger}$ in (ref)). Given that the sample size is extremely large, we use a random sub-sample of observations to implement the warm-start algorithm and to construct the initial weight for all methods, as well as $\Phi_n$ for the second-order stochastic approximation. For the first-order stochastic approximation without the refinement step, we do not update the weighting matrix with the optimal weight; thus, it corresponds to one-step (possibly inefficient) GMM estimation. For inference, we use the RS method, which is computationally efficient.
The results are reported in Table (ref). We first examine the performance in estimating the Engel curves. As shown in the table, the first-order SLIM approach without refinement produces accurate estimates with $N=300,000$ iterations, requiring slightly less than 5 hours of computation time. Moreover, the second-order approach provides additional refinement, yielding more accurate estimates within 1.11-1.38 hours. For the price effect on the rent budget share, the first-order SLIM without refinement exhibits over-rejection with $N=$ 300,000, while more reliable inference is achieved with $N=$ 500,000. The first-order approximation with refinement yields more accurate inference within 2.09-2.38 hours. Further improvement, along with reduced computation time, is achieved using the second-order SLIM approach.
This paper has proposed SLIM, a stochastic approximation method for nonlinear generalized method of moments that does not require a consistent initial estimator. The approach exploits U-statistics for both the moment vector and its derivative matrices and incorporates multi-pass mini-batch updates to improve stability and accuracy. Monte Carlo experiments show that SLIM can solve large-scale nonlinear GMM problems with hundreds of parameters and samples of size $n=100{,}000$ or more in under two hours, whereas full-sample GMM requires an order of magnitude more computation time. These results demonstrate that SLIM provides a practical and scalable tool for econometric applications that are otherwise computationally infeasible with standard methods.
Several extensions merit further study. One promising direction is to extend SLIM to GMM with nonsmooth moment functions, following Chen:Liao:2015, which would broaden its applicability to models with quantile restrictions and other non-differentiable structures. Another is to examine the case of model misspecification, for example by embedding stochastic approximation into the framework of hansen2021inference, in order to assess the robustness of SLIM. It would also be valuable to develop principled strategies for tuning hyperparameters such as mini-batch sizes, stopping rules, and learning rates, as well as to explore applications in high-dimensional settings where the number of parameters grows with the sample size. We leave these topics to future research.
As emphasized in the introduction and earlier sections, SLIM delivers a consistent estimator without the need to initialize the algorithm with a consistent preliminary estimator. Establishing the consistency of SLIM is therefore of fundamental importance. The proof of Theorem (ref) is provided below. The auxiliary lemmas required for this proof, together with the proofs of the remaining theorems, are presented in Appendix (ref).