SGMM: Stochastic Approximation to Generalized Method of Moments

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

78,391 characters · 22 sections · 58 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

abstractWe introduce a new class of algorithms, Stochastic Generalized Method of Moments (SGMM), for estimation and inference on (overidentified) moment restriction models. Our SGMM is a novel stochastic approximation alternative to the popular GMM (offline) GMM, and offers fast and scalable implementation with the ability to handle streaming datasets in real time. We establish the almost sure convergence, and the (functional) central limit theorem for the inefficient online 2SLS and the efficient SGMM. Moreover, we propose online versions of the Durbin-Wu-Hausman and Sargan-Hansen tests that can be seamlessly integrated within the SGMM framework. Extensive Monte Carlo simulations show that as the sample size increases, the SGMM matches the standard (offline) GMM in terms of estimation accuracy and gains over computational efficiency, indicating its practical value for both large-scale and online datasets. We demonstrate the efficacy of our approach by a proof of concept using two well known empirical examples with large sample sizes.

\onehalfspacing

Introduction

Machine learning techniques have revolutionized the analysis of vast and unconventional datasets. Among them, stochastic approximation (SA) or more commonly called stochastic gradient descent (SGD) pioneered by Robbins-Monro-1951 has proven highly valuable due to its computational simplicity and scalable online implementation. In econometrics, Halbert White was a great trailblazer of SGD. For example, white1989jasa applied earlier general theory on the almost sure consistency and asymptotic normality of recursive nonlinear least squares (NLS) to parametric single-hidden layer artificial neural network (ANN) regression models with independent and identically distributed (iid) data; kuanwhite1994er developed asymptotic theory for general nonlinear models of weakly dependent processes, including applications to nonlinear regression via neural networks;\footnote{Pastorello:Patilea:Renault applied the results of kuanwhite1994er to obtain the consistency and asymptotic normality for their recursive latent backfitting procedure in a just-identified moment problem.} Chen:White:1998:JET applied stochastic approximation to bounded rationality learning; chen2002asymptotic established asymptotic theory of SGD for Hilbert space-valued mixingale, dependent error processes.

While traditionally used for computational purposes, such as optimizing objective functions bottou2018optimization, SGD has also received attention for its statistical properties. As an early path-breaking work, polyak1992acceleration obtained conditions under which an average of the SGD sequence is asymptotically normal with mean zero and an efficient variance matrix in parametric regressions. A more recent literature on SGD covers diverse topics: regularized methods for high-dimensional M-estimators agarwal2010fast; implicit SGD Toulis2017,Won:lab:2022; moment-adjusted SGD liang2017statistical; non-asymptotic results for the averaged SGD pmlr-v99-anastasiou19a,pmlr-v125-mou20a; among many others. A branch of the recent literature is concerned with online statistical inference: bootstrap fang2018online; batch-means chen2020statistical,zhu2021online; random scaling ChenLai2021,lee2021fast,li2021statistical among other possible modes of inference Chee-Kim-Toulis:23. The studies by date, however, have mainly focused on M-estimation. That is, the SGD has been mainly used for estimating a parameter of interest $\beta^*$ that is identified as the unique minimizer of a population loss function $\min_{\beta} \mathbb{E}\left[ \ell_i(\beta) \right]$, where $\ell_i(\beta)$ is a known real-valued function of the $i$-th observation and a parameter $\beta$. We therefore refer to the usual SGD-type estimators as “M-type SGD”, which takes the following basic form: \[ \beta_{i}=\beta_{i-1}-\gamma_{i}\frac{\partial}{\partial \beta} \ell_i \left(\beta_{i-1}\right),~~\text{for some~}\gamma_{i}\searrow 0~. \]

In applications in economics and finance, we often encounter a different type of estimation problems, the so-called “Z-estimation,” where the parameter of interest is identified as a unique solution to a set of moment conditions, i.e., $\mathbb{E}\left[ g_i(\beta_*) \right]=0$, where $g_i(\beta)$ is a known function of the $i$-th observation and a parameter $\beta$. Here, $\beta_*$ denotes the unique solution. These moment conditions, under just (or exact) identification $\dim(\mathbb{E}\left[ g_i(\beta) \right]) = \dim (\beta)$, would yield the “estimating equations” for the parameter of interest $\beta_*$. Most importantly, the popular (offline) generalized methods of moments (GMM) of GMM allows for overidentified moment restrictions in the sense that $\dim(g_i) > \dim (\beta)$, and that efficient estimation of $\beta_*$ and model-specification test can be carried out using the same optimally weighted GMM loss function $\min_{\beta} \left(\mathbb{E}[ g_i(\beta)]'\left[\operatorname{var}(g_i(\beta^*))\right]^{-1}\mathbb{E}[ g_i(\beta)]\right)$. Unlike M-type SGD, it is unclear how to obtain a SGD alternative to the optimally weighted GMM and to establish its statistic properties. This is especially the case for the overidentified moment restriction models.

In this paper, we develop new stochastic approximation methods for GMM, allowing for possibly overidentified moment restriction models. As a premier example, we focus on linear instrumental variable (IV) regression, where the moment restrictions are linear in the parameters of interest. Despite of being restricted to linearity, this type of models are widely applicable in economics and finance applications. We argue that aside from the more traditional IV estimators (e.g., two-stage least squares), the SA-based estimation is a natural option for IV regression because of the following reasons. First, it is fully capable of handling problems of very large datasets. Because by nature of SA, the estimation is updated one-observation-at-a-time, so it is suitable for online learning. Second, it is convenient to work with the moment conditions of econometric models, the essence of possibly overidentified Z-estimation. Hence, we view our approach as a highly scalable estimation and inference method for the moment restriction models.

We first propose a stochastic approximation to the two-stage least squares (2SLS) and analyze its stochastic properties. We provide conditions under which our SA-based estimator is first-order asymptotically equivalent to the standard (offline) 2SLS estimator. The inference problems studied in this paper based on the SA-based 2SLS estimator include: obtaining confidence interval for $\beta_*$ as well as testing for the validity of the specified instruments. For the former problem, we employ the recently developed random scaling inference of lee2021fast, which is fast, suitable for online learning, and easily adaptable for subvector inference. For the latter problem, we develop an “online” version of the Durbin-Wu-Hausman test by comparing the probability limits of the OLS and 2SLS estimators, both obtained using the SA-based methods. In both problems, because the pivotal statistics are scaled by a random matrix similar to the “fixed-b” smoothing, the asymptotic distributions are mixed normal, whose critical values have been tabulated in the literature, and are readily available for statistical inference. See, e.g., kiefer2000simple,velasco2001edgeworth,sun2008optimal,sun2013ej,Chen:Liao:Sun:2014,lazarus2018har,gupta2021robust for related papers in the time series literature.

As in the regular GMM-estimation, one of the central problems is the efficient estimation of $\beta_*$ using optimally weighted moment conditions. We show that the optimal weighting is also naturally incorporated by the overidentified Z-type SA algorithm, where we sequentially update the optimal weighting matrix along the path of the SA iteration. Despite of sequentially updating an inverse of a covariance matrix, we show that implementation is still fast because it is based on the Sherman–Morrison–Woodbury (SMW) formula; in other words, our implementation does not involve an actual high-dimensional matrix inversion. In theory, we show that the optimally weighted SA-based estimator, termed stochastic GMM (SGMM), is first-order asymptotically equivalent to the well known (offline) two-step efficient GMM estimator of the IV regression model. As by-products, we provide online plug-in optimal inference on $\beta_*$ as well as an online version of the Sargan-Hansen specification test using the efficient SGMM estimator.

The literature on stochastic approximation to 2SLS or GMM is almost non-existent. Some exceptions are venkatraman2016online and della2023online. venkatraman2016online proposed to use fitted values from the first stage to build an online algorithm; della2023online considered the just-identified IV estimator to study online regression as well as the bandit problem. The aforementioned papers carried out some sort of regret analyses but none of them focused on statistical properties of their proposed methods.

The remainder of the paper is organized as follows. Section (ref) outlines our basic algorithm and its theoretical properties. Section (ref) provides an efficient online algorithm that is first-order asymptotically equivalent to efficient GMM estimators. In Section (ref), we show that online versions of the Durbin-Wu-Hausman and Sargan-Hansen tests can be seamlessly integrated within the SGMM framework. Section (ref) reports the results of extensive Monte Carlo experiments and Section (ref) provides empirical examples based on two well known studies: angrist1991does and angrist1998children. Section (ref) discusses extensions, including an extension to Nonlinear SGMM. Section (ref) provides all the proofs of the theoretical results in the main text.

Notation. We denote the $\ell^2$-vector norm of $x$ in $\mathbb R^n$ by $\|x\|_n = (\sum_{i=1}^n |x_i|^2)^{1/2}$, and the $\ell^2$-operator norm of an $n$ by $m$ matrix by $\|{A}\|_{\operatorname{op}} := \sup_{\|x\|_m\le 1}\|A x\|_n$. We will occasionally write the $\ell^2$-vector and operator norm as $\|x\|$ or $\|A\|$ suppressing the subscripts when there is no possibility of confusion. The $\ell^2$-vector norm is equal to the operator norm when vectors are seen as $n\times 1$ matrices. Let $(S, d)$ be a complete metric space. We denote the weak convergence of $S$-valued random variables $X_n$ by $X_n \rightsquigarrow X$, where $X$ denotes the weak limit. In addition, $\overset{d}{\to}$ refers to convergence in distribution. For a real sequence $a_n$ and positive numbers $b_n$, we write $a_n = O(b_n)$ or $a_n \lesssim b_n$ if there exists a uniform constant $C>0$ such that $|a_n|\le C b_n$ holds for all $n$. For a sequence of real r.v.'s, we also denote $X_n = o_{a.s.}(b_n)$ or $X_n = O_{a.s.}(b_n)$ if $\lim_{n\to \infty}|X_n|/b_n = 0$ or $\limsup_{n\to\infty} |X_n|/b_n < \infty$ with probability 1, respectively. We use $X_n = O_P(1)$ and $X_n = o_P(1)$ to indicate that $X_n$ is a tight sequence of random variables and converges to $0$ in probability, respectively.

commentMachine learning techniques have been tremendously successful in economic research, which provides new tools for analyzing big and complex economic data. In particular, stochastic gradient descent (SGD) has been one of the very useful modern techniques for data analytics due to its simplicity and convenient implementations. While SGD has been traditionally used mainly for computational purposes (e.g., optimizing objective functions), in the recent decade, extensive studies have been conducted in the machine learning literature to analyze the actual statistical estimators that are defined by the SGD-type algorithms, and these studies showed the associated statistical properties. The literature is broad and and covers a large range of statistical models such as M-estimations and high-dimensional regularized methods, e.g., agarwal2010fast,liang2017statistical,Toulis2017,pmlr-v99-anastasiou19a. Meanwhile, recent studies also achieved statistical inference for estimators defined by SGD- type algorithms, see polyak1992acceleration,chen2002asymptotic,lee2021fast,zhu2021online,fang2018online,ChenLai2021,li2021statistical,Su:Zhu,Godichon-Baggioni,jin2022federated,li2019convergence. The studies by date, however, have been focused on M-estimations. That is, the SGD has been mainly used for estimating parameter identified as the optimum of a population objective function. We therefore refer to the usual SGD-type estimations as “M-type SGD". Economic studies also often encounter a different type of estimation problems, the so-called “Z-estimations", where the parameter of interest is identified as a unique solution to a set of moment conditions. These moment conditions, under exact identifications, would yield the “estimating equations" of the parameter of interest. Meanwhile, the popular generalized methods of moments (GMM) of GMM is an important type of Z-estimations that can address the case when the parameter of interest is over identified. While SGD, as an algorithm, has still been practically applied to solve the Z-estimation problems, it is mainly used as a implementation tool. The statistical properties of SGD algorithms for Z-estimations have been an open question by date, which is the main object of interest in this paper. This paper, to our best knowledge, provides the first study of statistical inference for Z-estimations using SGD, and we refer to this type of estimators as “Z-type SGD" to distinguish from the usual M-type methods. We show how Z-type SGD can be implemented to estimate models identified by econometric moment conditions. As a premier exploration, this paper focuses on linear instrumental variable (IV) regressions, where the moment conditions are linear in the parameter of interest. Despite of restricting to linearity, this type of models are widely applicable in econometric studies. We argue that aside from the more traditional IV estimations (e.g., two-stage-least-squares), the SGD-based estimation is a natural option for IV regressions, because: first, it is fully capable to fast handle problems of very large datasets. Because by nature of SGD iteration, the estimation is updated one-observation-at-a-time, so is implemented very fast. Secondly, it is very convenient to work with “gradients" that are essentially the first-order-conditions (FOC) of the statistical model, the essence of Z-estimations. Hence we view Z-type SGD as the most potent option for scalable inferences for the moment condition based model. We develop a novel SGD-type algorithm to estimate linear IV regressions, allowing the parameter of interest to be either exactly or over identified. We analyze the statistical rate of convergence of the Z-type SGD estimator in such models, and also develope inferential theories of the “averaged SGD" estimator, defined as the average of the iterated estimators obtained along the SGD-path. The sufficient conditions of statistical properties are analogous to the convexity of the objective functions. The inference problems studied in this paper include: obtaining confidence interval based on the Z-type SGD as well as testing for the validity of the specified instruments. For the former problem, we employ the recently developed “fixed-b SGD inference" of lee2021fast, which is fast and also works well for subvector inference. For the latter problem, we develop an “online" version of the Durbin-Wu-Hausman test by comparing the probability limits of the OLS estimatior and the IV estimator, both obtained using the SGD-based method. In both problems, because the pivotal statistics are scaled by a random variable similar to the “fixed-b" smoothing, the asymptotic distributions are mixed normal, whose critical values have been tabulated in the literature, and are readily available for inferences. As in the regular GMM-estimations, one of the central problems is the efficient estimation using optimally weighted moment conditions. We show that the optimal weighting is also naturally incorporated by the Z-type SGD algorithm, where we iteratively update the optimal weighting matrix, an inverse-covariance, along the path of the SGD iteration. Despite of iterating an “inverse covariance", we show that the implementation is still very fast because the actual iteration is based on the Sherman–Morrison–Woodbury formula, so does not involve an actual matrix inversion in the implementation. In theory, we show that the optimally weighting SGD is asymptotically equivalent to the well known efficient GMM estimation of the IV regression. The instruments considered in the paper are assume to be valid and “strong". Allowing many weak IVs is an important question that has been fruitfully studied in the past two decades. From the implementation perspective, incorporating weak IVs and even more general nonlinear models is straightforward, because the FOC of these models are readily available. And we expect that the Z-type SGD is even more appealing in this case, because it can be helpful to deal with the problem of many local minimums. For instance, being fast algorithms, the Z-type SGD gives us an avenue of attempting many different starting values. However, we expect that the theoretical studies could be technically challenging, so we leave these extensions for future research.

Stochastic Approximation for Instrumental Variable Regression

Model

Consider a linear instrumental variables regression model

align[align omitted — 88 chars of source]

where $y_i \in \mathbb{R}$ is the dependent variable, $x_i \in \mathbb{R}^{d_\beta}$ a vector of covariates and some of which are endogenous in the sense that $\mathbb{E}\left[x_i u_i\right]\neq 0$, $u_i \in \mathbb{R}$ is the regression error, $z_i \in \mathbb{R}^{d_g}$ is a vector of instrumental variables, and $\beta_* \in \mathbb{R}^{d_\beta}$ is a vector of unknown true parameters of interest. Let $g_i(\beta) \equiv z_i \left( x_i'\beta - y_i \right)$ and assume $d_g \geq d_{\beta}$. We focus on estimation of $\beta_*$ using the following linear moment restriction models:

align*[align* omitted — 54 chars of source]

Throughout the paper, we let $G_i \equiv z_i x_i'$ and $H_i \equiv -z_i y_i$; hence, $g_i(\beta) = G_i \beta + H_i$. We also let $G \equiv \mathbb E [G_i]$ and $H \equiv \mathbb E [H_i]$. That is, a letter without subscript denotes its expectation. The linear instrumental variable regression model (ref) becomes $G \beta_* + H = 0$.

S2SLS Algorithm

In this section, we propose a new stochastic approximation algorithm to estimate $\beta_*$. Let $\mathbb S \equiv \{\mathcal{D}_i = (x_i, z_i, y_i)\}_{i=1}^{n}$ be an i.i.d. sample of size $n$ drawn from a population distribution satisfying model (ref). Let $\mathbb S_0 \equiv \{\mathcal{D}_{0j}= (x_{0j}, z_{0j}, y_{0j})\}_{j=1}^{n_0}$ be an initialization random sample of size $n_0$ drawn from model (ref), with $n_0 \ll n$. Denote $\mathcal F_i = \sigma(\mathbb S_0 \cup \{\mathcal{D}_j\}_{j=1}^{i})$ for $i \ge 0$. Compute the initial estimator $\beta_0 \in \sigma(\mathbb S_0)=\mathcal F_0$ using 2SLS, GMM or any other estimation methods.\footnote{In fact, our asymptotic theory allows for any arbitrary choice of the initial estimator, including $\beta_0 = 0$. The finite sample performance depends on the quality of the initial estimator, however.} Let $\Phi_0 \equiv \frac 1 {n_0} \sum_{j=1}^ {n_0} z_{0j} x_{0j}' \in \mathcal F_0$ and $W_0 \equiv \left( \frac{1}{n_0} \sum_{j=1}^{n_0} z_{0j} z_{0j}' + \eta_0 I \right)^{-1}$ for a fixed constant $\eta_0 \geq 0$.\footnote{In many applications, the choice of $\eta_0 = 0$ will suffice, provided that $n_0$ is large enough and $z_i$ is linearly independent.} Starting from $\left(\beta_0, \Phi_0, W_0\right)$, we update the stochastic process $(\beta_i, \Phi_i, W_i)_{i=1}^\infty$ sequentially as

subequations\begin{align} \beta_i &= \beta_{i-1} - \gamma_i (\Phi_{i-1}' W_{i-1} \Phi_{i-1})^{\dagger} \Phi_{i-1}' W_{i-1} g_i(\beta_{i-1}), \\ \Phi_i &= \frac{n_0 + i - 1}{ n_0 + i} \Phi_{i-1} + \frac{1} { n_0 + i} G_i, \\ m_i&=n_0+i-1+z_i'W_{i-1}z_i, \\ W_i&= \frac{n_0+i}{n_0+i-1} W_{i-1} \left[I- m_i^{-1} z_i z_i'W_{i-1}\right], \\ \bar \beta_i &= \frac{i-1 }{ i} \bar \beta_{i-1} + \frac{1} { i}\beta_i, \end{align}

where $g_i(\beta_{i-1}) = G_i\beta_{i-1} + H_i$, $\bar \beta_0 = \beta_0$ and $\gamma_i \equiv \gamma_0 i^ {-a}$ is a learning rate with some predetermined constants $\gamma_0 > 0$ and $a \in (1/2, 1)$. Here, $A^{\dagger}$ denotes the generalized inverse of a matrix $A$. We propose to use $\bar\beta_n$, which is called the Polyak1990-ruppert1988efficient average, as an estimator of $\beta_*$.

remark[Sherman-Morrison-Woodbury matrix inversion] Note that we update the $d_g \times d_g$- weighting matrix $W_{i}$ sequentially in the S2SLS algorithm (ref). To convey the idea behind this updating rule, let $Q \equiv \mathbb{E}[ z_i z_i']=W^{-1}$, and consider an updating rule for $Q_i$: \begin{align*} Q_i &= \frac{n_0 + i - 1}{ n_0 + i} Q_{i-1} + \frac{1} { n_0 + i} z_i z_i' \end{align*} and update $\beta_i$ using the inverse of $Q_i$: \begin{align*} \beta_i &= \beta_{i-1} - \gamma_i (\Phi_{i-1}'Q_{i-1}^{-1}\Phi_{i-1})^{\dagger}\Phi_{i-1}'Q_{i-1}^{-1} g_i(\beta_{i-1}). \end{align*} The proposed algorithm in (ref) explicitly computes $Q_{i-1}^{-1}=W_{i-1}$ using Sherman-Morrison-Woodbury (SMW) formula, showing that it is unnecessary to compute the inverse of $Q_{i-1}$ each time but it suffices to update the scalar quantity $m_i$ and the matrix $W_i$ accordingly. The positive constant $\eta_0$ in the definition of $W_0$ ensures that $W_{i-1}$ is well defined for all $i$.

remark[Computation of $(\Phi_{i-1}' W_{i-1} \Phi_{i-1})^{\dagger}$] When the dimension $d_\beta$ of $\beta$ is high, it would be time-consuming to directly compute $(\Phi_{i-1}' W_{i-1} \Phi_{i-1})^{\dagger}$ in (ref). Under the identification assumption given below, $(\Phi_{i-1}' W_{i-1} \Phi_{i-1})$ is invertible with probability approaching one, and hence its inverse can be computed via SMW formula. Specifically, let $\mathcal{H}_i := (\Phi_{i}'W_{i} \Phi_{i})^{-1} \in \mathbb{R}^{d_{\beta} \times d_{\beta}}$, we have: \begin{align} \mathcal{H}_i = \frac {n_0 + i} {n_0 + i-1} \left( \mathcal{H}_{i-1} - \mathcal{H}_{i-1} \mathcal{U}_i (\mathcal{D}_i + \mathcal{U}_i'\mathcal{H}_{i-1}\mathcal{U}_i )^{-1} \mathcal{U}_i' \mathcal{H}_{i-1}\right), \end{align} where \begin{align*} \mathcal{D}_i &=\operatorname{diag}\left[-m_i, n_0 + i-1 \right] \in \mathbb{R}^{2 \times 2},\\ \mathcal{U}_i &= [x_i-\Phi_{i-1}' W_{i-1}z_i, x_i] \in \mathbb{R}^{d_{\beta} \times 2}. \end{align*} Because $\mathcal{D}_i + \mathcal{U}_i' \mathcal{H}_{i-1} \mathcal{U}_i$ is a 2 by 2 matrix, using (ref) would be computationally advantageous when $d_\beta \gg 2$.

Intuition Behind Our Algorithm

One difficulty in constructing an stochastic approximation algorithm for an IV regression is that the model (ref) is possibly an overidentified moment restrictions; therefore, there is no obvious form of stochastic gradient descent for an IV regression.

Suppose that $G = \mathbb{E} [ z_i x_i']$ has full rank $d_\beta$, which is the standard assumption for IV regression. Then $ \lambda_{\min}(G' G) > 0$, where $\lambda_{\min}(\cdot)$ is the smallest eigenvalue. We propose to build our algorithm based on a stochastic approximation of the ordinary differential equation:

align*[align* omitted — 88 chars of source]

where $W$ is a positive definite symmetric weighting matrix. Note that

align*[align* omitted — 407 chars of source]

Thus, our updating rule on average moves in the direction which reduces the difference between the current state $\beta_{i-1}$ and the true value $\beta_*$, provided that $\Phi_{i-1}$ and $W_{i-1}$ are close to $G$ and $W$ to make $\Phi_{i-1}'W_{i-1}G$ positive definite. The latter requirement is satisfied for sufficiently large $i$ due to the law of large numbers. When $n$ is large enough, $(\Phi_{i-1}' W_{i-1} \Phi_{i-1})^{\dagger} \Phi_{i-1}' W_{i-1} G$ is close to an identity matrix.

remarkIt is important to pre-multiply the scaling matrix $(\Phi_{i-1}' W_{i-1} \Phi_{i-1})^{\dagger}$ in (ref), so that the scale of observations can be automatically adjusted. From this perspective, our algorithm can be regarded as a sort of second-order method using the terminology of bottou2018optimization. They emphasize that SGD or the batch gradient method, which can be called first-order methods, are not scale invariant. Another perspective, which is more intimately tied to asymptotic theory, is that our algorithm is based on an influence function for 2SLS and GMM. In other words, our formulation of the second-order matrix is reverse-engineered to reproduce the same asymptotic variances of offline 2SLS and GMM.

remarkWe multiplied $(\Phi_{i-1}' W_{i-1} \Phi_{i-1})^{\dagger} \Phi_{i-1}' W_{i-1}$ before $g_i(\beta_{i-1})=G_i\beta_{i-1}+H_i$ in (ref), as opposed to $(\Phi_{i}' W_{i} \Phi_{i})^{\dagger} \Phi_i' W_{i}$. If the latter were multiplied instead of the former, it would have introduced $O(\gamma_i /i)$-bias in each update, whose exact magnitude depends on the data-generating process. Furthermore, asymptotic analysis would have been more complicated.

remarkIt is noteworthy that the updating rule depends only on the relative location to the true value $\Delta_{i-1} = \beta_{i-1} - \beta_*$, but not directly through the location of $\beta_{i-1}$. This guarantees that we can develop asymptotic theory by assuming that $\beta_* = 0$ without loss of generality.

remark[IV Clustered Dependence] We can extend IV model (ref) to a cluster-dependent setting: \begin{align} y_{i,t} = x_{i,t}'\beta_* + u_{i,t}, \mathbb{E}\left[u_{i,t} z_{i,t} \right]=0, \quad t=1, \ldots, T_i, i=1, \ldots, n, \end{align} where there are $n\to\infty$ clusters, and within each cluster we have finite many ($T_i$) observations. It is straightforward to accommodate clustered dependence. Specifically, we now update the estimator at the $i$ level: run a modified version of (ref) with \begin{align*} G_i = \frac{1}{T_i} \sum_{t=1}^{T_i} z_{i,t} x_{i,t}' \; and \; H_i = - \frac{1}{T_i} \sum_{t=1}^{T_i} z_{i,t} y_{i,t}. \end{align*} In words, instead of updating the estimate for each observation, we treat all individuals within a cluster as a “mini-batch” and update the estimate in batches. This allows for arbitrary dependence within clusters.

remark[Two Sample IV] It is easy to accommodate our algorithm for the case when $G_i$ and $H_i$ are from two different datasets.

Asymptotic Properties

We first state basic regularity conditions.

asm*\begin{enumerate}[({A}1)] • $\mathbb S = \{\mathcal{D}_i = (x_i, z_i, y_i) \in \mathbb R^{d_\beta} \times \mathbb R^{d_g} \times \mathbb R: i=1,2,\ldots,n \}$ is a random sample of size $n$ drawn from model (ref), with $d_\beta \leq d_g <\infty$. • For some fixed $n_0 \in \mathbb N$, $\mathbb S_0 = \{ \mathcal{D}_{0j} = (x_{0j}, z_{0j}, y_{0j}) : j=1,\ldots, n_0\}$ is an initialization random sample drawn from model (ref), set aside from $\mathbb S$. • $\lambda_{\min}(G' G) > 0$ with $G = \mathbb E [G_i]= \mathbb E [z_i x_i']$, and $\lambda_{\min}(\mathbb{E} [ z_i z_i']) > 0$. • $\beta_* \in \mathbb R^ {d_\beta}$ is the (unique) solution to $\mathbb E[g_i(\beta)] = \mathbb E[ z_i \left( x_i'\beta - y_i \right)] = 0$. • $\gamma_i = \gamma_0 i^ {-a}$ for some $\gamma_0 > 0$ and $a \in (1/2, 1)$. • $\mathbb E \|G_i\|^ {2} < \infty$, and $\mathbb{E}[\|g_i(\beta_*)\|^2] < \infty$. • $\mathbb E\|G_i\|^{2p} < \infty$ and $\mathbb E \|g_i(\beta_*)\|^{2p}<\infty$ for some integer $p > (1 - a)^{-1}$. \end{enumerate}

Assumption (ref) is an initialization sample that is used to construct $(\beta_0,\Phi_0,W_0)$. Condition (ref) amounts to identification conditions and equivalent to the conditions that $G$ has full rank $d_\beta$ and $\mathbb{E} [ z_i z_i']$ is non-singular. Condition (ref) defines $\beta_*$ and its uniqueness is guaranteed by (ref). Assumption (ref) is the standard condition for the learning rate in the literature polyak1992acceleration. Conditions (ref)-(ref) impose moment conditions: (ref) is a less stringent assumption that ensures that the non-averaged estimator $\beta_n$ is strongly consistent for $\beta_*$ and its $L_2$ convergence rate is $O(\gamma_n)$. It is also used to obtain asymptotic normality of the averaged estimator $\bar \beta_n$, which converges faster than $\beta_n$; (ref) is a more stringent condition under which we obtain the functional central limit theorem (FCLT) for the sequence of S2SLS estimators $\beta_i$.

The following lemma establishes strong consistency of $\beta_n$.

lem[Strong consistency] Let Assumptions (ref)-(ref), (ref) hold, and $\gamma_i=\gamma_0 i^{-a}$ for $a\in (1/2,1]$. Then, as $n \to \infty$, $\beta_n \to \beta_*$ almost surely.

Lemma (ref) is non-trivial to prove because our proposed algorithm is based on the Z-estimator, not on the M-estimator. The proof of Lemma (ref) relies on martingale techniques: in particular, Robbins and Siegmund robbins1971convergence, which provides a convergence theorem for non-negative “almost supermartingales.”\footnote{ See, e.g., Chapter 5 of benveniste2012adaptive for an application of the Robbins-Siegmund theorem to the Robbins-Monro algorithm Robbins-Monro-1951.} Moreover, as for the original Robbins-Monro algorithm, the almost sure convergence of $\beta_n$ allows for the learning rate of $\gamma_i =1/i$.

We now present asymptotic normality of the averaged estimator $\bar\beta_n$.

thmLet Assumptions (ref)-(ref) hold. Then, as $n\to \infty$, \begin{equation*} \sqrt {n} (\bar \beta_{n} - \beta_*) \overset{d}{\to} N(0, \operatorname{Avar}(\bar \beta)), \end{equation*} where $\Omega \equiv \operatorname{var}(g_i(\beta_*))$, $W = ( \mathbb{E} [ z_i z_i'] )^{-1}$ and $\operatorname{Avar}(\bar \beta) \equiv (G' W G)^{-1}G' W \Omega W G(G'W G)^{-1}$.

To prove Theorem (ref), it is necessary to extend Theorems 1 and 2 of polyak1992acceleration to accommodate the additional dynamics due to $\Phi_{i}$ and $W_i$. Since $\Phi_{i} \ne G$ in general, we must carefully consider the error $\Phi_{i} - G$. It is central to control this error to obtain asymptotic normality.

remarkThe limiting distribution of our averaged S2SLS is the same as that of the standard (offline) 2SLS estimator. We will propose an efficient estimator in Section (ref).

remarkIf $x_i$ is exogenous in model (ref), then we can take $z_i = x_i$, $G = \mathbb{E}[ x_i x_i']$, $W = (\mathbb{E}[ x_i x_i'])^{-1}$, and $\Omega = \mathbb{E}[ u_i^2 x_i x_i']$. This implies that $\operatorname{Avar}(\bar \beta) = \mathbb{E}[ x_i x_i']^{-1} \mathbb{E}[ u_i^2 x_i x_i'] \mathbb{E}[ x_i x_i']^{-1}$, which is exactly identical to the asymptotic variance for the standard (offline) OLS estimator. In other words, when $x_i$ is exogenous, our S2SLS estimator is not algebraically equivalent to the standard SGD-OLS estimator, but it is first-order asymptotically equivalent to it.

We now strengthen Theorem (ref) to the following functional central limit theorem (FCLT).

thmLet Assumptions (ref)-(ref) and (ref) hold. Then: as $n\to \infty$, \begin{equation*} \left\{\frac 1 {\sqrt {n}} \sum_{i=1}^{\lfloor n r\rfloor}(\beta_i - \beta_*) \right\}_{r\in[0,1]}\rightsquigarrow \operatorname{Avar}(\bar \beta)^{1/2}\{ \mathbb{W}_{d_\beta}(r)\}_{r\in [0,1]}, \end{equation*} where $\{\mathbb{W}_{d_\beta}(r)\}_{r\in [0,1]}$ denotes the $d_\beta$-dimensional standard Wiener process.

The FCLT in Theorem (ref) states that the partial sum of the sequentially updated estimates $\beta_i$ converges weakly to a rescaled Wiener process, with the scaling matrix equal to a square root of the asymptotic variance of $\bar \beta_{n}$. Note that Theorem (ref) is a special case of Theorem (ref) with $r=1$ (albeit Theorem (ref) is derived under milder moment conditions). Theorem (ref) allows us to construct robust online confidence regions for $\beta_*$; see Subsection (ref) below.

comment\subsection{Further Discussion} One might consider improving the asymptotic performance of the estimation by applying \begin{equation*} \beta_i = \beta_{i-1} - \gamma_i \Phi_{i-1}' W_{i-1} (G_i\beta_{i-1} + H_i) \end{equation*} where $W_{i-1}$ is presumably either a good approximation of the optimal weighting matrix $\Omega^{-1}$ or achieves better balancing across the moment conditions. As a simple case, suppose $W_{i-1} \equiv W_0 \in \mathbb R^{d_g \times d_g}$ is an $\mathcal F_0$-measurable symmetric matrix satisfying $\mathbb P(C^{-1} \le \lambda_{\min}(W_0)\le \lambda_{\max}(W_0) \le C) = 1$ for some $C >1$. Since our procedure is conditioning on $\mathcal F_0$, the method allows an updating rule pre-conditioned on $W_0$. We may apply our method to the redefined data set $\check G_i = W_0^{1/2} G_i$, $\check H_i = W_0^{1/2} H_i$ while maintaining the validity of the assumptions and that of the theorems by conditioning on $\mathcal F_0$. For instance, this embraces the use of studentized moment conditions or a pre-estimate of $\Omega^{-1}$ with possible regularization based on the initial sample $\mathbb S_0$. If we could obtain a good online estimator $W_{i-1}$ of $\Omega^{-1}$ with a guarantee of regularity, it may lead to more precise or near-optimal estimation. We did not pursue this direction explicitly in this paper for the following reasons. First, since one of our core motivations is fast computation of large-scale models (in an online fashion), it may not be computationally feasible to update a large matrix $W_{i-1}$ along the way, which could potentially cause a bottleneck in computation. Second, if $d_g$ is too large, it is not straightforward to get a regularized estimate of a large covariance matrix online. This poses another concern of different nature, which might elicit future research. As a result, it is ambiguous whether the inclusion of online optimal weighting matrix could lead to a substantial improvement outweighing these potential issues. Lastly, it may be worth noting that $\Phi$ is much more essential ingredient to operationalize our procedure because the moment conditions themselves are directionless and may be overidentifying the model. $\Phi$ reduces the dimension by picking up $d_\beta$ linear combinations out of $d_g$ moment conditions and informs the direction to the true value of the parameter, while $W$ only orthogonalizes \todo[color=red!20!white,fancyline,size=\tiny]{We may experiment on this. e.g. we may check the degree of slowdown of the MSE of the path as we increase the correlation between the IVs while maintaining the conditional homoskedasticity.} the moment conditions possibly intertwined. {\color{blue} [Matt's ideas] I am curious if we may use something like $W_{i+1} = (\frac{1} { n_0 + i} \sum_j g_j (\bar{\beta_i}) g_j (\bar{\beta_i})' + a_i I )^{-1} $ for $a_i \to 0$ suitably? We may set the rate for $a_i$ quite slow to reduce the computational burden? Or a solution path of $W$ solving $Evech(I-g_i g_i' W) $=0 with a diagonal matrix set as the initial? }

Efficient Estimation

SGMM Algorithm

In general, the S2SLS estimator in (ref) is not efficient for $\beta_*$, just like the standard 2SLS estimator is inefficient. To obtain an efficient estimator of $\beta_*$, we now propose to implement the following procedure. First, we randomly partition the main sample into two subsamples: $\mathbb{S} = \mathbb{S}_1 \cup \mathbb{S}_2$, where the sample size of $\mathbb{S}_j$ is denoted by $n_j$, where $j=1,2$. Thus, $n = n_1 + n_2$. Using $\mathbb{S}_1$, we run (ref) until $i=n_1$, and then using $\mathbb{S}_2$, we sequentially update from $i = n_1 + 1$ until $i = n$:

subequations\begin{align} \beta_i &= \beta_{i-1} - \gamma_i (\Phi_{i-1}' W_{i-1} \Phi_{i-1})^{\dagger} \Phi_{i-1}' W_{i-1} g_i(\beta_{i-1}), \\ \Phi_i &= \frac{n_0 + i - 1}{ n_0 + i} \Phi_{i-1} + \frac{1} { n_0 + i} G_i, \\ m_i&=n_0+i-1+g_i(\bar\beta_{n_1})'W_{i-1} g_i(\bar\beta_{n_1}), \\ W_{i}&= \frac{n_0+i}{n_0+i-1} W_{i-1} \left[I- m_i^{-1}g_i(\bar\beta_{n_1})g_i(\bar\beta_{n_1})'W_{i-1} \right], \\ \bar \beta_i &= \frac{i-1 }{ i} \bar \beta_{i-1} + \frac{1} { i}\beta_i. \end{align}

To achieve efficiency, we assume that $n_1 \rightarrow \infty$ but $n_1 / n \rightarrow 0$. In practice, the iterations up to $n_1$ can be viewed as a “warm-up” stage to avoid any too abrupt path in $\beta_i$.

remarkNote that our efficient algorithm (ref) is virtually the same as the inefficient algorithm (ref), except that we now update the weighting matrix $W$ differently, aiming for the optimal weighting $W=\Omega^{-1} \equiv [\operatorname{var}(g_i (\beta_*))]^{-1}$. As in Remark (ref), we apply SMW formula to sequentially update $W=\Omega^{-1}$ in (ref). Also, note that we keep the same $\bar\beta_{n_1}$ in (ref) and (ref), which is a consistent estimator for $\beta_*$.

remark[Computation of $\mathcal{V}_{i-1} =(\Phi_{i-1}' W_{i-1} \Phi_{i-1})^{\dagger}$ for efficient estimation] As in Remark (ref), it would be demanding to directly compute $\mathcal{V}_{i-1} =(\Phi_{i-1}' W_{i-1} \Phi_{i-1})^{\dagger}$ in (ref) when $d_\beta$ is large. Analogous to (ref), we use SMW formula: for $\mathcal{V}_i := (\Phi_{i}'W_{i} \Phi_{i})^{-1} \in \mathbb{R}^{d_{\beta} \times d_{\beta}}$, \begin{equation*} \mathcal{V}_i = \frac {n_0 + i} {n_0 + i-1} \left( \mathcal{V}_{i-1} - \mathcal{V}_{i-1} \mathcal{U}_i (\mathcal{D}_i + \mathcal{U}_i'\mathcal{V}_{i-1}\mathcal{U}_i )^{-1} \mathcal{U}_i' \mathcal{V}_{i-1}\right), \end{equation*} where \begin{align*} \mathcal{D}_i &=\operatorname{diag}\left[-z_i'W_{i-1} z_i, z_i'W_{i-1} z_i, -m_i\right] \in \mathbb{R}^{3 \times 3},\\ \mathcal{U}_i &= \left[\Phi_{i-1}' W_{i-1} z_i, \Phi_{i-1}' W_{i-1} z_i + \frac{z_i'W_{i-1}z_i}{n_0 + i-1}x_i, \Phi_{i-1}' W_{i-1} g_i(\bar\beta_{n_1}) + \frac{z_i'W_{i-1}g_i(\bar\beta_{n_1})}{n_0 + i-1}x_i\right] \in \mathbb{R}^{d_{\beta} \times 3}. \end{align*} Since $\mathcal{D}_i + \mathcal{U}_i' \mathcal{V}_{i-1} \mathcal{U}_i$ is a 3 by 3 matrix, inverting it would require much less computation compared to directly inverting $\Phi_{i-1}' W_{i-1}\Phi_{i-1}$ when $d_\beta \gg 3$.

Asymptotic Efficiency

We make the following additional regularity condition.

asm*\begin{enumerate}[({A}8)] • $n_1 \to \infty$, $n_1/n \to 0$, $\mathbb{E}[\|\beta_0\|^{2\tilde{p}} \,] < \infty$ for some constant $\tilde{p}$, and $\inf_{\beta \in \mathcal{K}} \lambda_{\min}(\mathbb{E}[g_i(\beta)g_i(\beta)']) > 0$ for some compact set $\mathcal{K}$ that contains $\beta_*$ in its interior. \end{enumerate}

The following theorems establish asymptotic properties of SGMM.

thmLet Assumptions (ref) -- (ref) and (ref) hold with $\tilde{p}=1$. Then, as $n \to \infty$, $\beta_n$ and $\bar \beta_n$ are weakly consistent for $\beta_*$ and \begin{equation*} \sqrt{n}(\bar \beta_n - \beta_*) \overset{d}{\to} N(0, (G' \Omega^{-1} G)^{-1}). \end{equation*}

Theorem (ref) shows that SGMM is asymptotically first-order equivalent to the standard (offline) efficient GMM estimator. Theorem (ref) below strengthens Theorem (ref) by establishing the FCLT under extra moment condition.

thmLet Assumptions (ref) - (ref), (ref) and (ref) hold with $\tilde{p}=p$, where $p$ is defined in Assumption (ref). Then, it holds \begin{equation*} \left\{ \frac 1 {\sqrt n}\sum_{i=1}^{\lfloor nr \rfloor} (\beta_i - \beta_*) \right\}_{r\in[0,1]}\rightsquigarrow (G'\Omega^{-1}G)^{-1/2} \{ \mathbb{W}_{d_\beta}(r) \}_{r \in [0, 1]}. \end{equation*}

Theorem (ref) suggests robust online confidence regions as those presented in Subsection (ref) below.

comment{\color{blue}[MS] The case with $\zeta > 0$ seems a bit hard to deal with because we have a time-varying dynamics \begin{align*} \bar W(t)&:=\operatorname{plim}_{n\to\infty}W_{[nt]} \\ &= \begin{cases} Q^{-1} & t \le \zeta, \\ t(\zeta Q + (t - \zeta) \Omega)^{-1} & \zeta < t \le 1. \end{cases} \end{align*} The asymptotic variance of $\bar\beta_n$ would be an average of $$\{ (G'\bar W(t)G)^{-1} G'\bar W(t) \Omega \bar W(t) G (G'\bar W(t) G)^{-1}\}_{0\le t\le 1}$$ along the time path $0\le t\le 1$, but as it varies with time, it seems hard to establish this with the current argument. It would also affect Theorem (ref) because $\sqrt n(\bar{\beta}_{\lfloor nr \rfloor} - \beta_*)$ would have a time-varying variance factor. Perhaps it would be more tractable if we assume $n_1 = o(n)$. } {\color{blue} [Idea for proof] Assume that the smallest eigenvalue of $\mathbb{E} [g_i(\beta) g_i(\beta)' ]$ is bounded from below by a strictly positive constant in a neighborhood of $\beta_*$. Then, consistency of $\bar{\beta}_{n_1}$ can be used to ensure that the smallest eigenvalue of the weighting matrix is bounded away from zero. We expect that the asymptotic variance has the standard form as in Theorem 1 but the weighting matrix is a weighted average between the 2SLS weighting matrix and the efficient weighting matrix [to be checked] }

Inference

In this section, we first present two simple methods to construct fast online confidence regions. We then show that a couple of well-known statistical tests can be seamlessly integrated within the SGMM framework.\footnote{The purpose of this section is to showcase the usefulness of our approach. It is a topic for future research to investigate a variety of inference problems more extensively.}

Online Confidence Regions

We propose two simple online confidence regions: the plug-in (PI) base approach and the random scaling (RS) approach.

Plug-in consistent online confidence regions

As a by-product of the efficient algorithm, $W_n$ defined in (ref) consistently estimates $\Omega^{-1}=(\operatorname{var}[g_i(\beta_*)])^{-1}$, and the $\mathcal{V}_n=(\Phi_{n}' W_{n} \Phi_{n})^{\dagger}$ defined in Remark (ref) consistently estimate the asymptotic efficient variance $(G' \Omega^{-1} G)^{-1}$ in Theorem (ref). Hence, we can conduct asymptotic optimal inference using $\mathcal{V}_n =(\Phi_{n}' W_{n} \Phi_{n})^{\dagger}$, and the resulting inference will be called “plug-in inference”. In particular, we can bulid the optimal plug-in online confidence regions based on

equation[equation omitted — 170 chars of source]

Random scaling robust online confidence regions

For both inefficient S2SLS and efficient SGMM, given the FCLT Theorems (ref) and (ref), we can apply the random scaling method proposed in lee2021fast, lee2022fastqr, which is based on the following robust, inconsistent long-run variance (LRV) estimate idea of kiefer2000simple,velasco2001edgeworth,gupta2021robust for $\operatorname{Avar}(\bar \beta)$:

equation[equation omitted — 206 chars of source]

See lee2021fast for simple online version to compute $\widehat V_{rs,n}$ sequentially. Then a robust online confidence region for $\beta_{*}$ can be constructed using the following statistic:

align[align omitted — 406 chars of source]

where the critical values can be simulated as in kiefer2000simple.

We have implemented online confidence sets using $\mathcal W_{eff}$ and $\mathcal W_{rs}$ versions in Monte Carlo experiments and real data applications below.

Online Endogeneity Tests

To further illustrate the usefulness of our inference method, we now consider an endogeneity test focusing on only a subset $\beta_{\text{sub}}$ of $\beta_*$. Under the null, the probability limits of OLS and IV estimators are the same; under the alternative, the IV estimator is still consistent but OLS is not. Let $\alpha_{\operatorname{sub}}$ denote the probability limit of OLS for the subvector. The null hypothesis is then

align[align omitted — 103 chars of source]

We propose an online algorithm to implement the Durbin-Wu-Hausman (DWH) test. Let $\beta_i$ and $\alpha_i$ respectively denote the stochastic sequences of the IV-estimator and OLS. They are jointly updated as follows:

equation[equation omitted — 283 chars of source]

remarkNote that $\alpha_i$ follows the usual SGD path for M-estimation. To improve the finite-sample performance, we could have multiplied $x_i(x_i'\alpha_{i-1}-y_i)$ by $(\frac{1}{n_0+i-1}\sum_{j\le n_0} x_{0j}x_{0j}' + \frac{1}{n_0+i-1}\sum_{j\le i-1} x_j x_j')^{-1}$.

Let $\mbox{\boldmath $\beta$}_i=(\beta_i',\alpha_i')'$ denote the vector stacking all elements of the updating sequences and let $\bar\mbox{\boldmath $\beta$}_n=\frac{1}{n}\sum_{i=1}^n\mbox{\boldmath $\beta$}_i$. Let $$\mbox{\boldmath $\beta$}_{*}=

pmatrix[pmatrix omitted — 65 chars of source]

. $$ We show that under either the null or the alternative, for some covariance matrix $\Gamma$,

equation*[equation* omitted — 238 chars of source]

The above FCLT allows us to construct a simple online DWH test for endogeneity. The test statistic has to be properly scaled using the asymptotic variance. While the scaling asymptotic variance in the DWH test stems from the idea of estimation efficiency, it is computationally demanding to implement its exact form via stochastic approximation in the online context. We adopt the above random scaling robust LRV estimate $\widehat{V}_{rs,n}$ in (ref) instead. In particular we use $$ \widehat{V}_{n}:=\frac{1}{n}\sum_{s=1}^n\left(\frac{1}{\sqrt{n}}\sum_{i=1}^s(\mbox{\boldmath $\beta$}_i-\bar\mbox{\boldmath $\beta$}_n)\right)\left(\frac{1}{\sqrt{n}}\sum_{i=1}^s(\mbox{\boldmath $\beta$}_i-\bar\mbox{\boldmath $\beta$}_n)\right)'. $$ See lee2021fast for updating ($\bar\mbox{\boldmath $\beta$}_n,~\widehat V_{n}$) sequentially.

Let $\mbox{\boldmath $\beta$}_{\operatorname{sub},i}$ denote the subvector of $\mbox{\boldmath $\beta$}_i$, corresponding to $\beta_{\operatorname{sub}}$ and $\alpha_{\operatorname{sub}}$. In the algorithm, $\mbox{\boldmath $\beta$}_i$, $\Phi_i$ and $W_i$ are potentially high-dimensional objects. Instead of sequentially update the full vector/matrix $\bar\mbox{\boldmath $\beta$}_i$ and $\widehat{V}_{i}$, we just need to update the subvector $\bar\mbox{\boldmath $\beta$}_{\operatorname{sub},i}$ and its corresponding submatrix $\widehat{V}_{\operatorname{sub},i}$.

Note that we can express $\bar\mbox{\boldmath $\beta$}_{\operatorname{sub},n}=(\bar\beta_{\operatorname{sub},n}',\bar\alpha_{\operatorname{sub},n}')'$ corresponding to the online IV and OLS estimators. The online DWH test is then conducted by comparing $\bar\beta_{\operatorname{sub},n}$ and $\bar\alpha_{\operatorname{sub}, n}$, which can be expressed as $$ \bar\beta_{\operatorname{sub},n}-\bar\alpha_{\operatorname{sub}, n}= \Xi \bar\mbox{\boldmath $\beta$}_{\operatorname{sub},n},\text{ where } \Xi =(I, -I). $$ Let $q$ denote the number of restrictions in the null hypothesis (ref). The pivotal statistic is now defined as $$ \mathcal{S}_{rs}:= \frac{n}{q}(\bar\beta_{\operatorname{sub},n}-\bar\alpha_{\operatorname{sub}, n})'(\Xi\widehat{V}_{\operatorname{sub},n}\Xi')^{-1} (\bar\beta_{\operatorname{sub},n}-\bar\alpha_{\operatorname{sub}, n}). $$ The asymptotic distribution of the pivotal statistic can be derived using the FCLT of the stacked vector. This implies the asymptotic null distribution of the pivotal statistic, stated as follows.

corLet assumptions for Theorem (ref) and Theorem (ref) hold, respectively. Under the null hypothesis that $\alpha_{\operatorname{sub}}=\beta_{\operatorname{sub}}$, we have: $$ \mathcal{S}_{rs} \rightarrow\!_{d} ~\frac{1}{q} \mathbb{W}_{q}(1)'\left(\int_0^1[\mathbb{W}_{q}(r) -r \mathbb{W}_{q}(1)][\mathbb{W}_{q}(r) -r \mathbb{W}_{q}(1)]'dr\right)^{-1}\mathbb{W}_{q}(1), $$ where $q = \mathrm{dim}(\alpha_{\operatorname{sub}})=\mathrm{dim}(\beta_{\operatorname{sub}})$.

Critical values for testing linear restrictions are given in kiefer2000simple.

Online Sargan-Hansen Tests

As a straightforward corollary to the main result in Section (ref), which proposes an online efficient GMM estimation, we also implement the test for the overidentifying moment restrictions. Let $\hat{g}_{n_{1}}=\frac{1}{n_{1}}\sum_{i=1}^{n_{1}}g_{i}\left(\bar{\beta}_{n_{1}}\right)$ and add the following to the end of the efficient algorithm (ref): from $i = n_1 + 1$ until $i = n$,

eqnarray*[eqnarray* omitted — 107 chars of source]

Then, we obtain the conventional chi-squared test for the overidentifying restrictions.

corLet Assumptions in Theorem (ref) hold with $d_g>d_{\beta}$. Then we have: \[ \mathcal {J}_{n}:=n\hat{g}_{n}'W_{n}\hat{g}_{n}\overset{d}{\to}\mathcal{X}_{d_{g}-d_{\beta}}, \] as $n\to\infty$, where $W_n$ is defined in (ref).

Monte Carlo Experiments

In this section, we investigate the numerical performance of the SGMM estimator via Monte Carlo experiments. Initially, we discuss the process of selecting the learning rate $\gamma_1$, which will be useful when working with a real data set.

Selection of the Learning Rate in Applications

In this section, we describe a rule of thumb regarding how to choose $\gamma_i = \gamma_0 i^ {-a}$. Suppose that $a \in (1/2,1)$ is fixed at a given constant (in the examples reported below, we set $a=0.501$). Then, it remains to choose $\gamma_0 > 0$, that is, the initial value of the learning rate. Recall that we have the initialization sample $\mathbb S_0$ with sample size $n_0$ to compute $\beta_0$, $\Phi_0$ and $W_0$ and start with the first update as

align[align omitted — 140 chars of source]

where $\gamma_1 = \gamma_0$. We first define

align[align omitted — 222 chars of source]

where $\| \cdot \|_2$ is the spectral norm. We propose to use

align[align omitted — 76 chars of source]

where $\alpha$ is a predetermined quantile level (e.g., $\alpha = 0.5$). The rational behind this rule of thumb is that we choose $\gamma_0$ small enough such that it is likely that the $\beta_i$ path is not explosive when $i$ is relatively small.

Simultation Results

We consider the following data generating process as a baseline model:

align[align omitted — 98 chars of source]

where $x_i$ is a $p$-dimensional vector of regressors, with the first element being endogenous. There exists a $q$-dimensional vector of exogenous variables $z_i$, which follows a multivariate normal distribution $N(0,\Sigma)$. The $(i,j)$ element of $\Sigma$ is set to be $\Sigma_{i,j}:=\rho^{|i-j|}$. The endogenous regressor $x_{i,1}$ is generated as follows: for some $\underline{p} \le p$ and $\underline{q} \le q$,

align[align omitted — 141 chars of source]

where $x_{i,j} = z_{i,j-1}$ for $j=2,\ldots,\underline{p}$ and $\nu_i \sim N(0,1)$. Finally, the error term in (ref) is generated by

align[align omitted — 65 chars of source]

where $\sigma_i = 5 \cdot \exp(z_{i,\underline{q}})$ and $\eta_i\sim N(0,1)$. Therefore, the model allows for both heteroskedasticity and endogeneity.

We consider four different sample sizes $n=\{10^4, 10^5, 10^6,10^7\}$. We set the correlation coefficient of $z_i$ as $\rho=0.5$ and the true regression coefficients as $\beta_*=(1,\ldots,1)$. The dimensions of $x_i$ and $z_i$ are set to $(p,q)=(5,20)$ and $(10,25)$ with $(\underline{p},\underline{q})= (5,20)$. Therefore, we conduct the Monte Carlo experiments over 8 different designs. We replicate each design 1,000 times to compute the performance statistics.

The simulations are conducted using the Graham cluster of the Digital Research Alliance of Canada, which consists of several Intel CPUs (Broadwell, Skylake, and Cascade Lake) operating at frequencies between 2.1GHz and 2.5GHz. The memory budget is set to 64 gigabytes of RAM.

Tables (ref)--(ref) summarize the simulation results. We estimate the model using two different weight schemes, as described in Sections (ref) and (ref). We denote them as S2SLS and SGMM, respectively. To compare the performance, we also estimate the model using the offline counterparts: 2SLS and GMM through R packages ivreg (CRAN version 0.6.2) and gmm (CRAN version 1.7.0), respectively.

For S2SLS and SGMM, we need to set some tuning parameters and initial values. The learning rate $\gamma_i \equiv \gamma_0 i^{-a}$ is set with $a=0.501$ and $\gamma_0$ as the rule of thumb method described in section (ref). This size of an initialization sample is set to $n_0=1000$. Using the initialization sample, we estimate the initial value $\beta_0$ by 2SLS, $\Phi_0=n_0^{-1}\sum_{j=1}^{n_0} G_{0j}$, and $W_0 = (n_0^{-1}\sum_{j=1}^{n_0} z_{0j}z_{0j}')^{-1}$. Finally, we fix $n_1=10\sqrt{n}$ for SGMM. Two alternative methods for inference are considered: SGMM RS (specifically, as in $\mathcal W_{rs}$ in (ref)) and SGMM PI, respectively, refer to random scaling (RS) and plug-in (PI) inference with the same point estimate SGMM.

In the tables, we focus on the coefficient of the endogenous regressor $x_{i,1}$ and report the following performance statistics: root mean square error (RMSE), average bias (Bias), standard deviation (SD), coverage probability of the 95% confidence interval (Coverage Prob), the average confidence interval length (CI Length), and the average computation time in seconds (Time).

table[table omitted — 2,196 chars of source]

table[table omitted — 1,663 chars of source]

Overall, the numerical performance of S2SLS and SGMM is satisfactory. First, both S2SLS and SGMM demonstrate good coverage probabilities across all designs. Additionally, other measures such as RMSE, Bias and SD also indicate good performance. When we examine RMSEs specifically, they are slightly larger than those of the offline estimators when $n=10^4$ and $10^5$. However, for sample sizes $n\ge10^6$, RMSEs become comparable to those of the offline estimators, aligning with the asymptotic theory in the previous sections.

Second, both S2SLS and SGMM shows substantial gains in computation time as the sample size increases. In the model of $(p,q)=(5,20)$, 2SLS takes 1.65 times longer computation time than S2SLS, and GMM does 7.6 times more than SGMM when $n$ is bigger than $10^6$. We observe a similar pattern in the model of $(p,q)=(10,25)$. Note that we compute the whole matrix $\widehat{V}_n$ in these simulations. If we are interested in the inference of a single parameter, we can improve the result further by focusing on a single element of $\widehat{V}_n$.

Third, SGMM demonstrates efficiency gains over S2SLS across all designs as predicted by asymptotic theory. As we discuss earlier, RMSEs of SGMM are comparable to those of GMM when the sample size is $10^6$. If the sample size is as large as $10^7$, GMM exceeds the memory budget of 64 gigabytes, resulting in a value of `NA' in the tables. These results highlight the computational advantage of SGMM over GMM while maintaining its efficiency property.

Finally, we observe that the average confidence interval length (CI Length) is larger for S2SLS and SGMM RS than their offline counterparts. We observe similar phenomena in linear mean and quantile regression models (see lee2021fast, lee2022fastqr.

comment\section{Estimation of a Dynamic Panel Data Model} Consider the following panel autoregressive model with regressors and a one-way error component structure: \begin{align} y_{i, t}=\alpha_* y_{i, t-1}+x_{i, t}^{\prime} \theta_* + u_i+\varepsilon_{i, t}, \quad t=1, \ldots, T, i=1, \ldots, n. \end{align} where $\beta_* \equiv (\alpha_*, \theta_*')'$ are unknown parameters, $x_{i, t}$ is a vector of regressors, $u_i$ is an individual effect, and $\varepsilon_{i, t}$ is an idiosyncratic error. Further, assume that the errors $u_i$ and $\varepsilon_{i, t}$ are mutually independent, the $\varepsilon_{i, t}$ are serially uncorrelated and mean zero, and the regressors $x_{i, t}$ are pre-determined in the sense that $ \mathbb{E}\left[x_{i, t-s} \varepsilon_{i, t}\right]=0 $ for all $s \geq 0$. Apply the differencing operator to obtain \begin{align} \Delta y_{i, t}=\alpha_* \Delta y_{i, t-1}+ \Delta x_{i, t}^{\prime} \theta_* + \Delta\varepsilon_{i, t} \end{align} for $t=3,\ldots,T$. Arellano and Bond (1991) \begin{align*} G_i = \begin{bmatrix} \begin{pmatrix} y_{i,1} \\ x_{i, 1} \\ x_{i,2} \end{pmatrix} & 0 & & 0 \\ 0 & \begin{pmatrix} y_{i,1} \\ y_{i,2} \\ x_{i, 1} \\ x_{i,2} \\ x_{i, 3} \end{pmatrix} & & 0 \\ & & \ddots & \\ 0 & 0 & & \begin{pmatrix} y_{i,1} \\ \vdots \\ y_{i,T-2} \\ x_{i, 1} \\ \vdots \\ x_{i, T-1} \end{pmatrix} \end{bmatrix} \begin{bmatrix} \begin{pmatrix} \Delta y_{i, 2} & \Delta x_{i, 3}^{\prime} \end{pmatrix} \\ \\ \begin{pmatrix} \Delta y_{i, 2} & \Delta x_{i, 4}^{\prime} \end{pmatrix} \\ \\ \vdots \\ \\ \begin{pmatrix} \Delta y_{i, T-1} & \Delta x_{i, T}^{\prime} \end{pmatrix} \end{bmatrix} \end{align*} and \begin{align*} H_i = \begin{bmatrix} \begin{pmatrix} y_{i,1} \\ x_{i, 1} \\ x_{i,2} \end{pmatrix} & 0 & & 0 \\ 0 & \begin{pmatrix} y_{i,1} \\ y_{i,2} \\ x_{i, 1} \\ x_{i,2} \\ x_{i, 3} \end{pmatrix} & & 0 \\ & & \ddots & \\ 0 & 0 & & \begin{pmatrix} y_{i,1} \\ \vdots \\ y_{i,T-2} \\ x_{i, 1} \\ \vdots \\ x_{i, T-1} \end{pmatrix} \end{bmatrix} \begin{bmatrix} \Delta y_{i, 3} \\ \\ \Delta y_{i, 4} \\ \\ \vdots \\ \\ \Delta y_{i, T} \end{bmatrix} \end{align*} {\color{blue} We can add Blundell-Bond moments as well. We may build a Monte Carlo based on this example....}

Empirical Examples

In this section, we explore two empirical applications to demonstrate the effectiveness of the SGMM estimator. Specifically, we revisit the empirical findings presented in angrist1991does and angrist1998children.

angrist1991does

We re-visit the 2SLS estimate of return to education in column (2) of Table IV in angrist1991does:

align*[align* omitted — 84 chars of source]

where $wage_i$ denotes a weekly wage, $educ_i$ denotes the years of education, and $x_i$ is a vector of 9 cohort dummies. The object of interest is $\beta_1$ representing returns to schooling. A vector of instruments, $z_i$, is constructed by the interaction of quarter of birth and cohort dummies, where $d_z=30$. The model is overidentified in this application, as we have only 11 regressors.

Table (ref) provides a summary of the estimation results. Similar to the simulation studies, we employ five different estimators: 2SLS, GMM, S2SLS, SGMM RS, and SGMM PI. Among the total 247,199 observations, we allocate $n_0=20,000$ into the initialization sample, resulting in $n=227,199$.\footnote{Because of exclusion of the initialization sample, the 2SLS estimate in Table (ref) is slightly different from one reported in column (2) of Table IV in angrist1991does. The latter is $0.0769$ with the standard error of $0.0150$.} Similar to the simulations studies, we set $n_1=10\sqrt{n}$ for SGMM. Finally, we adopt the method described in Subsection (ref) and set $\gamma_0=0.200$. In the table, we present the point estimates of $\beta_1$, along with their corresponding 95% confidence intervals, the lengths of these confidence intervals, and the computation times.

table[table omitted — 740 chars of source]

We observe that both SGMM estimators are computed nearly 12.6 times faster than the corresponding offline GMM estimator, and S2SLS performs slightly than 2SLS in this application. Additionally, the confidence intervals of S2SLS and SGMM are wider than those of their offline counterparts, as we confirmed in the simulation studies.

To explore the gap between the point estimates, we next consider implementing a multi-epoch algorithm within this application. Note that the point estimates of the stochastic methods (S2SLS and SGMM) are around 0.111, while the offline (2SLS and GMM) estimates are around 0.076. In order to assure a stable estimation result, we embark on a multi-epoch approach for the stochastic estimations.\footnote{In machine learning, an epoch refers to a complete pass through a training dataset in an iterative optimization algorithm. In our setting, an epoch corresponds to the use of the all observations to run S2SLS or SGMM. In practice, training a machine learning model typically involves running through multiple epochs. This is because a single pass through the training data might not be sufficient for the estimate to converge to an optimal or near-optimal value. This seems the case with the angrist1991does example.} Precisely, we shuffle the order of observations within the sample during each epoch and subsequently compute S2SLS or SGMM over a series of epochs.\footnote{In other words, within each epoch, we implement uniform sampling without replacement with sample size $n$.}

Figure (ref) demonstrates the estimation path of SGMM PI over 10 epochs along with the pointwise 95% confidence band. The pointwise SGMM plug-in confidence interval was calculated using a sample size of $\min\{i,n\}$, i.e.\ the variance was initially divided by $\sqrt{i}$ for $i=1,\ldots, n$, in the first epoch, and subsequently by $\sqrt{n}$ thereafter. In the graph, we can observe that the estimation process stabilizes after the fifth epoch. The estimation results after the 10th epoch are presented at the bottom of Table (ref), designated as SGMM ME. The point estimate stands at 0.0790, closely aligned with the offline GMM estimate, and the length of the confidence interval is considerably reduced (0.632). Furthermore, SGMM ME necessitates only approximately half the computation time of GMM. Therefore, when the dataset size allows for the storage of all observations, as is the case in both examples, we recommend opting for a multi-epoch algorithm and an assessment of SGMM stability in empirical applications.\footnote{However, it is an interesting open question for future research to formally extend our theory to a multi-epoch (or multi-pass) setting and develop an early stopping rule for the number of epochs.}

figure[figure omitted — 456 chars of source]

angrist1998children

In angrist1998children, they study the effect of childbearing on female labor supply. In our application, we use data consisting of 394,840 observations from the 1980 U.S. census. The dependent variable is the number of working weeks divided by 52; the endogenous regressor is a binary variable that takes value 1 if the number of children is greater than 2; the instrument is a binary variable that takes value 1 if siblings are of the same sex. To be consistent between two applications, we repeat the same exercises as in the previous subsection. Specifically, we set $n_0=20,000$, $n = 374,840$, and $n_1=10\sqrt{n}$. The initial-data-dependent choice $\gamma_0$ was that $\gamma_0 = 0.058$.

Table (ref) and Figure (ref) present the estimation results. As this is a just-identified case, we expect little difference across S2SLS and SGMM, which was empirically verified. It is interesting to notice that the SGMM estimate here basically converges only after one or two epochs, unlike the previous application. This is likely due to the fact that the number of parameters is just two, including the intercept term, and the model is just-identified in the second example.

table[table omitted — 786 chars of source]

figure[figure omitted — 301 chars of source]

Extensions

We conclude the paper by mentioning two possible extensions. First, recall the standard (nonlinear) GMM estimator for the general GMM model, $\mathbb{E}[g_i (\beta_*)]=0$ with $\dim(g_i)\ge d_{\beta}$:

align[align omitted — 88 chars of source]

where $\bar{g}_n (\beta) = n^{-1} \sum_{i=1}^n g_i(\beta)$, and $W_n$ is a weighting matrix that may depend on an initial estimator of $\beta_*$. The first-order condition to (ref) is

align*[align* omitted — 95 chars of source]

The following is a natural efficient online algorithm for nonlinear GMM. We assume that the parameter space for $\beta$ is bounded in this case. For simplicity, we drop the step using the subsample of size $n_1$. Specifically, for $ n_0 $ as before, compute an initial estimate

subequations\begin{align} \beta_0 &= \arg \min_{\beta}\bar{g}_{n_0} (\beta)' \bar{g}_{n_0} (\beta)\\ \Phi_0 &= \frac{\partial}{\partial \beta} \bar{g}_{n_0} (\beta_0) \\ W_0 &= \left( \frac{1}{n_0} \sum_{j=1}^{n_0} g_{j}(\beta_0) g_{j}(\beta_0)' +\eta_0 I \right)^{-1} . \end{align}

Let $ G_i (\beta) = \frac{\partial}{\partial \beta} g_i (\beta) $. We sequentially update from $i = 1$ until $i = n$:

subequations\begin{align} \beta_i &= \beta_{i-1} - \gamma_i ( \Phi_{i-1}'W_{i-1}\Phi_{i-1})^{\dagger} \Phi_{i-1}' W_{i-1} g_i(\beta_{i-1}), \\ \Phi_i &= \Phi_{i-1} - \frac{1}{ n_0 + i}(\Phi_{i-1} - G_i (\bar\beta_{i-1})), \\ m_i &= n_0+i-1+g_i(\bar\beta_{i-1})'W_{i-1} g_i(\bar\beta_{i-1}), \\ W_{i} &= \frac{n_0+i}{n_0+i-1} W_{i-1} \left[ I- m_i^{-1}g_i(\bar\beta_{i-1})g_i(\bar\beta_{i-1})'W_{i-1} \right], \\ \bar \beta_i &= \bar{\beta}_{i-1} - \frac{1}{i} (\bar{\beta}_{i-1}-\beta_i ). \end{align}

We leave it to future research to derive the asymptotic properties of the above nonlinear SGMM. The online Sargan-Hansen test can be computed the same way as that in subsection (ref).

Second, the instruments considered in the paper are assumed to be valid and strong. Allowing many weak IVs in general nonlinear GMM models is an important question that has been fruitfully studied in the past two decades. We expect that overidentified Z-type stochastic approximation is appealing with many weak IVs because it can be helpful to deal with the problem of many local minima. For instance, being a fast algorithm, overidentified Z-type stochastic approximation gives us an avenue of attempting many different starting values. However, we expect that the theoretical studies could be technically challenging, so we leave this extension for future research.