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.
54,225 characters · 11 sections · 30 citation commands
Least Squares Estimation using Sketched Data with Heteroskedastic Errors
\thispagestyle{empty} \setcounter{page}{0}
\onehalfspacing
Big data sets can be costly to store and analyze, and one approach around the data bottlenecks is to work with a randomly chosen subset, or a sketch, of the data. Data privacy may also dictate that a sketch of the data be made available for public use. The early works of sarlos-06, dmm:06 and dmms consider the algorithmic properties of the least squares estimator using sketched data. Subsequent work extends the analysis to ridge regression wang2017sketched,Liu2020Ridge, and logistic regression wang2019more. See, e.g., woodruff2014sketching, drineas2018lectures and martinsson2020randomized for a review. However, ma2015statistical, raskutti2016statistical, dobriban2019asymptotics and ma2020asymptotic have found that an optimal worse-case (algorithmic) error may not yield an optimal mean-squared (statistical) error. This has led to interest in better understanding the sketched least squares estimates in a Bayesian setting as in geppert-etal:17, or its asymptotic distribution as in ahfock-astle-richardson and ma2020asymptotic. LN:2020 highlights the tension between a large $m$ required for accurate inference, and a small $m$ for computation efficiency. To date, these results have been derived under the assumptions that the errors are homoskedastic and that the regressors are exogenous. But these assumptions are not innocuous. The estimates will be biased when the regressors are not exogenous, as would normally be the case in causal inference. And if the errors are heteroskedastic, test statistics must use standard errors robust to heteroskedasticity, or else inference will not be accurate even if the regressors are exogenous.
In this paper, we obtain the surprising result that when sketching is based on random projections, robust standard errors will not be needed, meaning that inference using the sketched estimates can proceed as though the errors were homoskedastic. The proof is obtained by analyzing the difference between the full sample and the sketched estimates in terms of degenerate $U$-statistics. However, the result does not hold when sketching is based on random sampling. Our analysis of the least squares estimator and two-stage least squares estimator shows that these findings hold both when the regressors are exogenous and endogenous. An implication is that inference, including first-stage F tests for instrument strength, may not require heteroskedasticity robust standard errors if the estimates are based on appropriately sketched data.
The following notation will be used. Let $\left\lVert a \right\rVert$ denote the Euclidean norm of any vector $a$. Let $A_{ij}$ or $[A]_{ij}$ denote the $(i,j)$ element of a matrix $A$. For $k = 1,\ldots,d$, let $\sigma_k(A)$ be a singular value of $A$. Let $\left\lVert A \right\rVert_2= \sigma_{\text{max}}(A)$ denote its spectral norm, where $\sigma_{\text{max}}(A)$, and $\sigma_{\text{min}}(A)$ are the largest and smallest singular values of $A$. The superscript $T$ denotes the transpose of a matrix. For an integer $n \geq 1$, $[n]$ is the set of positive integers from 1 to $n$. Let $\rightarrow_p$ and $\rightarrow_d$, respectively, denote convergence in probability and in distribution. For a sequence of random variables $A_n$ and a sequence of positive real numbers $a_n$, $A_n = o_p(a_n)$ iff $ a_n^{-1} A_n \rightarrow_p 0$; $A_n = O_p(a_n)$ iff $a_n^{-1} A_n$ is bounded in probability.
An accompanying R package is available on the Comprehensive R Archive Network (CRAN) at \url{https://CRAN.R-project.org/package=sketching} and all replication files are available at \url{https://github.com/sokbae/replication-LeeNg-2022-ICML}.
Given $n$ observations $\{ (y_i, X_i,Z_i): i=1,\ldots,n \}$, we consider a linear regression model:
where $y_i$ is the scalar dependent variable, $X_i$ is a $p \times 1$ vector of regressors, $\beta_0$ is a $p \times 1$ vector of unknown parameters. The innovation $e_i$ is said to be (conditionally) homoskedastic if $E[e_i^2|X_i]=E[e_i^2]$. Otherwise, $e_i$ is said to be heteroskedastic. The regressors are said to be exogenous if $E[e_iX_i]=0$. Otherwise it is endogenous. In that case, we assume a $q\times 1$ vector of instrumental variables, $Z_i$, satisfying $E[e_iZ_i]=0$ are available. In matrix form, the model given in (ref) can be written as
where $y$ and $e$ are $n \times 1$ vectors whose $i$-th rows are $y_i$ and $e_i$, respectively, and $X$ is the $n \times p$ matrix of regressors whose $i$-th row is $X_i^T$.
We first study the exogenous regressor case when $\mathbb E(e_iX_i)=0$. The least squares estimator $\widehat\beta_{OLS}:=(X^TX)^{-1} X^T y$ is $\sqrt{n}$ consistent and asymptotically normal, i.e., $ \sqrt{n}(\widehat\beta_{OLS}-\beta_0)\rightarrow_d N(0,V_1)$ as $n\rightarrow\infty$, where \[V_1 := [\mathbb E(X_iX_i^T) ]^{-1} \mathbb E(e_i^2X_iX_i^T) [\mathbb E(X_iX_i^T)]^{-1}\] is the heteroskedasticity-robust asymptotic variance. Under homoskedasticity, $V_1$ becomes \[ V_0:= \mathbb E(e_i^2) [\mathbb E(X_iX_i^T)]^{-1}.\] The point estimates $\widehat\beta$ can be used to test hypothesis, say, $H_0:\beta_2=\bar \beta_2$ using the $t$ test $\frac{\sqrt{n}(\widehat\beta_2-\bar \beta_2)}{\sqrt{[\widehat V]_{22}}}$, where $\widehat V$ is an estimate of either $V_1$ or $V_0$, $\beta_2$ is a specific element of $\beta_2$, $\bar \beta_2$ is the null value, and $[\widehat V]_{22}$ the (2,2) diagonal element of $\widehat V$. The distribution of this test under the null hypothesis crucially depends on the correct standard error $\sqrt{[\widehat V]_{22}}$ being used. Using $\widehat V_0$ when the robust estimator $\widehat V_1$ should have been used would lead to inaccurate inference, in the sense of rejecting the null hypothesis too often or not enough.
A sketch of the data $(y, X)$ is $(\widetilde{y}, \widetilde{X})$, where $\widetilde{{y}} = \Pi {y}$, $\widetilde{{X}} = \Pi {X}$, and $\Pi$ is usually an $m \times n$ random matrix. The sketched least squares estimator is $\widetilde\beta_{OLS}:=(\widetilde X^T\widetilde X)^{-1}\widetilde X^T \widetilde y$. Even though the sketched regression is based on a sample of size $m$, $\widetilde X^T\widetilde X = X^T \Pi^T \Pi X$ and $\widetilde X^T \widetilde y = X^T \Pi^T \Pi y$ can be seen as weighted moments in a sample of size $n$. Thus let $\widetilde{ g}_n := \widetilde{{X}}^T \widetilde{{e}}/n$, $\widehat g_n := {X}^T {e}/n$, $\widetilde{ A}_n := (\widetilde X^T\widetilde X/n)^{-1}$, and $\widehat{ A}_n := ( X^T X/n)^{-1}.$ Then
By the law of large numbers, $\widehat A_n - A = o_p (1)$, where $ A := [ \mathbb{E} ( X_i X_i^T ) ]^{-1}$, and by the central limit theorem, $\widehat g_n=O_p(n^{-1/2})$. We show in Section (ref) that for $\Pi$ with subspace embedding property
We study $\widetilde \beta_{OLS}$ under the following regularity conditions.
Assumptions (i) and (ii) are standard. For (iii), note that for a general random $\Pi$ whose $(k,i)$ element is $\Pi_{ki}$, the difference between the full and the sketched moments such as $\widetilde g_n-\widehat g_n$ and $\widetilde A_n-\widehat A_n$ are of the form
where $U \in \mathbb{R}^n$ and $V \in \mathbb{R}^n$ are vectors of certain i.i.d. random variables $(U_i, V_i) \in \mathbb{R}^2$ that are independent of $\Pi$,
and $\textrm{r.dim}(\Pi) \in \{m, n\}$ denotes the row dimension of $\Pi$.
There are two classes of sketching schemes to consider. Random sampling schemes have $\varphi_{ij} = 0$ for all $i \neq j$ because there is only one non-zero entry in each row of $\Pi$. In such cases, $T_{2n}$ is negligible and $T_{1n}$ is the leading term. The second class is random projection schemes with which $T_{1n}$ is asymptotically negligible and $T_{2n}$ is the leading term.
To gain intuition, we first provide results for Bernoulli sampling (BS) from the first type and countsketch (CS) from the second type.
Though Theorem (ref) indicates that both sampling schemes yield asymptotically normal estimates, their variances are different, and normality holds for different reasons. The proof is given in the Appendix. Here, we sketch the main arguments.
First, under BS, the sampling probability is determined by i.i.d. Bernoulli random variables with success probability $m/n$. Thus, $\Pi = \sqrt{\frac{n}{m}} B$ is an $n\times n$ matrix (not $m\times n$), where $B$ is a diagonal sampling matrix. We have
Since the summands are i.i.d. with mean zero with variance $( \frac{n}{m} - 1 ) \mathbb{E} ( e_i^2 X_i X_i^T)$, applying central limit theorem for i.i.d. observations yields the sandwich variance $V_1$.
Consider now CS. Each column of its $\Pi$ has one non-zero entry taking on value $\{ +1, -1\}$ randomly drawn with equal probability and located uniformly at random. For such $\Pi$ and every nonzero $c \in \mathbb{R}^p$,
where $\bar{X}_i(c) := \sum_{u}^q c_u X_{iu}$ is a weighted sum of elements of the $i$-th row of $X$. The term $T_{n1}$ is identically zero because there is only one non-zero entry per column of $\Pi$.
To analyze $T_{n2}$, let $W_i = (Y_i, X_i^T, \Pi_{1i},\ldots, \Pi_{mi})^T$. Since the columns of $\Pi$ are i.i.d., $\{W_i : i=1,\ldots,n \}$ are i.i.d. Now let $w = (y, x^T, \pi_1, \ldots, \pi_m)^T$ be a non-random index. Define $ \widetilde{H} (w_1, w_2) : = \sum_{k=1}^{m} \bar{x}_{1}(c) \pi_{k1} S\pi_{k2} e_{i} $ and $ H (w_1, w_2):= \widetilde{H} (w_1, w_2) + \widetilde{H} (w_2, w_1). $ We can write
noting that $H (w_1, w_2) = H (w_2, w_1)$, and $\mathbb{E} ( H (W_1, W_2) | W_1) = \mathbb{E} ( H (W_1, W_2) | W_2) = 0$. Importantly, $T_{n2} $ has now been represented as a degenerate $U$-statistic. In general, the asymptotic distribution of such statistics is either a weighted average of independent, centered chi-square random variables with complex weights, or a centered normal distribution. But if the conditions given in hall1984central are satisfied, the latter holds. Precisely,
A sufficient condition for this result which we verify in the Appendix is
where $ G (w_1, w_2) := \mathbb{E} [ H (W_1, w_1) H (W_1, w_2) ]$. Furthermore, we also verify that for $W_i \neq W_j$, $$ \frac{1}{2} \mathbb{E} [ H^2 (W_i, W_j) ] = \frac{1}{m} \mathbb{E} ( \bar{X}^2_i(c)) \mathbb{E} (e_i^2 ). $$
Note that $\mathbb E (e_i^2)$ appears separately from $\mathbb E(\bar X_i^2(c))$. This is key to the claim in Theorem (ref) that when $\widetilde \beta$ is based on CS, $m^{1/2}(\widetilde\beta-\widehat\beta)\rightarrow_d N(0,E[e_i^2] A)=N(0,V_0) $. Analogous arguments show that each entry of $(\widetilde X^T \widetilde X- X^T X)$ can also be written as a degenerate $U$-statistic and $\widetilde A_n-\widehat A=o_p(1)$, which is needed for consistent estimation of $V_0$ and $V_1$.
As discussed in charikar2004finding and clarkson2013low,clarkson2017low, the sparsity of $\Pi$ significantly reduces the run time required of the countsketch to compute $\Pi A$ to O(nnz(A)), where nnz(A) is the number of non-zero entries of $A$. Another appeal of countsketch is that the sketches can be obtained by streaming without constructing $\Pi$. Here, we show that countsketch removes heteroskedasticity which is appealing because it simplifies inference. In the next section, we study the mean-squared sketching error and show that part (i) of Theorem (ref) also holds for other $\Pi$s in the first class, while part (ii) holds for other $\Pi$s in the second class. Section (ref) then shows that these results also hold when the regressors are not exogenous.
For a random variable $G$, let $\mathrm{MSE} ( G ) = [\mathbb{E} (G)]^2 + \mathrm{Var} (G)$ denote the mean squared error. We now analyze the asymptotic behavior of mean squared sketching errors of $(U^T \Pi^T \Pi V - U^T V)/n$, where $U \in \mathbb{R}^n$ and $V \in \mathbb{R}^n$ denote vectors of i.i.d. random variables $(U_i, V_i) \in \mathbb{R}^2$ that are independent of $\Pi$, with $\mathbb{E} (U_i^4) < \infty$ and $\mathbb{E} (V_i^4) < \infty$. Recall that $m = m_n \rightarrow \infty$ but $m/n \rightarrow 0$ as $n \rightarrow \infty$.
For sketching by random sampling with replacement (RS), we suppose that for each $t = 1,\ldots,m$, we sample $k_t$ from $[n]$ with probability $p_i := \mathrm{Pr}(k_t = i)$ independently and with replacement. The random matrix $\Pi \in \mathbb{R}^{m \times n}$ is then
where $\iota_k$ is the $k$-th column vector of the $n \times n$ dimensional identity matrix. Sketching schemes of the RS class have properties characterized by Lemma (ref) in the Appendix. Importantly, each $\Pi$ in this class has $T_{2n}=0$ and as a result, $T_{n1}$ is the only term we need to consider. An important example in the RS class is uniform sampling with replacement with $p_i = n^{-1}$.
The mean-squared errors for RS and BS are the same if $U_i V_i$ is mean zero.
Theorem (ref) is useful in two ways. First, let $U = c^T X^T$ and $V = X c$. Under Assumption (ref) and for all nonzero $c \in \mathbb{R}^p$, Theorem (ref) yields
By Chebyshev's Inequality, $( \widetilde X^T\widetilde X - X^T X ) /n = O_p ( m^{-1/2} )$. Similarly, for $U=X c$ and $V=e$, the theorem implies $\widetilde{ g}_n - \widehat g_n = O_p ( m^{-1/2} )$. Applying continuous mapping theorem gives the result stated earlier that $\widetilde{\beta}_{OLS} - \widehat{\beta}_{OLS} = A ( \widetilde{ g}_n - \widehat g_n) + o_p ( m^{-1/2} )$.
It is known that the efficient estimator under heteroskedasticity is the generalized least squares (GLS), defined as \[\widehat\beta_{GLS}=(X^T\Omega^{-1} X)^{-1} X^T\Omega^{-1} y,\] where $\Omega$ is $n\times n$ a diagonal matrix with $\Omega_{ii}=\sigma^2_i$. GLS weights each observation with $\sigma_i^{-1}$ so that the errors in the weighted regression are homoskedastic. Now the OLS estimator applied to sketched data can be written as \[ \widetilde\beta_{OLS}= (X^T\Pi^T\Pi X)^{-1} X^T\Pi^T \Pi y.\] A question of interest is whether $\Pi^T\Pi$ can play the role of $\Omega^{-1}$. The theorem sheds light on this problem as its second use is to obtain the asymptotic variance of the sketched estimator. To this end, let again $U_i=c^T X_i$ and $V_i=e_i$. Assuming $\mathbb E(g_i)=0$ where $g_i=e_i X_i$, RS and BS imply:
The asymptotic standard error is generally the expectation of a product of $e_i^2$ and $(c^TX_i)^2$ and becomes the product of two expectations only under homoskedasticity when $\mathbb E[e_i^2 c^T X_iX_i^T c]=\mathbb E[e_i^2] \mathbb E[c^T X_iX_i^Tc]$. Thus, under RS and BS, the asymptotic variance of $\widetilde\beta_{OLS}$ is $V_0$ only if homoskedasticity is explicitly imposed, implying that $\Pi^T\Pi$ corresponding to random sampling will not homogenize error variance in the same way that $\Omega^{-1}$ can.
It is noteworthy that even under homoskedsaticity, we cannot always use a central limit theorem for i.i.d. data even if the full sample of data are i.i.d. because the sampling scheme may induce dependence in the sketched data. Thus the asymptotic normality result can only be analyzed on a case by case basis. ma2020asymptotic confronts a similar problem when studying the asymptotic distribution of estimators in linear regressions under random sampling with replacement and homoskedastic errors. Let $K_{i}$ and $p_i$, respectively, denote the number of times and the probability that $i^{t h}$ observation is sampled. Their estimator has $W =\operatorname{diag}\left\{K_{i} / ( m p_i ) \right\}_{i=1}^{n}$ playing the role of $\Pi^T\Pi$. Our Theorem 3.1 applies to their setup with uniform sampling where $p_i=1/n, n\ge m$, but it would not apply when $p_i$ is data dependent. In this case, ma2020asymptotic also cannot use a central limit theorem for i.i.d. data. Instead, they apply Hay\'ek-Sidak central limit theorem and use Poissonization to account for dependence in the sketched data that arises from sampling.
Sketching schemes in the RP class have properties characterized by Lemma (ref) in the Appendix if $\Pi \in \mathbb{R}^{m \times n}$ is a random matrix with the following properties:
Under Assumption (ref), $T_{n1} = O_p( n^{-1/2})$ is asymptotically negligible, and $T_{n2}$ becomes the leading term for RP. As discussed above, the $\Pi$ for the CS only has one non-zero entry in each column. Since $\Pi_{ki} \Pi_{kj} = 0$ for all $k, i \neq j$, it is straightforward to check that the above conditions are satisfied. For Gaussian random projections, \[ \text{GP}: \Pi_{ki}\sim N(0,m^{-1}).\] Since all elements of $\Pi$ are i.i.d. with mean zero, variance $m^{-1}$ and the fourth moment $O(m^{-1})$, the conditions are also satisfied. The SRHT has \[ \text{SRHT}: \Pi = \sqrt{\frac{n}{m} } S H D,\] where $S \in \mathbb{R}^{m \times n}$ is a uniform sampling matrix with replacement, $H \in \mathbb{R}^{n \times n}$ is a normalized Walsh-Hadamard transform matrix, and $D \in \mathbb{R}^{n \times n}$ is a diagonal Rademacher matrix with i.i.d. entries of $\pm 1$. The Appendix shows that the conditions for RP hold for SRHT.
The following theorem gives the asymptotic mean squared sketching errors of RP schemes.
The limiting MSE of {RP} is simply the product between two marginal expectations when $\mathbb{E} ( U_i V_i ) = 0$ (and not the expectation of the product). It implies \[\mathrm{MSE} \left[ \sqrt{m} \, c^T( \widetilde{ g}_n -\widehat g_n) \right]= \mathbb{E} ( e_i^2 ) \mathbb{E} ( c^T X_i X_i^T c )\] and is the reason why the asymptotic variance for $\widetilde\beta_{OLS}$ for RP schemes is $V_0$.
If $e_i^2$ and $(c^T X_i)^2$ are positively (respectively, negatively) correlated, the limiting MSE of {RP} is smaller (respectively, larger) than that of {RS} and {BS}. The limiting MSE is the same if $e_i^2$ and $(c^T X_i)^2$ are uncorrelated.
Asymptotic normality of $\widetilde\beta$ can be established by applying a central limit theorem for degenerate $U$-statistic if the columns of $\Pi$ are i.i.d., as reported in part (ii) of Theorem (ref). The SRHT and SRFT are not covered by this result because the columns of their $\Pi$ matrix are not i.i.d. and requires a limit theorem for a particular type of mixing data. In general, establishing asymptotic normality of $\widetilde\beta$ based on SRHT or SRFT require different proof techniques. The approach taken in ahfock-astle-richardson is to condition on the data $\mathcal{D}_n$ and apply a central limit theorem for a triangular array of random variables. We do not condition on the data and appeal to the theory of degenerate $U$-statistics. Though deriving distribution theory for the SRHT and SRFT estimates is not straightforward, we will show in simulations that their finite sample properties are similar to those of CS.
The 2SLS estimator is appropriate when $\mathbb E(X_ie_i) \ne 0$ but exongeous instruments $Z_i$ satisfying $\mathbb E(Z_ie_i)=0$ are available. The 2SLS estimator is \[\widehat\beta_{2SLS}=(X^TP_Z X)^{-1} X^T P_Z y,\] where $P_Z:=Z(Z^T Z)^{-1}Z^T$ is the projection matrix. The estimator first projects on $Z$ to purge the variations in $X$ correlated with $e$, and in the second step replaces $X$ with $\widehat X=P_Z X$. Let $\widehat{ g}_n := {{Z}}^T {{e}}/n$ and $\widehat{ A}_n := [( X^T Z/n) ( Z^T Z/n)^{-1} ( Z^T X/n)]^{-1} ( X^T Z/n) ( Z^T Z/n)^{-1}$. Analyzing $\widehat\beta_{2SLS}-\beta_0=\widehat A_n \widehat g_n$ under Assumption (ref) given below, as $n \rightarrow \infty$, we have
where $W_1 := A \,\mathbb{E} (e_i^2 Z_i Z_i^T ) \,A^T$ with
Under homoskedasticity, $\mathbb E(e_i^2|Z_i)=\sigma^2$ and $W_1$ reduces to \[ W_0 := \mathbb{E} ( e_i^2 ) [ \mathbb{E} ( X_i Z_i^T ) [ \mathbb{E} ( Z_i Z_i^T ) ]^{-1} \mathbb{E} ( Z_i X_i^T ) ]^{-1}. \] A sketched version of the 2SLS estimator is
We now provide some algorithmic results not previously documented in the literature.
Low level conditions for Assumption (ref)(i)-(iii) are given in cohen_et_al:ICALP2016, among others. Assumption (ref)(i) is equivalent to the statement that the all eigenvalues of $U_Z^T \Pi^T \Pi U_Z$ are bounded between $[1 - \varepsilon_1, 1 + \varepsilon_1]$. This ensures that $\widetilde{{Z}}^T \widetilde{{Z}}$ is non-singular with probability at least $1-\delta$. Part (iv) strengthens non-singularity of $X^T P_Z X$ to require that $\sigma_{\text{min}}^2 ( U_Z^T U_X )$ is strictly positive and bounded below by the constant $2 f_1 (\varepsilon_1, \varepsilon_2)$.
The sketched estimator $\widetilde\beta_{2SLS}$ involves, firstly, a regression of $\widetilde X$ on $\widetilde Z$, and then a regression of $\widetilde y$ on the fitted values in the first step. The estimator thus depends on adequacy of subspace approximation in both steps. Theorem (ref) provides a worst-case bound for $\widetilde\beta_{2SLS}-\widehat\beta_{2SLS}$ with the data $\mathcal{D}_n$ being fixed. It depends on (i) $\varepsilon_j, j=1,2,3$, (ii) variability of $\left\| \widehat{e} \right\|$, (iii) the signal from $X$ as given by $\sigma_{\text{min}}(X)$, and (iv) instrument strength as given by $\sigma_{\text{min}} ( U_Z^T U_X )$. The sketched estimator can be arbitrarily close to the full sample estimate with high probability, provided that the subsample size $m$ is sufficiently large, $X$ is linearly independent, and the instrument $Z$ is sufficiently relevant for $X$.
Though 2SLS is a two step estimator, we can still write
as in the OLS case, but now $\widetilde{ g}_n := \widetilde{{Z}}^T \widetilde{{e}}/n$, and $$ \widetilde{ A}_n := [(\widetilde X^T \widetilde Z/n) (\widetilde Z^T\widetilde Z/n)^{-1} (\widetilde Z^T\widetilde X/n)]^{-1} (\widetilde X^T \widetilde Z/n) (\widetilde Z^T \widetilde Z/n)^{-1}. $$ A statistical analysis of $\widetilde \beta_{2SLS}$ requires additional assumptions.
Arguments similar to those used to prove Theorem (ref) lead to the following.
Theorem (ref) provides statistical properties of the sketched 2SLS estimator in Theorem (ref) to complement the algorithmic results.
Theorem (ref) states that when the data are sketched by RP, $\widetilde{\beta}$ is asymptotically normally distributed with mean $\widehat{\beta}$ and variance $W_0/m$. Under our assumptions, $W_0$ can be consistently estimated by
where $\widehat{e} := {y} - {X} \widehat{\beta}$ (not the residuals from the second step).
Interestingly, the asymptotic variance $W_0$ is the same as if the errors in the full sample regression were homoskedastic. But the result follows from estimation using sketched data rather than by assumption. This is not the case when inference is based on the full sample estimates, or the estimates computed from sketched data of the RS type. In such cases, a homoskedastic covariance weighting matrix would be inefficient since $\mathbb{E} ( e_i^2 | Z_i ) \neq \mathbb{E} ( e_i^2 )$.
Our analysis can be extended to the two-sample 2SLS estimator analyzed in AK:92,AK:95 and Inoue:Solon:10. However, it is not pursued for brevity of the paper.
In the econometrics literature, the instruments are said to be relevant if $\mathbb E[Z_iX_i^T]\ne 0$. The latter is formalized by the rank condition in Assumption (ref)(ii). Tests for instrument relevance usually require robust standard errors corresponding to the parameter estimates in a regression of $X$ on $Z$ unless heteroskedasticity can be ruled out. An implication of our preceding analysis is that this is not necessary when the regression is estimated on data sketched by RP, as will be illustrated below.
In applications, researchers would like to test a hypothesis about $\beta_0$ using a sketched estimate, and our results provide all the quantities required for inference. In the exogenous regressor case, we generically have
where the form of $\widetilde V_{m}$ depends on $\Pi$. For any $\Pi$ in BS or RP class, we can use white1980's heteroskedasticity-consistent estimator: \[\widetilde V_m=\widetilde V_{1,m}= (\widetilde X^T\widetilde X)^{-1} (\sum_{i=1}^m \widetilde X_i\widetilde X_i^T \widetilde e_i^2) (\widetilde X^T\widetilde X)^{-1}.\] For $\Pi$ in the RP class, we can let $\widetilde s^2_{OLS}:=\frac{1}{m}\sum_{i=1}^m (\widetilde y_i-\widetilde X_i^T \widetilde\beta_{OLS})^2$. Then without assuming homoskedasticity, \[\widetilde V_m=\widetilde V_{0,m}= \widetilde s_{OLS}^2(\widetilde X^T\widetilde X)^{-1},\]
In the endogenous regressor case, $\widetilde W_m^{-1/2}(\widetilde \beta_{2SLS}-\beta_0)\approx N(0,I_p)$. For RP, we let $\widetilde s^2_{2SLS}:=\frac{1}{m}\sum_{i=1}^m ( \widetilde{y}_i - \widetilde{X}_i^T \widetilde{\beta}_{2SLS})^2$ and $ \widetilde{W}_{0,m} := \widetilde s^2_{2SLS} \left( \widetilde{X}^T \widetilde{Z} ( \widetilde{Z}^T \widetilde{Z} )^{-1} \widetilde{Z}^T \widetilde{X} \right)^{-1}. $ For BS, we define $\widetilde W_{1m}$ from $\widetilde W_1$.
Sketching estimators require a choice of $m$. From the algorithmic perspective, $m$ needs to be chosen as small as possible to achieve computational efficiency. If $\Pi$ is constructed from SRHT, the size of $m$ is roughly (ignoring the log factors) of order $q$ in the best case. The requirement for countsketch is more stringent and is proved in the appendix (see Theorem (ref)). In view of this, we may set
where $C_m$ is a constant that needs to be chosen by a researcher. However, statistical analysis often cares about the variability of the estimates in repeated sampling and a larger $m$ may be desirable from the perspective of statistical efficiency. An inference-conscious guide $m_2$ can be obtained in LN:2020 by targeting the power at $\bar\gamma$ of a one-sided $t$-test for given nominal size $\bar\alpha$. In particular, let $\widetilde\beta$ be $\widetilde\beta_{OLS}$ if $\mathbb E[X_i e_i]=0$ and let $\widetilde\beta$ be $\widetilde \beta_{2SLS}$ when $\mathbb E[X_i e_i]\ne 0$ but $\mathbb E[Z_i e_i]=0$. For pre-specified effect size $c^T(\beta^0-\beta_0)$,
where $S(\alpha,\gamma):=\Phi^{-1}(\gamma) +\Phi^{-1}(1-\alpha)$ and $\textsc{se}(c^T \widetilde\beta)$ is the standard error of $c^T \widetilde\beta$.
Alternatively, a data-oblivious sketch size for a pre-specified $\tau_2(\infty)$ is defined as
Note that $m_3$ only requires the choice of $\bar \alpha, \bar\gamma,$ and $\tau_2(\infty)$ which, unlike $m_2$, can be computed without a preliminary sketch. The condition $m / n \rightarrow 0$ can be viewed as $\tau_2 (\infty) \rightarrow \infty$ as $n \rightarrow \infty$.
In this section, we use Monte Carlo experiments to establish that when the errors are homoskedastic, estimates based on data sketched by random sampling or random projections will yield accurate inference. However, when the errors are heteroskedastic, sketching by random sampling will yield tests with size distortions, rejecting with much higher probability than the nominal size, unless robust standard errors are used.
We first consider the simulation design for which all the regressors are exogenous. The regressors $X_i = (1, X_{2,i}, \ldots, X_{p,i})^T$ consist of a constant term and a $(p-1)$-dimensional random vector $(X_{2,i}, \ldots, X_{p,i})^T$ generated from a multivariate normal distribution with mean zero vector and the variance covariance matrix $\Sigma$, whose $(i,j)$ component is $\Sigma_{ij} = \rho^{|i-j|}$ with $\rho = 0.5$. The dependent variable is generated by \[ y_i = X_i^T \beta_0 + \sigma (X_i) e_i, \] where $\beta_0 = (0, 1, \ldots, 1)^T$, and $e_i$ is generated from $N(0,1)$ independently from $X_i$. We consider two designs for $\sigma(X_i)$: (i) homoskedastic design $\sigma(X_i) = 1$ for all $i$ and (ii) heteroskedastic design $\sigma(X_i) = \exp ( X_{p,i} )$, where $X_{p,i}$ is the $p$-th element of $X_i$. Throughout the Monte Carlo experiment, we set $n = 10^6$, $m = 500$, and $p = 6$. There were 5,000 replications for each experiment. Six sketching methods are considered: (i) Bernoulli sampling, (ii) uniform sampling, (iii) leverage score sampling and reweighted regression as in ma2020asymptotic; (iv) countsketch, (v) SRHT, (vi) subsampled randomized Fourier transforms using the real part of fast discrete Fourier transform (SRFT). Table (ref) reports the empirical size and power of the $t$-test. The null and alternative hypotheses are that $H_0: c^T \beta_0 = 1$ vs. $H_1: c^T \beta_0 \neq 1$ with $c^T = (0,\ldots, 0, 1)$. Equivalently, the null hypothesis is $\beta_p=1$. The power is obtained for the null value $c^T \beta_0 = 1.1$ for the homoskedastic design and $c^T \beta_0 = 1.4$ for the heteroskedastic design, respectively. The nominal size is 0.05.
In column (1) in Table (ref), we report the size of the test, namely, the probability of rejecting $H_0$ when the null value is true. In this column, the $t$-statistic is constructed using homoskedasticity-only standard errors S.E.0. Though many methods perform well, both Bernoulli and uniform sampling show substantial size distortions for the heteroskedastic design. Leverage score sampling combined with reweighted regression seems to account for heteroskedasticity to some extent, but not enough to remove all size distortions. In column (2) which reports results using robust standard errors S.E.1, all methods have satisfactory size. In column (3), we report the power of the test, i.e., the probability of rejecting $H_0$ when the null value is false. For the heteroskedastic design, the powers of the tests using homoskedastic standard errors S.E.0 are worse for Bernoulli, uniform and leverage samplings than those for countsketch, SRHT, and SRFT. The power loss of the RS schemes is much more pronounced when the robust standard errors are used in column (4). This efficiency loss is consistent with asymptotic theory developed in the paper because the squared regression error is positively correlated with one of the elements of squared $X_i$ under the heteroskedastic design. All RP schemes perform similarly, hinting that even though a formal proof awaits future research, asymptotic normality may also hold for both SRHT and SRFT, and not just countsketch.
We now move to the case when the regressors are $X_i=(1,X_{2,i}, \ldots, X_{p-1,i}, X_{p,i})^T$, and $y_i$ is generated by
where $\epsilon_i\sim N(0,1)$ is randomly drawn independently from $X_i$ and $\eta_i$. The first $p-1$ regressors, including the intercept term, are exogenous, but
where $\zeta_0=(\zeta_{1,0}, \ldots, \zeta_{q,0})^T$, $\eta_{i}\sim N(0,1)$ independently from $Z_i=(1, Z_{2,i}, \ldots, Z_{q,i})^T$. The presence of $\eta_i$ in both (ref) and (ref) induces endogeneity of $X_{p,i}$.
In each of the 1000 replications, $(X_{2,i},\ldots, X_{p-1,i})^T=(Z_{2,i},\ldots, Z_{p-1,i})^T$, while the $(q-1)$-dimensional $Z_i$ is multivariate normal with mean zero and the variance $\Sigma$, whose $(i,j)$ component is $\Sigma_{ij} = \rho^{|i-j|}$ with $\rho = 0.5$. We consider two designs for $\sigma_1(Z_i)$: (i) homoskedastic design $\sigma_1(Z_i) = 1$ for all $i$ and (ii) heteroskedastic design $\sigma_1(Z_i) = \exp \left( \frac{5}{q} \sum_{j=2}^{q} |Z_{j,i}| \right) /100$. As in the previous section, we set $n = 10^6$, $m = 500$, $p = 6$, and $q=21$. We consider five sketching schemes and no longer include leverage score sampling since it is unclear how to implement it in the case of 2SLS. The nominal size is 0.05. Throughout, $(\zeta_{1,0},\ldots,\zeta_{p-1,0})= (0,0.1,\ldots,0.1)^T$, but values of $\zeta_{j,0}$ for $j \geq p$ depend on the context as explained below.
We first examine the so-called first-stage F-test for instrument relevance. In this case of a scalar endogenous regressor, the null hypothesis of irrelevant instruments amounts to a joint test of $H_0: \zeta_{j,0} =0$ for every $j=p,\ldots, q$ in (ref). The size of the test is evaluated at $\zeta_{j,0}= 0$ and the power at $\zeta_{j,0}=0.1$ for $j=p,\ldots, q$. The F-test statistic is constructed as $$ F = \frac{1}{q-p+1} \widehat{\zeta}_{-(p-1)}^T \left( [\widehat V]_{-(p-1),-(p-1)} \right)^{-1} \widehat{\zeta}_{-(p-1)}, $$ where $\widehat{\zeta}_{-(p-1)}$ is a $(q-p+1)$-dimensional vector of the OLS estimate $\widehat{\zeta}$ of regressing $X_{p,i}$ on $Z_i$, excluding the first $(p-1)$ elements, and $[\widehat V]_{-(p-1),-(p-1)}$ is the corresponding submatrix of $\widehat V$.
In Table (ref), we report the size and power of the F-test for $H_0: \zeta_{p,0} = \zeta_{p+1,0} = \ldots = \zeta_{q,0} = 0$ using homoskedasticity-only (V.0) and heteroskedasticity-consistent (V.1) asymptotic variance estimates, respectively. As in the previous subsection, Bernoulli and uniform sampling sketches suffer from size distortions in the heteroskedastic design but V.0 is used. Tests based on V.1 have good size without sacrificing much power when the $F$ test is constructed from data sketched by RP.
We now turn to 2SLS estimation of $\beta_0$. To ensure that the instruments are powerful enough to estimate $\beta_0$ well, we now set $\zeta_{j,0}=0.5$ for $j=p,\ldots, q$ with $\sigma_1 (Z_i) = 1$ for all $i$. We set $\beta_0 = (0, 1, \ldots, 1)^T$ and consider two designs for $\sigma_2(Z_i)$: (i) homoskedastic design $\sigma_2(Z_i) = 1$ for all $i$ and (ii) heteroskedastic design $\sigma_2(Z_i) = \exp \left( \frac{5}{q} \sum_{j=2}^{q} |Z_{j,i}| \right) /100$.
As in the previous subsection, we test $H_0: \beta_p=1$ against $H_1: \beta_p\ne 1$, or equivalently, $c^T = (0,\ldots, 0, 1)$. The power is obtained for $ \beta_p = 1.05$ in the homoskedastic design and $ \beta_p = 1.10$ for the heteroskedastic design, respectively. Table (ref) reports results for nominal size of 0.05. Basically, the same patterns are observed as in the previous section. Thus, simulations support the theoretical result that robust standard errors are not needed for inference when estimation is based on sketched data using sketching schemes in the RP class.
An exemplary application of the 2SLS in economics is causal inference, such as to estimate the return to education. Suppose that $y_i$ is the wages for worker $i$ (typically in logs) and $X_i$ contains educational attainment $\texttt{edu}_i$ (say, years of schooling completed). Here, the unobserved random variable $e_i$ includes worker $i$'s unobserved ability among other things. Then, $\texttt{edu}_i$ will be correlated with $e_i$ if workers with higher ability tends to attain higher levels of education. The least-squares estimator may not provide a consistent estimate of the return to schooling. To overcome this problem, economists use an instrumental variable that is uncorrelated with $e_i$ but correlated with $\texttt{edu}_i$.
We re-examine the OLS and 2SLS estimates of return to education in columns (1) and (2) of Table IV in AK1991. The dependent variable $y$ is the log weekly wages, the covariates $X$ include years of education, the intercept term and 9 year-of-birth dummies $(p=11)$. Following AK1991, the instruments $Z$ are the exogenous regressors (i.e., the intercept and year-of-birth dummies) and a full set of quarter-of-birth (one quarter omitted) times year-of-birth interactions $(q= 1 + 9 + 3 \times 10 = 40)$. Their idea was that season of birth is unlikely to be correlated with workers' ability but can affect educational attainment because of compulsory schooling laws. The full sample size is $n = 247,199$.
To construct sketched data, we need to choose the sketch size $m$. We use a data-oblivious sketch size $m_3$ defined in (ref) with target size set to $\alpha= 0.05$ and target power to $\gamma = 0.8$, giving $S^2(\bar\alpha,\bar\gamma) = 6.18$. It remains to specify $\tau_2(\infty)$, which can be interpreted as the value of $t$-statistic when the sample size is really large.
In the OLS case, we take $\tau_2(\infty) = 10$ resulting in $m = 15,283$ (about 6% of $n$). Table (ref) reports empirical results for the OLS estimates. For each sketching scheme, only one random sketch is drawn; hence, the results can change if we redraw sketches. Remarkably, all sketched estimates are 0.08, reproducing the full sample estimate up to the second digit. The sketched homoskedasticity-only standard errors are also very much the same across different methods. The Eicker-Huber-White standard error S.E.1 is a bit larger than the homoskedastic standard error S.E.0 with the full sample. As expected, the same pattern is observed for Bernoulli and uniform sampling, as these sampling schemes preserve conditional heteroskedasticity.
For 2SLS, as it is more demanding to achieve good precision, we take $\tau_2(\infty) = 5$, resulting in $m = 61,132$ (about 25% of $n$). Table (ref) reports empirical results for the 2SLS estimates. The sketched estimates vary from 0.053 to 0.115, reflecting that the 2SLS estimates are less precisely estimated than the OLS estimates. Both types of standard errors are almost identical across all sketches for 2SLS, suggesting that heteroskedasticity is not an issue in this data.