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.
70,988 characters · 20 sections · 79 citation commands
Fast Inference for Quantile Regression with Tens of Millions of Observations
\onehalfspacing
Quantile regression (QR) has been increasingly popular in economics since the groundbreaking work of Koenker1978. Many significant developments have been made to this important statistical methodology. Meanwhile, economists are now able to analyze datasets that contain tens of millions of observations. For instance, one of the important economic applications of quantile regression is to study the wage structure. The 5% sample in the 2000 U.S. Census contains more than 14 million observations and more than 4 million even after restricting the sample to a subpopulation of working \textcolor{black}{White} adults. One motivation of using datasets with such “ultra-large\textcolor{black}{”} sample sizes is to deal with econometric models with a large number of parameters within the framework of simple parametric models. Suppose that there are $d$ parameters to estimate with a sample of size $n$. A parametric regression model can be as large as $d \sim 1,000$ because a large number of controls might be necessary. In the setting of “ultra-large\textcolor{black}{”} sample size with $n\sim 10^7$, straightforward parametric estimators can be used with $d \sim 1,000$ under standard textbook assumptions without demanding conditions such as sparsity assumptions.\footnote{\textcolor{black}{A classical condition for a uniform normal approximation of an $M$-estimator is $ (d \log n)^{3/2}/n \rightarrow 0$ Portnoy:85.}}
With datasets containing tens of millions of observations, however, one major challenge of QR is that it would require huge computing powers and memory that are often not accessible using ordinary personal computers. While inference is desirable in empirical studies to quantify statistical uncertainty, it often \textcolor{black}{demands more computational resources} than point estimation. Typically, asymptotic normal inference with QR estimators involves solving optimization problems and computing asymptotic covariance matrices. Computationally efficient implementation of both tasks is crucial with ultra-large datasets.
In this paper, we focus on this challenge, by studying the QR inference problem in the context of $n\sim 10^7$. We propose a \textcolor{black}{fast statistical inference framework} to analyze cross-sectional data with millions of observations\textcolor{black}{, a typical sample size for the census data}, embracing stochastic gradient descent techniques, one of the most active research areas in machine learning. This literature dates back to Robbins-Monro-1951. While we focus on the QR framework, whose inference has been well known to be a hard problem, our proposed method can be \textcolor{black}{naturally} generalized to other econometric frameworks.
Empirically, we use the data from IPUMS USA IPUMS, a common data source for the U.S. wage structure, to study the trends in the gender gap in terms of the college wage premium. In our empirical application, we aim to control for work experience by flexibly interacting workers' age with gender and state fixed effects, which creates over one thousand regressors, thereby motivating the need of ultra-large datasets. We find that existing inference methods are not applicable due to time or memory constraints, while our inference method provides new insights into trends in the gender gap. There has been a large literature on college wage premium and wage structure. See, for instance, chamberlain1994quantile,Buchinsky:1994,Buchinsky:1998,Gosling:2000,angrist2006quantile,CFM:2013 among others. To the best of knowledge, our empirical illustration is the first to obtain confidence intervals based on quantile regression with more than 4 million observations and more than 1000 regressors.
Consider the setting of a linear quantile regression model for which data $\{ Y_i \equiv (y_i,x_i) \in \mathbb{R}^{1+d} : i=1,\ldots,n\}$ are generated from
where $\beta^*$ is a vector of unknown parameters and the unobserved error $\varepsilon_i$ satisfies $P(\varepsilon_i \leq 0|x_i)= \tau$ for a fixed quantile $\tau \in(0,1)$. Note that $\beta^{*}$ is characterized by \[ \beta^{*}:=\arg\min_{\beta\in\mathbb{R}^{d}}Q\left(\beta\right), \] where $ Q(\beta):= \mathbb E [q(\beta, Y_i)]$ with the check function $ q(\beta, Y_i):= (y_i-x_i'\beta)(\tau-I\{y_i-x_i'\beta\leq 0\}). $ Here, $I(\cdot)$ is the usual indicator function. The standard approach for estimating $\beta^{*}$ is to use the following M-estimator originally proposed by Koenker1978:
See Koenker2005 and Koenker17 for a monograph and a review of recent developments, respectively.
There have been two potential challenges for the inference of standard quantile regression. \textcolor{black}{The first} is the computational problem. The optimization problem is typically reformulated to a linear programming problem, and solved using interior-point algorithms Portnoy:Koenker:97. \textcolor{black}{The second} is associated with statistical inference. As was shown by Koenker1978, the asymptotic distribution of the M estimator is as follows:
where $f_{\varepsilon}(\cdot|x_i)$ is the conditional distribution of $\varepsilon_i$ given $x_i$, assuming data are independent and identically distributed (i.i.d.). Hence, in the heteroskedastic setting, standard inference based on the “plug-in” method would require \textcolor{black}{estimating matrix $H$, which involves the conditional density function $f_{\varepsilon}(\cdot|x_i)$}. An alternative inference would be based on bootstrap.
One of the attempts to solving the computational/inference difficulties is to rely on the smoothing idea, via either smoothing the estimating equation Horowitz:1998 or the convolution-type smoothing FGH,TanWangZhou,HE2021. Both require a choice of the smoothing bandwidth. In particular, the convolution-type smoothing (conquer as coined by HE2021) has received recent attention in the literature because the optimization problem is convex, so it is more scalable. Figure 1 of HE2021 shows that conquer works well for point estimations when the sample size ranges from $10^3\sim 10^6$ with the number of regressors $d \approx n^{1/2}$.
Meanwhile, in terms of inference, both types of smoothed estimators are first-order asymptotically equivalent to the unsmoothed estimator $\widehat{\beta}_n$\textcolor{black}{, and the scale of the numerical studies of these solutions is less ambitious.} For instance, the largest model considered for inference in HE2021 is merely $(n, d) = (4000, 100)$ (see Figure 7 in their paper). Therefore, there is a scalability gap between point estimation and inference in the literature. In this paper, we aim to bridge this gap by proposing a method for fast inference.
We focus on the inference problem at the scale up to \textcolor{black}{$(n, d) \sim (10^7, 10^3)$, } which we refer to as an “ultra-large\textcolor{black}{”} quantile regression problem. This is a possible scale with IPUMS USA, but very difficult to deal with using benchmark inference procedures for quantile regression. To tackle this large-scale problem, we estimate $\beta^*$ via stochastic (sub)gradient descent (S-subGD). Our assumption on the asymptotic regime is that $d$ is fixed but $n$ grows to infinity.
Suppose that we have i.i.d. data from a large cross-sectional survey: $$Y_1,..., Y_n,\quad Y_i=(x_i, y_i),$$ where the ordering of these observations is randomized. Our method produces a sequence of \textcolor{black}{estimators}, denoted by $\beta_i$, which is a solution path being updated as $$ \beta_{i}=\beta_{i-1}-\gamma_{i}\nabla q\left(\beta_{i-1},Y_{i}\right). $$ Here $Y_i$ is a \textcolor{black}{new} observation in the randomized data sequence. Also, $\nabla q\left(\beta_{i-1},Y_i\right)$ is a subgradient of the check function with respect to the current update, and $\gamma_i$ is a predetermined learning rate. Then\textcolor{black}{,} we take the average of the sequence $$ \bar\beta_n:=\frac{1}{n}\sum_{i=1}^n\beta_i, $$ which is called as the Polyak1990-ruppert1988efficient average in the machine learning literature. It will be shown in this paper that $\sqrt{n}(\bar{\beta}_{n} - \beta^*)$ is asymptotically normal with the asymptotic variance the same as that of the standard Koenker1978 estimator $\widehat{\beta}_n$. \textcolor{black}{Therefore,} it achieves the same first-order asymptotic efficiency.
The main novelty of this paper comes from how to conduct inference based on $\bar{\beta}_{n}$. We use a sequential transformation of $\beta_i$'s \textcolor{black}{in a suitable way for on-line updating, and construct asymptotically pivotal statistics.} That is, we studentize $\sqrt{n}\left(\bar{\beta}_{n}-\beta^{*}\right)$ via a random scaling matrix $\widehat V_n$ whose exact form will be given later. The resulting statistic is not asymptotically normal but asymptotically pivotal in the sense that its asymptotic distribution is free of any unknown nuisance parameters; thus, its critical values are easily available. Furthermore, the random scaling quantity $\widehat{V}_{n}$ does not require any additional input other than stochastic subgradient path $\beta_i$, and can be computed \textcolor{black}{fast}.
The main contribution of this paper is computational. We combine the idea of stochastic subgradient descent with random scaling and make large-scale inference for quantile regression much more practical. We provide a discussion of computational complexity and show that our algorithm is comparable to the existing ones in terms of point estimation and is superior to those in terms of inference (especially, subvector inference). See Section (ref) for details. In addition, we demonstrate the usefulness of our method via Monte Carlo experiments. Empirically, we apply it to studying the gender gap in college wage premiums using the IPUMS dataset. The proposed inference method of quantile regression to the ultra-big dataset reveals some interesting new features. First, it shows heterogeneous effects over different quantile levels. Second, we find that the female college wage premium is significantly higher \textcolor{black}{than} that of male workers at the median, which is different from what have been concluded in the literature. In fact, although the S-subGD inference tends to be conservative than the standard normal approximation, it reveals statistically significant results while controlling over $10^3$ covariates to mitigate confounding effects. With more availability of such a large dataset, the S-subGD method will make it possible to obtain convincing empirical evidence for other analyses. Another topic is to examine the relationship between seasonality and birth weight Buckles:Hungerman:2013,currie2013within,Chatterji:22 because the quantiles of birth weight are of interest and birth data from vital statistics contain tens of millions of observations.\footnote{For example, the National Bureau of Economic Research has public use data at \url{https://www.nber.org/research/data/vital-statistics-natality-birth-data}.}
Our estimator is motivated from the literature on online learning, where data were collected in a streaming fashion, and as a new observation $Y_i$ arrives, the estimator is updated to $\beta_i$. Under the M-estimation framework with strongly convex, differentiable loss functions, polyak1992acceleration showed that the Polyak1990-ruppert1988efficient average is asymptotically normal. More recently, fang2018online, chen2020statistical, and zhu2021online studied the statistical inference problem and proposed methods that require the implementation of bootstrap, consistent estimation of the asymptotic variance matrix, or the use of “batch-means” approach, respectively. Besides the smoothness assumption that excludes quantile regression, \textcolor{black}{neither consistent estimation of the asymptotic variance matrix nor implementation of bootstrap is} computationally attractive for datasets as large as those with which this paper is concerned. \textcolor{black}{The batch-means approach does not seem to provide adequate finite sample coverage even if $n$ is sufficiently large lee2021fast.}
The idea of random scaling is borrowed from the time-series literature on fixed bandwidth heteroskedasticity and autocorrelation robust (HAR) inference kiefer2000simple,sun2008optimal,sun2014fixed,lazarus2018har and has been recently adopted in the context of stochastic gradient descent (SGD) in machine-learning problems: online inference for linear mean regression lee2021fast, federated learning li2021statistical, and Kiefer-Wolfowitz methods ChenLai2021 among others.\footnote{\textcolor{black}{li2021statistical and ChenLai2021 applied the idea of random scaling after our previous work lee2021fast appeared on arXiv.}} The current paper shares the basic inference idea with our previous work lee2021fast; however, quantile regression is sufficiently different from mean regression and it requires further theoretical development that is not covered in the existing work. In particular, our previous work is based on polyak1992acceleration, which limits its analysis to differentiable and strongly convex objective functions; however, the check function $q(\beta, Y_i)$ for quantile regression is non-differentiable and is not strongly convex. We overcome these difficulties by using the results given in gadat2022optimal. One could have adopted convolution-type smoothing for S-subGD as the corresponding objective function is differentiable and locally strongly convex. We do not go down this alternative route in this paper as the existence of a sequence of smoothing bandwidths converging to zero complicates theoretical analysis and it would be difficult to optimally choose the sequence of smoothing bandwidths along with the sequence of learning rates. Previously, Pesme:Flammarion:20 studied online robust regression in the context of a Gaussian linear model with an adversarial noise; in fact, their algorithm corresponds to our median regression case (see equation (2) in their paper) but they did not study quantile regression models. gadat2022optimal considered, as their example, the recursive quantile estimator, which corresponds to the quantile regression estimator with an intercept being only the regressor. Neither of these two papers delved into the issue of inference.
In the econometric literature, forneron2021estimation and forneron2022estimation developed SGD-based resampling schemes that deliver both point estimates and standard errors within the same optimization framework. The major difference is that their stochastic path is a Newton-Raphson type which requires computing the Hessian matrix or its approximation. Hence, their framework is not applicable to quantile regression because the check function is not twice differentiable. In other words, for quantile regression inferences with over millions of observations, estimating the inverse Hessian matrix \textcolor{black}{$H^{-1}$ in (ref)} is a task that we particularly would like to avoid.
As an alternative to full sample estimation, one may consider sketching Lee:Ng:2020. For example, Portnoy:Koenker:97 developed a preprocessing algorithm that can be viewed as a sketching method and Yang:QR:ICML proposed a fast randomized algorithm for large-scale quantile regression and solved the problem of size $n \sim 10^{10}$ and $d = 12$ by randomly creating a subsample of about $n = 10^5$. However, the theoretical analysis carried out in Yang:QR:ICML is limited to approximations of the optimal value of the check function and does not cover the issue of inference. As an alternative sketching method, one may employ wang2021optimal and construct a random subsample using data-dependent, nonuniform weights to improve the efficiency of the estimator. In mean regression models, the precision of a sketched estimator depends on the subsample size and is typically worse than that of the full sample estimator Lee:Ng:2020,Lee:Ng:2022. We expect a similar phenomenon for quantile regression models.
Let $a'$ and $A'$, respectively, denote the transpose of vector $a$ and matrix $A$. Let $\| a \|$ denote the Euclidean norm of vector $a$ and $\| A \|$ the Frobenius norm of matrix $A$. Also, let $\ell^{\infty}\left[0,1\right]$ denote the set of bounded continuous functions on $[0,1].$ Let $I(A)$ be the usual indicator function, that is $I(A) = 1$ if $A$ is true and 0 otherwise. For a symmetric, positive definite matrix $S$, let $\lambda_{\min}(S)$ denote its smallest eigenvalue.
In this section, we \textcolor{black}{propose our} inference algorithm. Suppose that we have cross-sectional data $\{Y_1,..., Y_n \}$, where $Y_i=(x_i, y_i)$.
First, we randomize the ordering of the observed data: $$ Y_1=(x_1, y_1),..., Y_n=(x_n, y_n). $$ We start with an \textcolor{black}{initial} estimator, denoted by $\beta_0$, using a method which we shall discuss later. Then produce \textcolor{black}{an S-subGD solution path} which updates according to the rule:
where the updating subgradient depends on a single “next” observation $Y_i$:
This is the subgradient of $q\left(\beta,Y_i\right)$ with respect to $\beta$. Here $\beta_{i-1}$ is the “current” estimate, which is updated to $\beta_{i}$ when the next observation $Y_i$ comes into play. Also $\gamma_i$ is a pre-determined step size, which is also called a learning rate and is assumed to have the form $\gamma_{i}:=\gamma_0 i^{-a} $ for some constants $\gamma_0 > 0$ and $a \in (1/2,1)$. Then, we take the average of the sequence $ \bar\beta_n:=\frac{1}{n}\sum_{i=1}^n\beta_i $ as the final estimator.
Computing this estimator does not require storing the historical data; thus, it is very fast and memory-efficient even when $n\sim 10^7$. In fact, as we shall see below, the computational gain is much more substantial for inference.
It is most common to make inference based on a consistent estimator of the asymptotic variance for $\sqrt{n}(\bar\beta_n-\beta^*)$, which is \textcolor{black}{computationally demanding.} For instance, one \textcolor{black}{needs} to estimate \textcolor{black}{the inverse Hessian matrix $H^{-1}$ in (ref)}. \textcolor{black}{Instead of pursuing asymptotic normal inference using a consistent estimator of the asymptotic variance, we employ asymptotic mixed-normal inference based on a suitable random scaling method through the partial sum process of the solution path as done in the fixed-b standardization.} Motivated by the fact that the solution path $\beta_i$ is sequentially updated, we apply the random scaling approach by defining:
Inference will be conducted based on the standardization using $\widehat V_n.$ Computing $\widehat V_n$ can be \textcolor{black}{efficient} even if $n\sim 10^7$ or larger, since it can be sequentially computed as detailed later. Furthermore, its subvector inference involves only the associated elements, substantially reducing the computational burden.
Once $\bar{\beta}_{n}$ and $\widehat{V}_{n}$ are obtained, the inference is straightforward. For example, for the $j$ th component of $\bar \beta_n$, let $\widehat V_{n,jj}$ denote the $(j,j)$ th diagonal entry of $\widehat V_n$. The t-statistic is then
whose asymptotic distribution is mixed normal and symmetric around zero\textcolor{black}{. We will} formally derive it in the next section. The mixed normal asymptotic distribution for the t-statistic in (ref) is the same as the distribution of the statistics observed in the estimation of the cointegration vector by johansen1991estimation and Abadir:Paruolo:97 and in the heteroskedasticity and autocorrelation robust inference in kiefer2000simple. They are different statistics but \textcolor{black}{converge to} the identical distribution\textcolor{black}{, a functional} of the standard Wiener process as shown by Abadir:Paruolo:02. We can use the t-statistic to construct the $(1-\alpha)$ asymptotic confidence interval for the $j$-th element $\beta_{j}^{*}$ of $\beta^{*}$ by \[ \left[ \bar{\beta}_{n,j} - \textrm{cv} (1- \alpha/2) \sqrt{\frac{\widehat{V}_{n,jj}}{n}}, \; \bar{\beta}_{n,j} + \textrm{cv} (1- \alpha/2) \sqrt{\frac{\widehat{V}_{n,jj}}{n}} \; \right], \] where the critical value $\textrm{cv} (1- \alpha/2)$ is tabulated in Abadir:Paruolo:97. For easy reference, the critical values are reproduced in Table (ref). For instance, the critical value for $\alpha = 0.05$ is 6.747. Critical values for testing linear restrictions $H_{0}: R\beta^{*} = c$ can be found in kiefer2000simple.
While computing $(\widehat V_n,\bar\beta_n)$ is straightforward for datasets of medium sample sizes, it is not practical when $n\sim 10^7$ or even larger size, as storing the entire solution path \textcolor{black}{$\beta_i$ } can be infeasible for memory of usual personal computing devices.
We can construct a solution path for $(\widehat V_n,\bar\beta_n)$ along the way of updating $\beta_i$, so all relevant quantities of the proposed inference can be obtained sequentially. Specifically, let $(\bar\beta_{i-1},\beta_{i-1})$ be the current update when we use up to $i-1$ observations. We then update using the new observation $Y_i$ by:
where the intermediate quantities $A_i$ and $b_i$ are also updated \textcolor{black}{sequentially}: $$ A_i=A_{i-1} +i^2\bar\beta_i\bar\beta_i',\quad b_i= b_{i-1} + i^2\bar\beta_i. $$ Therefore computing $(\bar\beta_n, \widehat V_n)$ can be cast sequentially until $i=n$, and does not require saving the entire solution path of $\beta_i$.
A good initial value could help achieve the computational stability in practice, though it does not matter in the theoretical results. We recommend applying conquer in HE2021 to the subsample (e.g. 5% or 10%) to estimate the initial value, which we adopted in our numerical experiments. Specifically, we start with a smooth quantile regression proposed in FGH and HE2021: $$ \beta_0=\arg\min\sum_{i\in \mathcal S}q_{n}(\beta, Y_i), $$ where $q_n$ is a convolution-type-smoothed checked function, and $\mathcal S$ is a randomized subsample. The conquer algorithm is very fast to compute a point estimate of $\beta_0$ but it is much less scalable with respect to inference (see the Monte Carlo results later in the paper). For the intermediate quantities $A_i$ and $b_i$, we set $A_0 = 0$ and $b_0 = 0$.
Another approach is the burn-in method that applies the S-subGD to some initial observations, and throw them away for the main inference procedure. In our experiments, initial values estimated by conquer outperform those by the burn-in method because conquer \textcolor{black}{tends} to produce high-quality initial solutions.
We introduce a practical guideline for determining the value of $\gamma_i = \gamma_0 i^{-a}$. Let $\widehat{\sigma}$ denote the sample standard deviation of $y_i$ using subsample $\mathcal{S}$. Then, we suggest to fix $a = 0.501$ and select the initial learning rate $\gamma_0$ via
where $\phi(\cdot)$ and $\Phi(\cdot)$ are standard normal density and distribution functions, respectively. For example, if $\tau = 0.5$, then $\gamma_0 \approx 0.798/\widehat{\sigma}$; if $\tau = 0.1$ or $0.9$, then $\gamma_0 \approx 0.585/\widehat{\sigma}$. The rule-of-thumb selection of $\gamma_0$ is heuristically motivated from the asymptotic distribution of the sample quantile when the regression errors are normally distributed.
In most empirical studies, while many covariates are being controlled in the model, of interest are often just one or two key independent variables that are policy-related. In such cases, we are particularly interested in making inference for subvectors, say $\beta^*_1$, of the original vector $\beta^*$. Subvector inference is not often direct target of interest in the quantile regression literature, partially because inference regarding the full vector $\beta^*$ and then converting to the subvector is often scalable up to the medium sample size. But this is no longer the case for the ultra-large sample size that is being considered here.
It is straightforward to tailor our updating algorithm to focusing on subvectors. Continue denoting $\beta_i$ as the $i$-th update of the full vector. While the full vector $\beta_i$ still needs to be computed in the S-subGD, both the Polyak-Ruppert average and the random scale can be updated only up to the scale of the subvector. Specifically, let $\bar\beta_{i, sub}$ denote the subvector of $\bar\beta_i$, corresponding to the subvector of interest. Also, let $\widehat V_{i,sub}$ denote the sub-matrix of $\widehat V_i$; both are updated according to the following rule:
In most cases, steps ((ref))-((ref)) only involve small dimensional objects.
Portnoy:Koenker:97 demonstrated that the computational complexity of the interior point algorithm for quantile regression can be on the order of $\mathcal{O} (n d^3 \log^2 n)$ (see discussions after Theorem 5.1 of Portnoy:Koenker:97). They further proposed an initial phase of preprocessing that effectively reduces computational complexity to $\mathcal{O} (n^{2/3} d^3 \log^2 n) + \mathcal{O} (nd)$ Portnoy:Koenker:97. The basic idea behind preprocessing is that using a preliminary estimate of $\beta^*$, one can remove observations that are above and below the quantile regression line with high probability and add two pseudo-observations called “globs”, resulting in the reduction of the sample size from $n$ to $\mathcal{O} (n^{2/3} )$.\footnote{It is an interesting future research topic to examine whether and how preprocessing can be used in our inference procedure.} These results indicate that the interior point algorithm works well with a large $n$ but not with a large $d$. In addition, Chen:07:JCGS developed a smoothing algorithm that has the computational complexity of $\mathcal{O} (n d^{2.5})$, which reduces the dependence on $d$ from $d^3$ to $d^{2.5}$. In HE2021, conquer is computed using gradient descent (GD), which has the computational complexity of $\mathcal{O}(n d k)$, where $k$ is the number of iterations in GD.
Table (ref) summarizes the computational cost of estimation and inference. Importantly, the computational gain of our stochastic subgradient method is not in terms of point estimation, but of inference. Specifically, our S-subGD method has the computational complexity of $\mathcal{O}(n d)$. If the GD algorithm of conquer converges fast in a finite number of iterations, which seems to be the case empirically, we expect that our algorithm will be comparable---will not be superior---to conquer in terms of point estimation. However, note that GD method can suffer from memory constraints when $n d$ is very large.
However, the comparative advantage is more pronounced in statistical inference, which incurs additional computational costs. Our inference, based on random scaling, entails an additional cost of order $\mathcal{O} (n d^2)$. Nevertheless, it reduces to $\mathcal{O} (n d)$ when the object of inference is a subvector of $\beta^*$, with a dimension of $s = \mathcal{O} (d^{1/2})$. This reduction occurs because we only need to update a submatrix of the random scaling matrix, as discussed in the previous section. Similar gains are achieved when focusing on the inference of each element separately, as is common in applied practice. This is again because we only need to update the diagonal elements, referred to as diagonal random scaling in Table (ref). We emphasize that achieving the same order of computational cost for subvector inference is generally not possible with plug-in inference. This is due to the fact that the plug-in estimator of the sandwich formula $(H^{-1}\Sigma H^{-1})$ of asymptotic variance in (ref) cannot be computed using only diagonal elements or a submatrix.
As for the computational cost of competing methods for QR inference, first, the bootstrap inference will be time consuming as the computational cost will be a multiple of bootstrap draws. For instance, bootstrap inference with conquer has the computational complexity of $\mathcal{O}(b n d k)$, where $b$ is the number of bootstrap draws. Second, it requires access to the entire sample to conduct standard plug-in inference based on asymptotic normality. Thus, the plug-in inference has the additional computational complexity of $\mathcal{O} (n d^2)$ even with the conquer algorithm. This will be faster than bootstrap inference if $k > d$, but cannot be reduced to $\mathcal{O} (n s^2)$ or $\mathcal{O} (n d)$ even for subvector inference.
In short, our algorithm can be as competitive as conquer in terms of computational cost, while offering the advantages of lower memory requirements and faster computations for subvector inference.
In this section, we present asymptotic theory that underpins our inference method. We consider the following QR model: $$ y_{i}= x_i'\beta^*_\tau + \varepsilon_i, $$ where $P(\varepsilon_i\leq 0|x_i)= \tau$ for a fixed quantile $\tau \in(0,1)$. Denote by $\beta^*=\beta_{\tau}^*$ for simplicity.
Conditions (ii) imposes some moment conditions on the conditional density, and is satisfied if both $f_{\varepsilon}(\cdot|x_i)$ and $\frac{d}{d\varepsilon}f_{\varepsilon}(\cdot|x_i)$ are uniformly bounded along with a moment condition on $x_i$. Both conditions (ii) and (iii) are imposed here to establish the consistency of $\bar\beta_n$.
Condition (iii) can be viewed as a local identifiability condition, which ensures that the population loss function is locally strongly convex, and is reasonable as $\beta \mapsto q(\beta, Y_i)$ is convex.
Condition (iv) requires a tail condition on the regressors. It is stronger than the usual moment conditions for standard QR inference, but is required by both consistency and the asymptotic inference for S-subGD. In the literature, it is common to adopt stronger conditions on the regressors than on the error depending on the context. For example, zheng2018 and Chen:Lee:23 assumed the uniform boundedness of each of the covariates; the former paper emphasized that a linear quantile regression model with compact support for the regressors is most sensible to avoid the problem of quantile crossing.
One of the key technical steps in the proof of such one-observation-at-a-time updating, is to show that the nonlinear updating path is linearizable. Specifically, define $\Delta_i= \beta_i- \beta^*$, then the S-subGD updating can be rewritten using the estimation path: \[ \Delta_{i}=\Delta_{i-1}-\gamma_{i}\nabla Q\left(\beta_{i-1}\right)+\gamma_{i}\xi_{i}, \] where $\xi_i= \nabla Q(\beta_{i-1}) - \nabla q(\beta_{i-1}, Y_i)$ is a martingale difference sequence. We then show that the path of $\Delta_{i}$ can be approximated by a linear path: \[ \Delta_{i}^{1}:=\Delta_{i-1}^{1}-\gamma_{i}H\Delta_{i-1}^{1}+\gamma_{i}\xi_{i}\quad\text{and}\quad\Delta_{0}^{1}=\Delta_{0}, \] where $H=\nabla^2 Q(\beta^*)$ is the population Hessian matrix. Such an linearization, as one of the technical steps in the proof, transfers the study of nonlinear SGD/S-subGD problems into linear problems.
Condition (v) gives a requirement for the learning rate, which is standard for online learning. polyak1992acceleration showed that $a<1$ is needed for the learning rate to decrease sufficiently slow so that the asymptotic normality of the Polyak-Ruppert averaging estimator holds. Meanwhile, for nonlinear models as in the QR model, it should not degenerate too slowly so $a>1/2$ is also needed.
We extend lee2021fast (who assumed the sample loss function to be smooth and globally convex) to the quantile regression model, and establish a functional central limit theorem (FCLT) for the aggregated solution path of the S-subGD. Under Assumption (ref), we establish the following FCLT:
where $\Rightarrow$ stands for the weak convergence in $\ell^{\infty}\left[0,1\right]$, $W\left(r\right)$ stands for a vector of the independent standard Wiener processes on $\left[0,1\right]$, and $ \Upsilon := H^{-1}SH^{-1}$, with $ S := \mathbb E [x_i x_i'] \tau(1-\tau) $ and $H := \mathbb E [x_i x_i' f_{\varepsilon}(0|x_i)]$.
The FCLT in (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 the usual quantile regression estimator $\widehat{\beta}_n$. Building on the FCLT, we propose a large-scale inference procedure.
For any $\ell \leq d$ linear restrictions \[ H_{0}: R\beta^{*} = c, \] where $R$ is an $(\ell \times d)$-dimensional known matrix of rank $\ell$ and $c$ is an $\ell$-dimensional known vector, the conventional Wald test based on $\widehat{V}_{n}$ is asymptotically pivotal. We formally state the main theoretical \textcolor{black}{result} in the following theorem.
As an important special case of Theorem (ref), the t-statistic defined in (ref) converges in distribution to the following pivotal limiting distribution: for each $j = 1,\ldots,d$,
where $W_1$ is a one-dimensional standard Wiener process.
Next, consider two quantiles $\tau_1, \tau_2\in(0,1)$, and let $\beta_{\tau_1}^*$ and $\beta_{\tau_2}^*$ denote the corresponding coefficients at these two quantiles. We are interested in testing the heterogeneity at the two quantile levels: $$ H_0: \beta_{\tau_1, sub}^*=\beta_{\tau_2, sub}^*. $$ which is to test whether particular subvectors of the coefficients at the two quantile levels are the same. For instance, the subvector could be the slope of a particular regressor of interest, and we are interested in testing whether its quantile effect stays the same when the quantile level changes from $\tau_1$ to $\tau_2$.
Both subvector coefficients can be iteratively estimated using S-subGD, whose stochastic paths $(\beta_i^1, \bar\beta_{i,sub}^1)$ and $(\beta_i^2, \bar\beta_{i,sub}^2)$ are constructed following the Algorithm ((ref))-((ref)). In particular, the stacked random scaling matrix $\bar V_{i,sub}$ corresponding to the two subvectors $(\bar\beta_{i,sub}^1, \bar\beta_{i,sub}^2)$ is updated as:
In this case, we let $G:=(I, -I)$, and form a test statistic: $$ n(\bar\beta_{n,sub}^1-\bar\beta_{n,sub}^2)' (G\bar V_{n,sub}G')^{-1}(\bar\beta_{n,sub}^1-\bar\beta_{n,sub}^2) . $$ We reject the null if the test statistic is large. The following theorem establishes the asymptotic null distribution.
In this section we investigate the performance of the S-subGD random scaling method via Monte Carlo experiments. The main question of these numerical experiments is whether the proposed method is feasible for a large scale model.
The simulation is based on the following data generating process:
where $x_i$ is a $(d+1)$-dimensional covariate vector whose first element is 1 and the remaining $d$ elements are generated from $\mathcal{N} (0,I_d)$. The error term $\varepsilon_i$ is generated from $\mathcal N(0,1)$, and the parameter value $\beta^*$ is set to be $(1,\ldots,1)$. Without loss of generality, we focus on the second element, say $\beta^*_{(2)}$, of $\beta^*$. That is, we focus on subvector inference: the confidence interval for $\beta^*_{(2)}$ as well as testing the null hypothesis of $\beta^*_{(2)} = 1$. To accommodate a large scale model, the sample size varies in $n \in \{10^5, 10^6, 10^7\}$, and the dimension of $x$ varies in $d\in \{10, 20, 40, 80, 160, 320\}$ in the first set of experiments. Later, we extend the dimension of $x$ into $d\in \{500, 750, 1000, 1500\}$. The initial value $\beta_0$ is estimated by the convolution-type smoothed quantile regression (conquer) in HE2021 and we do not burn in any observations. The learning rate is set to be $\gamma=\gamma_0 t^{-a}$ with $\gamma_0=1$ and $a=0.501$. The simulation results are summarized from $1000$ replications of each design.
We compare the performance of the proposed method with four additional alternatives:
We use the following performance measures: the computation time, the coverage rate, and the length of the 95% confidence interval. Note that the nominal coverage probability is 0.95. The computation time is measured by 10 replications on a separate desktop computer, equipped with an AMD Ryzen Threadripper 3970X CPU (3.7GHz) and a 192GB RAM. Other measures are computed by 1000 replications on the cluster system composed of several types of CPUs. We set the time budget as 10 hours for a single estimation and the RAM budget as 192GB.
Figures (ref)--(ref) summarize the simulation result. To save space, we report the results of $d=20, 80,$ and $320$. We provide tables for all the simulation designs in the appendix. First, we observe that S-subGD is easily scalable to a large simulation design. On the contrary, some methods cannot complete the computation within the time/memory budget. For instance, CONQUER-plugin does not work for any $d$ when $n=10^6$ or $n=10^7$ because of the memory budget. QR shows a similar issue when $n=10^7$ and $d=320$. CONQUER-bootstrap cannot meet the 10-hour time budget for a single replication when $n=10^7$ and $d=320$. They are all denoted as `NA' in the figures.
Second, S-subGD outperforms the alternative in terms of computational efficiency. The bottom panel of Figure (ref) shows the relative computation time, and \textcolor{black}{the relative speed improvement is by a factor of 10 or 100.} In case of CONQUER-plugin, it is slower than S-subGD in the scale of 1000's. If we convert them into the actual computation time, S-subGD takes 16 seconds for the design of $n=10^6$ and $d=320$. However, QR, CONQUER-bootstrap, and SGD-bootstrap take 1374 seconds (22 minutes 54 seconds), 4113 seconds (1 hour 8 minutes 33 seconds), and 1643 seconds (27 minutes 23 seconds) on average of 10 replications.
Third, all methods shows satisfactory results in terms of the coverage rate and the length of the confidence interval. The top panel of Figure (ref) reports that coverage rates are all around 0.95 although S-subGD and SGD-bootstrap are slightly under-reject and over-reject when $d=320$. The average lengths of the confidence interval are reported on the bottom panel. The performance of S-subGD is comparable to QR, CONQUER-plugin, and CONQUER-bootstrap. SGD-bootstrap shows much larger CI lengths, especially when $d=320$.
In the second set of experiments, we stretch the computational burden by increasing the size of $d$ up to $d=1500$ when $n=10^7$. Because of the challenging model scale, we use the cluster system for these exercises and check only the performance of S-subGD focusing on inference of a single parameter. Table (ref) reports the simulation result along with different sizes of $d$. S-subGD still completes the computation within a reasonable time range. The average computation time is under 505 seconds (8 minutes 25 seconds) even for the most challenging design. The coverage rate and CI length becomes slightly worse than those of smaller $d$'s but they still fall in a satisfactory range.
In sum, S-subGD performs well when we conduct a large scale inference problem in quantile regression. It is scalable to the range where the existing alternatives cannot complete the computation within a time/memory budget. In addition to computational efficiency, the coverage rate and CI lengths are also satisfactory. Thus, S-subGD provides a solution to those who hesitate to apply quantile regression with a large scale data for the computational burden.
The study of the gender gap in the college wage premium has been a long-standing important question in labor economics. The literature has pointed out a stylized fact that the higher college wage premium for women as the major cause for attracting more women to attend and graduate from colleges than men (e.g., goldin2006homecoming, chiappori2009investment). Meanwhile, hubbard2011phantom pointed out a bias issue associated with the “topcoded” wage data, where the use of the sample from Current Population Survey (CPS) often censors the wage data at a maximum value. As a remedy for it, hubbard2011phantom proposed to use quantile regression that is robust to censoring. Analyzing the CPS data during 1970-2008, he found no gender difference in the college wage premium later years once the topcoded bias has been accounted for.
We revisit this problem by estimating and comparing the college wage premium between women and men. While the data also contains the topcoded issue, the quantile regression is less sensitive to the values of upper-tail wages. By applying the new S-subGD inference to the ultra-big dataset we are using, we aim at achieving the following goals in the empirical study: (1) \textcolor{black}{to} identify (if any) the heterogeneous effects across quantiles; (2) \textcolor{black}{to} understand the trends in the college wage premium respectively for female and male; and (3) \textcolor{black}{to} understand the gender difference in the college wage premium.
We use the data from IPUMS USA IPUMS that contains several millions of workers. The main motivation of using ultra-large datasets for wage regressions is to deal with high-dimensional regressors within a simple parametric quantile regression model. As has been pointed out in the literature, work experience is an important factor in wage regressions but is difficult to be measured precisely. \textcolor{black}{As career paths may be different across states and gender, one possible means of controlling for the experience is to flexibly interact workers' age, gender, and} state fixed effects, which would create over one thousand regressors\textcolor{black}{. Thus,} a very large sample size would be desirable to obtain precise estimates. But the ultra-large dataset makes most existing inference procedures for quantile regression fail to work.
We use the samples over six different years (1980, 1990, 2000-2015) from IPUMS USA at \url{https://usa.ipums.org/usa/}. In the years from 1980 to 2000, we use the 5% State sample which is a 1-in-20 national random sample of the population. In the remaining years, we use the American Community Survey (ACS) each year. The sampling ratio varies from 1-in-261 to 1-in-232 in 2001-2004, but it is set to a 1-in-100 national random sample of the population after 2005. To balance the sample size, we bunch the sample every 5 year after 2001. We also provide the estimation results using separate but smaller samples in 2005, 2010, and 2015 in the appendix, which is similar to those reported in this section.
We restrict our sample to $White$, $18 \le Age \le 65$, and $Wage \ge \$62$, which is a half of minimum wage earnings in 1980 ($\$3.10 \times 40 \mbox{hours} \times 1/2$). $Wage$ denotes the implied weekly wage that is computed by dividing yearly earnings by weeks worked last year. Since 2006, weeks worked last year are available only as an interval value and we use the midpoint of an interval. We only consider full-time workers who worked more than 30 hours per week. Then, we compute the real wage using the personal consumption expenditures price index (PCEPI) normalized in 1980. The data cleaning leaves us with 3.6-4.7 million observations besides 2001-2005, where we have around 2.5 million observations. Table (ref) reports some summary statistics on the key variables, where $Educ$ denotes an education dummy for some college or above. The table confirms that female college education has increased substantially over the years and female workers received more college education than male workers since 1990.
We use the following baseline model:
where $Wage_i$ is a real weekly wage in 1980 terms, $Female_i$ is a female dummy, $Educ_i$ is an education dummy for some college or above, and $X_i$ is a vector of additional control variables. For control variable $X_i$, we use 12 age group dummies with a four-year interval, 51 states dummies (including D.C.), and their interactions. Note that $(X_i\cdot Female_i)$ implies that there exist up to 3-way interactions. The model contains 1226 covariates in total. We also add 4 additional year dummies for the 5-year combined samples after 2001. We estimate the model using the proposed S-subGD method over 9 different quantiles: $\tau=0.1,0.2,\ldots,0.9$. We obtain the starting value of S-subGD inference by applying conquer to the 10% subsample.\footnote{Recall that $X_i$ contains a large number of dummy variables in this empirical application. If the subsample size is too small, conquer may face a singularity issue and cannot compute the initial value. Alternatively, one may use a randomazied algorithm of Yang:QR:ICML. } For the learning rate $\gamma_0 t^{-a}$, we set $a=0.501$ and $\gamma_0$ by the rule-of-thumb approach described in Section (ref). As an sensitivity analysis, we also vary $a$ over 17 equi-spaced points in $[0.501,0.661]$. The results are quite robust to the choice of $a$ as reported in the online appendix. We use the Compute Canada cluster system equipped with the AMD Milan CPUs whose clock speed is 2.65 GHz and with 650 GB of RAMs. For comparison, we also try to estimate the model using the standard QR method with the `quantreg' R package, but it fails to compute the result in the given cluster environments.
Figures (ref)--(ref) and Table (ref) summarize the estimation results. We also provide the full estimation results in the appendix.
Our results share several interesting features about the estimated college wage premium. First, some degrees of heterogeneity are found over different quantiles. At the upper tail quantile, $\tau=0.9$, male college premium is slightly higher than the female one since 2000, while female college premiums are higher over all years in other quantiles. This shows an interesting feature about high-income individuals that cannot be viewed by \textcolor{black}{mean regression}, though sensible economic interpretations might be subject to debates. Also, note that the 95% confidence intervals are slightly wider for tail quantiles ($\tau=0.1$ and 0.9, respectively).
Second, college premiums increase over time for both male and female. This result coincides with the overall findings in the literature. It is interesting that college premiums get flatter since 2010 for lower quantiles ($\tau \le 0.5$).
Third, when we focus on the median ($\tau=0.5$) as reported in Table (ref), we observe that the female-male college premium difference shows an inverse “U-shape” pattern. In addition, the difference is always significantly positive, which is different from the result in hubbard2011phantom. He found the gender difference to be insignificant starting from around the year 2000 from the CPS data. We also note that the computation time is within a reasonable range meanwhile it is not feasible to estimate the model with the alternative approaches.
We have proposed an inference method for large-scale quantile regression to analyze datasets whose sizes are of order $(n, p) \sim (10^7, 10^3)$. \textcolor{black}{Based on the stochastic subgradient descent updates, our method runs fast and constructs asymptotically pivotal statistics via random scaling.}
As an important extension, one can embed variable selections in the S-subGD framework. At the $i$ th update, define the t-statistic $ t_{i,j} = \bar\beta_{i,j}/ \sqrt{V_{i,jj}}. $ for $j=1,...,d. $ Then, variable $j$ is selected if $$ |t_{i,j}|>\lambda $$ for a tuning parameter $\lambda$. For instance, we can set $\lambda= 6.747$ which is the critical value for the 95% confidence level of the asymptotic distribution. For each variable $j$, we obtain a cumulative selection path along the iteration: $$ SP(j,i)= \frac{1}{i}\sum_{k=1}^i1\{|t_{k,j}|>\lambda\}. $$ The $SP$ path provides rich information about variable selections, which is particularly attractive when the signal-noise ratio is relatively weak. We would like to leave the theoretical justification of the selection path for future studies.
Besides the variable selection, there are a few more extensions worth pursuing in future research. First, we may build on chernozhukov2022fast to develop fast inference for the quantile regression process. Second, while we focus on the regular quantile regression models where the quantile is assumed to be bounded away from both zero and one, our framework is also potentially applicable to extreme quantile regression models chernozhukov2016extremal. The ability of handling ultra-large datasets is also appealing for extreme quantile regression models because the resulting sample size is considerably larger at the extreme quantiles. A formal treatment for these cases is out of the scope of this paper, and we leave it for future studies. Third, Shan:Yang:09 considered, among other things, an online setting for which time series data arrive sequentially and several quantile estimators are combined with weights being updated sequentially. It is an interesting topic to study the properties of our proposed method with time series data. Finally, Liu-et-al:ICML:2023 developed online inference on population quantiles with privacy concerns. It is another research topic to adapt our inference method for privacy-focused applications.