Heteroscedasticity-Robust Inference in Linear Regression Models With Many Covariates

Abstract

We consider inference in linear regression models that is robust to heteroscedasticity and the presence of many control variables. When the number of control variables increases at the same rate as the sample size the usual heteroscedasticity-robust estimators of the covariance matrix are inconsistent. Hence, tests based on these estimators are size distorted even in large samples. An alternative covariance-matrix estimator for such a setting is presented that complements recent work by Cattaneo, Jansson, and Newey. We provide high-level conditions for our approach to deliver (asymptotically) size-correct inference as well as more primitive conditions for three special cases. Simulation results and an empirical illustration to inference on the union premium are also provided. Supplementary materials for this article are available online.

Keywords:

• Heteroscedasticity
• Inference
• Many regressors
• Statistical leverage

1 Introduction

When performing inference in linear regression models it is common practice to safeguard against (conditional) heteroscedasticity of unknown form. The estimator of the covariance matrix of the least-squares estimator proposed by Eicker (1963, 1967) and White (1980) is known to be biased. The bias can be severe if the regressor design contains observations with high leverage (Chesher and Jewitt 1987). ¹ A necessary condition for the least-squares estimator to be asymptotically normal is that maximal leverage vanishes as the sample size grows large (Huber 1973). This condition, then, also implies consistency of the robust covariance-matrix estimator (under regularity conditions).

The requirement that maximal leverage vanishes is problematic when the regressors include a large set of control variables. Under asymptotics where their number, q_n , grows with the sample size, n, the robust covariance-matrix estimator will be inconsistent unless $q_n/n\to 0$, as demonstrated by Cattaneo, Jansson, and Newey (2018). They obtained the same result for the other members of the so-called HC-class of covariance-matrix estimators (see, e.g., Long and Ervin 2000; MacKinnon 2012 for reviews) and showed that the jackknife variance estimator of MacKinnon and White (1985), although inconsistent, can be used to perform asymptotically conservative inference under asymptotics where $q_n/n\nrightarrow 0$. On the other hand, a bias-corrected covariance-matrix estimator in the spirit of Hartley, Rao, and Kiefer (1969) and Bera, Suprayitno, and Premaratne (2002) is shown to be consistent under conditions that bound maximal leverage below $\frac{1}{2}$. One implication of these conditions is that $\lim \underset{n}{\sup }{q}_n/n<\frac{1}{2}$. Although not guaranteed to be positive semidefinite (see Bera, Suprayitno, and Premaratne 2002 for a discussion), this estimator is attractive as it is asymptotically equivalent to a minimum-norm unbiased estimator in the sense of Rao (1970).

In this article, we discuss an alternative estimator of the covariance matrix that can deal with designs where maximal leverage is bounded away from 1. As such it remains consistent when $\lim \underset{n}{\sup }{q}_n/n<1$. To show this, we need to impose additional conditions relative to Cattaneo, Jansson, and Newey (2018). A consistency result is first provided under high-level conditions. Primitive conditions are then given for three special cases; the partially linear regression model, the one-way model for short panel data, and the generic linear model with increasing dimension. Again, our covariance-matrix estimator need not be positive semidefinite in small samples. It is further not invariant to changes in the scale of the regression slopes. Achieving such invariance under asymptotics where $\lim \underset{n}{\sup }{q}_n/n>\frac{1}{2}$ appears difficult, as we discuss below.

The idea underlying our variance estimator can be traced back to work by Kline, Saggio, and Sølvsten (2019, Remark 4 and Lemma 5); see below. The chief difference lies in the conditions under which the consistency result is obtained. They considered settings where the observations are independent and the regressors are fixed, and entertained models that are correctly specified with regression functions that are uniformly bounded. This is reasonable in analysis-of-variance problems, which are the main focus of application in their work. Here, we maintain the framework of Cattaneo, Jansson, and Newey (2018). Some dependence between observations is allowed, regressors are stochastic, and the model can feature (vanishing) misspecification bias. Our consistency result can be understood to be an extension of Kline, Saggio, and Sølvsten (2019) to such settings. It allows for regressors to have unbounded support and for observations to depend on a growing number of parameters. Our conditions do rule out dynamic models, for example. Verdier (2020) provides estimation and inference results for two-way models estimated from short panel data that allow for dynamics over time.

In Section 2, we introduce the framework, present our covariance-matrix estimator, and provide a consistency result under a set of high-level conditions. In Section 3, we connect our work to the literature and, notably, to Cattaneo, Jansson, and Newey (2018) and Kline, Saggio, and Sølvsten (2019). In Section 4, we provide primitive conditions for three special cases of our setup. In Section 5, we present and discuss the results from Monte Carlo experiments and apply our variance estimator to perform inference on the union wage premium. A short conclusion ends the article. The supplemental appendix contains technical details and additional simulation results.

2 Inference With Many Regressors

2.1 Framework

Consider the linear model(1) $$y_{i,n}={\mathit{x}}_{i,n}^{\prime }\mathit{\beta }+{\mathit{w}}_{i,n}^{\prime }{\mathit{\gamma }}_n+u_{i,n},\hspace{1em}i=1,\dots ,n,$$(1) where $y_{i,n}$ is a scalar outcome, ${\mathit{x}}_{i,n}$ is a vector of regressors of fixed dimension r, ${\mathit{w}}_{i,n}$ is a vector of covariates whose dimension, q_n , may grow with n, and $u_{i,n}$ is an unobserved error term. Our aim is to perform asymptotically valid inference on $\mathit{\beta }$ that is robust to heteroscedasticity, when ${\mathit{\gamma }}_n$ is high-dimensional, in the sense that q_n is not a vanishing fraction of the sample size. In such a case, the nuisance parameter ${\mathit{\gamma }}_n$ is not consistently estimable, in general.

The (ordinary) least-squares estimator of $\mathit{\beta }$ is$${\widehat{\mathit{\beta }}}_n:={\left(\sum _{i=1}^n{\widehat{\mathit{v}}}_{i,n}{\widehat{\mathit{v}}}_{i,n}^{\prime }\right)}^{-1}\left(\sum _{i=1}^n{\widehat{\mathit{v}}}_{i,n}{y}_{i,n}\right),$$where$$\begin{array}{c}{\widehat{\mathit{v}}}_{i,n}:=\sum _{j=1}^n{({\mathit{M}}_n)}_{i,j} {\mathit{x}}_{j,n},\\ {({\mathit{M}}_n)}_{i,j}:=\{i=j\}-{\mathit{w}}_{i,n}^{\prime }{\left(\sum _{k=1}^n{\mathit{w}}_{k,n}{\mathit{w}}_{k,n}^{\prime }\right)}^{-1}{\mathit{w}}_{j,n},\end{array}$$ and $\{·\}$ denotes the indicator function. We will provide an inference result based on the limit distribution of ${\widehat{\mathit{\beta }}}_n$. We begin by stating a set of high-level conditions that guarantee this distribution to be Gaussian that cover the case where $q_n/n\nrightarrow 0$ as $n\to \infty$. Our point of departure for this is Cattaneo, Jansson, and Newey (2018, Theorem 1) and we closely follow their notation.

Let ${\mathcal{X}}_n:=({\mathit{x}}_{1,n},\dots ,{\mathit{x}}_{n,n})$ and let ${\mathcal{W}}_n$ denote a collection of random variables such that $\mathbb{E}({\mathit{w}}_{i,n}|{\mathcal{W}}_n)={\mathit{w}}_{i,n}$. We introduce$${\epsilon }_{i,n}:=u_{i,n}-e_{i,n},\hspace{1em}{\mathit{V}}_{i,n}:={\mathit{x}}_{i,n}-\mathbb{E}({\mathit{x}}_{i,n}|{\mathcal{W}}_n),$$where $e_{i,n}:=\mathbb{E}(u_{i,n}|{\mathcal{X}}_n,{\mathcal{W}}_n)$, to state our first assumption. We use $|·|$ to denote the cardinality of a set.

Assumption 1

(Sampling). The errors ${\epsilon }_{i,n}$ are uncorrelated across i conditional on ${\mathcal{X}}_n$ and ${\mathcal{W}}_n$, and the collections $\{{\epsilon }_{i,n},{\mathit{V}}_{i,n}:i\in N_g\}$ are independent across g conditional on ${\mathcal{W}}_n$, where $\{N_1,\dots ,N_{G_n}\}$ represents a partition of $\{1,\dots ,n\}$ into G_n sets such that ${\max }_g\left|N_g\right|=O(1)$.

This assumption covers standard randomly sampled data but also repeated-measurement data (such as short panel data) where strata are independent but dependence between observations within the strata is allowed, for example.

The second assumption contains regularity conditions. We let$${\sigma }_{i,n}^2:=\mathbb{E}({\epsilon }_{i,n}^2|{\mathcal{X}}_n,{\mathcal{W}}_n),\hspace{1em}{\tilde{\mathit{V}}}_{i,n}:=\sum _{j=1}^n{({\mathit{M}}_n)}_{i,j} {\mathit{V}}_{j,n},$$denote the Euclidean and Frobenius norms by $\left|\right|·\left|\right|$, and write ${\lambda }_{\min }(·)$ for the minimum eigenvalue of its argument.

Assumption 2

(Design). With probability approaching one ${\sum }_{i=1}^n{\mathit{w}}_{i,n}{\mathit{w}}_{i,n}^{\prime }$ has full rank,$$\begin{array}{c}\underset{i}{\max }\left(\mathbb{E}({\epsilon }_{i,n}^4|{\mathcal{X}}_n,{\mathcal{W}}_n)+\frac{1}{{\sigma }_{i,n}^2}+\mathbb{E}(\left|\right|{\mathit{V}}_{i,n}|{|}^4 |{\mathcal{W}}_n)\right)\\ \hspace{1em}+\frac{1}{{\lambda }_{\min }\left(\frac{\sum _{i=1}^n\mathbb{E}({\tilde{\mathit{V}}}_{i,n}{\tilde{\mathit{V}}}_{i,n}^{\prime }|{\mathcal{W}}_n)}{n}\right)}=O_p(1),\end{array}$$and $\lim \underset{n}{\sup }{q}_n/n<1$.

The rank condition on the design matrix ${\sum }_{i=1}^n{\mathit{w}}_{i,n}{\mathit{w}}_{i,n}^{\prime }$ is standard. Furthermore, given that the slope coefficients on ${\mathit{w}}_{i,n}$ are not of direct interest to us, dropping any covariates that are (perfectly) collinear is not an issue. The second condition contains conventional moment conditions. The third condition, finally, allows for q_n to grow at the same rate as the sample size.

Our setting covers situations where the regression in (1) is a linear-in-parameters mean-square approximation to the conditional expectation ${\mu }_{i,n}:=\mathbb{E}(y_{i,n}|{\mathcal{X}}_n,{\mathcal{W}}_n)$, in the sense that we allow that $e_{i,n}\ne 0$. The third assumption contains conditions on how fast such an approximation should improve. They are expressed in terms of the two constants$${\varrho }_n:=\frac{\sum _{i=1}^n\mathbb{E}(e_{i,n}^2)}{n},\hspace{1em}{\rho }_n:=\frac{\sum _{i=1}^n\mathbb{E}(\mathbb{E}{(e_{i,n}|{\mathcal{W}}_n)}^2)}{n}.$$

The assumption also contains a similar restriction on how well ${\mathit{V}}_{i,n}$ can be approximated by$${\mathit{v}}_{i,n}:={\mathit{x}}_{i,n}-\left(\sum _{j=1}^n\mathbb{E}({\mathit{x}}_{j,n}{\mathit{w}}_{j,n}^{\prime })\right){\left(\sum _{j=1}^n\mathbb{E}({\mathit{w}}_{j,n}{\mathit{w}}_{j,n}^{\prime })\right)}^{-1}{\mathit{w}}_{i,n},$$the deviation of ${\mathit{x}}_{i,n}$ from its population linear projection. This is expressed using the constant$${\chi }_n:=\frac{\sum _{i=1}^n\mathbb{E}(|\left|{\mathit{Q}}_{i,n}\right|{|}^2)}{n},$$ where ${\mathit{Q}}_{i,n}:=\mathbb{E}({\mathit{v}}_{i,n}|{\mathcal{W}}_n).$

Assumption 3

(Approximations). ${\chi }_n=O(1),{\varrho }_n+n({\varrho }_n-{\rho }_n)+n{\chi }_n{\varrho }_n=o(1)$, and ${\max }_i \left|\right|{\widehat{\mathit{v}}}_{i,n}\left|\right|/\sqrt{n}=o_p(1)$.

The last part of this assumption is a high-level negligibility condition on the residuals from the auxiliary regression of ${\mathit{x}}_{i,n}$ on ${\mathit{w}}_{i,n}$. Given that$$\begin{array}{c}\underset{i}{\max }\left({\widehat{\mathit{v}}}_{i,n}^{\prime }{\left(\sum _{j=1}^n{\widehat{\mathit{v}}}_{j,n}{\widehat{\mathit{v}}}_{j,n}^{\prime }\right)}^{-1} {\widehat{\mathit{v}}}_{i,n}\right)\\ \hspace{1em}\le \underset{i}{\max }\frac{\left|\right|{\widehat{\mathit{v}}}_{i,n}|{|}^2}{n} {\Vert {\left(\frac{\sum _{j=1}^n{\widehat{\mathit{v}}}_{j,n}{\widehat{\mathit{v}}}_{j,n}^{\prime }}{n}\right)}^{-1}\Vert }^2\end{array}$$and should vanish in large samples for ${\widehat{\mathit{\beta }}}_n$ to be asymptotically normal (Huber 1973), this requirement appears close to minimal.

Assumptions 1–3 coincide with Assumptions 1–3 in Cattaneo, Jansson, and Newey (2018). Consequently, by their Theorem 1,(2) $${\mathbf{\Omega }}_n^{-1/2}({\widehat{\mathit{\beta }}}_n-\mathit{\beta })\stackrel{d}{\to }\mathit{N}(0,{\mathit{I}}_r)$$(2) as $n\to \infty$, where$${\mathbf{\Omega }}_n:={\left(\sum _{i=1}^n{\widehat{\mathit{v}}}_{i,n}{\widehat{\mathit{v}}}_{i,n}^{\prime }\right)}^{-1}\left(\sum _{i=1}^n{\widehat{\mathit{v}}}_{i,n}{\widehat{\mathit{v}}}_{i,n}^{\prime } {\sigma }_{i,n}^2\right){\left(\sum _{i=1}^n{\widehat{\mathit{v}}}_{i,n}{\widehat{\mathit{v}}}_{i,n}^{\prime }\right)}^{-1},$$ and ${\mathit{I}}_r$ denotes the r × r identity matrix.

2.2 Variance Estimation

Constructing confidence intervals and test statistics based on (2) requires an estimator of ${\mathbf{\Omega }}_n$, and thus of$${\mathbf{\Sigma }}_n:=\sum _{i=1}^n{\widehat{\mathit{v}}}_{i,n}{\widehat{\mathit{v}}}_{i,n}^{\prime } {\sigma }_{i,n}^2.$$

When $q_n/n\nrightarrow 0$ and the errors are permitted to be heteroscedastic, the construction of a consistent estimator is nontrivial. To appreciate the problem, consider the estimator of Eicker (1963, 1967) and White (1980), which uses$${\widehat{\Sigma }}_n:=\sum _{i=1}^n{\widehat{\mathit{v}}}_{i,n}{\widehat{\mathit{v}}}_{i,n}^{\prime } {\widehat{u}}_{i,n}^2,$$where ${\widehat{u}}_{i,n}:={\sum }_{j=1}^n{({\mathit{M}}_n)}_{i,j} (y_{j,n}-{\mathit{x}}_{j,n}^{\prime }{\widehat{\mathit{\beta }}}_n)$ are the least-squares residuals. This estimator is well known to be (conditionally) biased. The bias arises from the sampling noise in the least-squares estimator and can be severe (Chesher and Jewitt 1987). Unless $q_n/n\to 0$ as $n\to \infty$, some observations will remain influential, in the sense that maximal leverage does not vanish. This causes the bias in ${\widehat{\Sigma }}_n$ to persist in large samples, implying that it is inconsistent.

The alternative to ${\widehat{\Sigma }}_n$ that we consider in this article is$${\acute{\Sigma }}_n:=\sum _{i=1}^n{\widehat{\mathit{v}}}_{i,n}{\widehat{\mathit{v}}}_{i,n}^{\prime } (y_{i,n}{\acute{u}}_{i,n}),\hspace{1em}{\acute{u}}_{i,n}:=\frac{{\widehat{u}}_{i,n}}{{({\mathit{M}}_n)}_{i,i}}.$$

As stated, this estimator is well defined provided that$$\underset{i}{\min }{({\mathit{M}}_n)}_{i,i}>0.$$

Notice that ${({\mathit{M}}_n)}_{i,i}=0$ means that the model reserves a parameter for this observation. This implies that the auxiliary regression of the regressors of interest on the other covariates yields a perfect prediction, in the sense that ${\widehat{\mathit{v}}}_{i,n}=0$. Consequently, such an observation does not carry information on $\mathit{\beta }$ and can be dropped. It does not affect the least-squares estimator ${\widehat{\mathit{\beta }}}_n$ and does not contribute to its covariance matrix ${\mathbf{\Omega }}_n$. This is important as perfect prediction of this form arises frequently in empirical work when many dummy variables are included.

Additional conditions are needed to show that ${\acute{\Sigma }}_n$ is consistent. We let$${\tilde{\mathit{Q}}}_{i,n}:=\sum _{j=1}^n{({\mathit{M}}_n)}_{i,j} {\mathit{Q}}_{i,n}$$and collect one such set of conditions in the following assumption.

Assumption 4

(Variance estimation). $n{\varrho }_n=O(1),\mathrm{Pr}({\min }_i{({\mathit{M}}_n)}_{i,i}>0)\to 1$,$$\frac{1}{\underset{i}{\min }{({\mathit{M}}_n)}_{i,i}}=O_p(1),\hspace{1em}\frac{\sum _{i=1}^n\left|\right|{\tilde{\mathit{Q}}}_{i,n}|{|}^4}{n}=O_p(1),$$and ${\max }_i\left|\right|{\mu }_{i,n}\left|\right|/\sqrt{n}=o_p(1)$.

The first part of Assumption 4 is a small-bias condition; it is relevant only when (1) is misspecified, in the sense that $e_{i,n}\ne 0$. In that case it is a strengthening of Assumption 3 only when ${\chi }_n=o(1)$. The conditions on the diagonal entries of the projection matrix are very weak. Providing primitive conditions for them in great generality appears to be difficult. However, when ${\mathit{w}}_{i,n}$ is multivariate normal they follow under $\lim \underset{n}{\sup }{q}_n/n<1$ as stated in Assumption 2 in the same way as in Cattaneo, Jansson, and Newey (2018). In the one-way panel model they hold automatically while Verdier (2020) gives sufficient conditions for them to be satisfied in the two-way model. To understand why the last part of Assumption 4 is needed, note that the (conditional) variance of $y_{i,n} {\acute{u}}_{i,n}$ depends on ${\mu }_{i,n}^2$. The requirement that ${\max }_i{\mu }_{i,n}^2=o_p(n)$ allows to control the variance of ${\acute{\Sigma }}_n.$ Weak moment requirements typically suffice for this condition to be satisfied. The condition on ${\tilde{\mathit{Q}}}_{i,n}$ is used in concordance with the condition on ${\mu }_{i,n}$. One simple sufficient condition for it is that $n{\chi }_n=O(1)$, but it can also be satisfied when ${\chi }_n=O(1)$. Primitive conditions for Assumption 4 in three special cases are given below.

We can now state our consistency result.

Theorem 1

(Inference). Let Assumptions 1–4 hold. Then$${\mathbf{\Sigma }}_n^{-1}{\acute{\Sigma }}_n\stackrel{p}{\to }{\mathit{I}}_r$$as $n\to \infty$.

Theorem 1, combined with the limit result in (2), implies that$${\acute{\Omega }}_n^{-1/2}({\widehat{\mathit{\beta }}}_n-\mathit{\beta })\stackrel{d}{\to }\mathit{N}(0,{\mathit{I}}_r)$$as $n\to \infty$, where$${\acute{\Omega }}_n:={\left(\sum _{i=1}^n{\widehat{\mathit{v}}}_{i,n}{\widehat{\mathit{v}}}_{i,n}^{\prime }\right)}^{-1}\left(\sum _{i=1}^n{\widehat{\mathit{v}}}_{i,n}{\widehat{\mathit{v}}}_{i,n}^{\prime } (y_{i,n} {\acute{u}}_{i,n})\right){\left(\sum _{i=1}^n{\widehat{\mathit{v}}}_{i,n}{\widehat{\mathit{v}}}_{i,n}^{\prime }\right)}^{-1}.$$

This result permits the construction of test statistics that (in large samples) will have correct size and of confidence regions that will exhibit correct coverage.

3 Connections to the Literature

3.1 HC-Class Estimators

The bias in the Eicker (1963, 1967) and White (1980) estimator has led to a variety of modifications to it being proposed that, following MacKinnon and White (1985), are often referred to as the HC-class of covariance-matrix estimators. These estimators are reviewed in Long and Ervin (2000) and MacKinnon (2012). Unfortunately, as shown by Cattaneo, Jansson, and Newey (2018, Theorem 3), none of these alternatives is consistent, in general, when $q_n/n\nrightarrow 0$. We briefly review their main findings on these estimators here.

The first variance estimator, HC1, differs from the conventional estimator, HC0, in that it performs a degrees-of-freedom correction (Hinkley 1977, eq. (2.11)). This estimator will be consistent in the special case where errors are homoscedastic and the covariate design is balanced, that is, when ${({\mathit{M}}_n)}_{1,1}=\cdots ={({\mathit{M}}_n)}_{n,n}$.

The second variance estimator, HC2, uses ² $$\sum _{i=1}^n{\widehat{\mathit{v}}}_{i,n}{\widehat{\mathit{v}}}_{i,n}^{\prime }({\widehat{u}}_{i,n} {\acute{u}}_{i,n})$$as an estimator of ${\mathbf{\Sigma }}_n$ (Horn, Horn, and Duncan 1975). This estimator will be consistent under homoscedasticity. Because$$y_{i,n}{\acute{u}}_{i,n}={\widehat{u}}_{i,n} {\acute{u}}_{i,n}+{\widehat{y}}_{i,n} {\acute{u}}_{i,n},$$ where ${\widehat{y}}_{i,n}:=y_{i,n}-{\widehat{u}}_{i,n}$ are fitted values, ${\acute{\Sigma }}_n$ can be interpreted as a bias-corrected version of the HC2 estimator.

The third variance estimator, HC3, is constructed with$$\sum _{i=1}^n{\widehat{\mathit{v}}}_{i,n}{\widehat{\mathit{v}}}_{i,n}^{\prime }({\acute{u}}_{i,n} {\acute{u}}_{i,n})$$

(MacKinnon and White 1985). While this estimator is inconsistent, its probability limit exceeds ${\mathbf{\Sigma }}_n$ (in the matrix sense). It follows that (under Assumptions 1–3) test procedures based on HC3 will be asymptotically conservative when $\lim \underset{n}{\sup }q/n\in (0,1)$, both under homoscedasticity and heteroscedasticity.

3.2 Bias-Corrected Estimation

Cattaneo, Jansson, and Newey (2018) also considered the estimator$${\stackrel{`}{\Sigma }}_n:=\sum _{i=1}^n{\widehat{\mathit{v}}}_{i,n}{\widehat{\mathit{v}}}_{i,n}^{\prime } \left(\sum _{j=1}^n(({\mathit{M}}_n*{\mathit{M}}_n{)}^{-1}{)}_{i,j} {\widehat{u}}_{j,n}^2\right),$$where ${\mathit{M}}_n*{\mathit{M}}_n$ denotes the elementwise product of the matrix ${\mathit{M}}_n$. This estimator has its origins in work by Hartley, Rao, and Kiefer (1969) and Rao (1970) and can be motivated through an (asymptotic) bias calculation; see also Bera, Suprayitno, and Premaratne (2002), and Anatolyev (2018) for a refinement under homoscedasticity. For ${\stackrel{`}{\Sigma }}_n$ to be well defined, ${\mathit{M}}_n*{\mathit{M}}_n$ needs to be nonsingular. Necessary and sufficient conditions for this to be the case are stated in Mallela (1972) but these are neither simple nor intuitive (Horn, Horn, and Duncan 1975). As noted by Horn and Horn (1975), a simple sufficient condition is that$$\underset{i}{\min }{({\mathit{M}}_n)}_{i,i}>\frac{1}{2}.$$

Depending on the problem at hand it may also be necessary; an example is the one-way panel model.

Cattaneo, Jansson, and Newey (2018, Theorem 4) show that ${\stackrel{`}{\Sigma }}_n$ is consistent for ${\mathbf{\Sigma }}_n$ if(3) $$\mathrm{Pr}\left(\underset{i}{\min }{({\mathit{M}}_n)}_{i,i}>\frac{1}{2}\right)\to 1,\hspace{1em}\frac{1}{\underset{i}{\min }{({\mathit{M}}_n)}_{i,i}-\frac{1}{2}}=O_p(1),$$(3) are added to Assumptions 1–3. Because ${\sum }_{i=1}^n{({\mathit{M}}_n)}_{i,i}=n-q_n$, ${\min }_i{({\mathit{M}}_n)}_{i,i}\le 1-q_n/n$, and so$$\lim \underset{n}{\sup }{q}_n/n<\frac{1}{2}$$ is required for (3) to be satisfied. This, in turn, is a strengthening of the condition that $\lim \underset{n}{\sup }{q}_n/n<1$ in Assumption 2.

Stock and Watson (2008) proposed a covariance-matrix estimator for linear fixed-effect models that is applicable to short panel data. It is based on an explicit calculation of the probability limit of ${\widehat{\Sigma }}_n-{\mathbf{\Sigma }}_n$, adjusting the inconsistent ${\widehat{\Sigma }}_n$ by subtracting from it a plug-in estimator of this limit quantity. As discussed in Cattaneo, Jansson, and Newey (2018), ${\stackrel{`}{\Sigma }}_n$ can be understood to be a generalization of this approach to the generic setting where $q_n/n\nrightarrow 0$.

Like ${\acute{\Sigma }}_n$, the estimator ${\stackrel{`}{\Sigma }}_n$ need not be positive semidefinite in small samples. Being based on estimators of the individual ${\sigma }_{i,n}^2$ that are linear combinations of ${\widehat{u}}_{1,n}^2,\dots ,{\widehat{u}}_{n,n}^2$ it does, however, retain the invariance to changes in ${\mu }_{i,n}$ that is inherent in the HC-class estimators. This is a desirable feature and explains why a restriction on the magnitude of ${\mu }_{i,n}$ is not needed for this estimator to be consistent. Furthermore, jointly estimating ${\sigma }_{1,n}^2,\dots ,{\sigma }_{n,n}^2$ in this way is attractive from an efficiency point of view, as it is asymptotically equivalent to a minimum-norm unbiased estimator (see Hartley, Rao, and Kiefer 1969; Rao 1970 for details).

3.3 Inference on Variance Components

The variance estimator ${\acute{\Sigma }}_n$ is closely related to the work of Kline, Saggio, and Sølvsten (2019). In our context, their proposal is to estimate each ${\sigma }_{i,n}^2$ by the cross-fit (Newey and Robins 2018) estimator$$y_{i,n} {\check{u}}_{i,n},$$where ${\check{u}}_{i,n}$ is the residual for observation i when the slope coefficients are estimated from the sample from which the ith observation has been omitted. When the regression model is correctly specified—that is, when $e_{i,n}=0$—then $y_{i,n} {\check{u}}_{i,n}$ is (conditionally) unbiased, provided that the leave-own-out regression slopes are well defined. This will be the case if maximal leverage is bounded away from one. ³ In this case, from the Sherman–Morrison formula (as in, e.g., Miller 1974),$$y_{i,n} {\check{u}}_{i,n}=y_{i,n}{\acute{u}}_{i,n}+o_p(1),$$ under our assumptions, connecting the cross-fit estimator of Kline, Saggio, and Sølvsten (2019) to our approach.

Kline, Saggio, and Sølvsten (2019) used this device to construct unbiased estimators of quadratic forms and present several applications. One of these (Kline, Saggio, and Sølvsten 2019, Remark 4 and Lemma 5) is a consistency result for the implied covariance-matrix estimator$$\sum _{i=1}^n{\widehat{\mathit{v}}}_{i,n}{\widehat{\mathit{v}}}_{i,n}^{\prime } (y_{i,n} {\check{u}}_{i,n})$$in designs with fixed regressors where errors are independent and maximal leverage is bounded away from unity. This result is established under the assumption that $e_{i,n}=0$ and that ${\max }_i\left|\right|{\mu }_{i,n}|{|}^2=O(1)$, together with standard regularity conditions as those in Assumption 2. It follows that our variance estimator can be seen as a modification of theirs that is targeted to a setting with many control variables. Theorem 1 may then be understood to be an extension of their Lemma 5 to settings with stochastic regressors and (vanishing) misspecification bias, where the regressors can have unbounded support and observations may depend on a growing number of parameters.

The implied interpretation of ${\acute{\Sigma }}_n$ as an (approximate) cross-fit estimator is also useful to highlight an apparent tension between invariance of a covariance-matrix estimator to changes in ${\mu }_{i,n}$ and the possibility for it to be consistent when we have $\lim \underset{n}{\sup }{q}_n/n>\frac{1}{2}$. A cross-fit estimator of ${\sigma }_{i,n}^2$ that is invariant is necessarily of the form ${\dot{u}}_{i,n}{\ddot{u}}_{i,n}$, where ${\dot{u}}_{i,n}$ and ${\ddot{u}}_{i,n}$ are two least-squares residuals, each one obtained from a (conditionally) independent subsample that excludes the ith observation. Because these auxiliary regressions can be based on at most $\lfloor (n-1)/2\rfloor$ and $\lceil (n-1)/2\rceil$ observations such an approach cannot accommodate situations where $q_n/n>\frac{1}{2}$. The alternative estimators $y_{i,n} {\check{u}}_{i,n}$ and $y_{i,n} {\acute{u}}_{i,n}$ circumvent the need for two independent estimators of .. by using the level of the outcome variable, $y_{i,n}$, as a proxy for $u_{i,n}$. This allows to deal with cases where $\lim \underset{n}{\sup }{q}_n/n>\frac{1}{2}$ but makes the variance estimator sensitive to ${\mu }_{i,n}$.

4 Examples

We next provide more primitive conditions in three special cases that fit out general setup. We focus on sufficient conditions for Assumption 4. Cattaneo, Jansson, and Newey (2018) already gave such conditions for the other assumptions—and, notably, for Assumption 3—to hold.

4.1 Partially Linear Model

Suppose that observations on $(y_i,{\mathit{x}}_i,{\mathit{z}}_i)$ are independent and identically distributed. The partially linear regression model states that$$y_i={\mathit{x}}_i^{\prime }\mathit{\beta }+\phi ({\mathit{z}}_i)+{\epsilon }_i,\hspace{1em}\mathbb{E}({\epsilon }_i|\mathit{x},{\mathit{z}}_i)=0,$$for an unknown function $\phi$. A series approximation of order κ_n of $\phi ({\mathit{z}}_i)$ takes the form ${\mathit{w}}_{i,n}^{\prime }{\mathit{\gamma }}_n$ for ${\mathit{w}}_{i,n}=(p_1({\mathit{z}}_i),\dots ,p_{{\kappa }_n}({\mathit{z}}_i))\prime$ and $\{p_1,\dots ,p_{{\kappa }_n}\}$ a collection of basis functions such as orthogonal polynomials. Our estimator ${\widehat{\mathit{\beta }}}_n$, then, is the least-squares estimator of $\mathit{\beta }$ in$$y_{i,n}={\mathit{x}}_i^{\prime }\mathit{\beta }+{\mathit{w}}_{i,n}^{\prime }{\mathit{\gamma }}_n+u_{i,n},\hspace{1em}{u}_{i,n}={\epsilon }_i+\phi ({\mathit{z}}_i)-{\mathit{w}}_{i,n}^{\prime }{\mathit{\gamma }}_n.$$

Note that $\mathbb{E}(u_{i,n}|{\mathit{x}}_i,{\mathit{z}}_i)\ne 0$, in general. Consistency of ${\widehat{\mathit{\beta }}}_n$ requires that ${\kappa }_n\to \infty$ (and, thus, that $q_n\to \infty$) as $n\to \infty$. Here,$$\begin{array}{c}{\varrho }_n=\underset{\mathit{\gamma }}{\min }\mathbb{E}\left(\left|\right|\phi ({\mathit{z}}_i)-{\mathit{w}}_{i,n}^{\prime }\mathit{\gamma }|{|}^2\right),\\ {\chi }_n=\underset{\mathit{\delta }}{\min }\mathbb{E}\left(\left|\right|\mathbb{E}({\mathit{x}}_{i,n}|{\mathit{w}}_{i,n})-\mathit{\delta }\prime {\mathit{w}}_{i,n}|{|}^2\right).\end{array}$$

In this example, the number of covariates can be large when the dimension of ${\mathit{z}}_i$ is large, so that many terms are included in the approximation even for small κ_n , or when the underlying functions are not (assumed to be) very smooth, so that a large κ_n needs to be used to control bias.

Standard smoothness conditions on $\phi$ imply that $n{\varrho }_n=O(1)$ (Newey 1997), yielding the first condition of Assumption 4. The fourth condition can be validated in the same way, by imposing sufficient smoothness on $\mathbb{E}({\mathit{x}}_{i,n}|{\mathit{w}}_{i,n})$ to get $n{\chi }_n=O(1)$, which implies that ${\sum }_{i=1}^n\left|\right|{\tilde{\mathit{Q}}}_{i,n}|{|}^4=O_p(n)$. This is different from (and stronger than) the first set of primitive conditions discussed by Cattaneo, Jansson, and Newey (2018) to validate Assumption 3, where more smoothness in $\mathbb{E}({\mathit{x}}_{i,n}|{\mathit{w}}_{i,n})$ can be used to compensate for less smoothness in $\phi$. Alternatively, if ${\chi }_n\le \mathbb{E}(||\mathbb{E}({\mathit{x}}_{i,n}|{\mathit{w}}_{i,n})|{|}^2)=O(1)$ and a partitioning estimator (Cattaneo and Farrell 2013) is used to approximate $\phi$ then ${\mathit{M}}_n$ is a band matrix. This is also sufficient to reach the desired result. Moreover, in this case, the first and fourth condition of Assumption 4 are implied by the rate requirements in Assumption 3, and by the second set of restrictions for this example given in Cattaneo, Jansson, and Newey (2018). Next, simple sufficient conditions for the requirement that ${\max }_i\left|\right|{\mu }_{i,n}\left|\right|/\sqrt{n}=o_p(1)$ are the moment conditions $\mathbb{E}(|\left|{\mathit{x}}_i\right|{|}^{2+\theta })=O(1)$ and $\mathbb{E}(||\phi ({\mathit{z}}_i)|{|}^{2+\theta })=O(1)$ for some $\theta >0$. While these two conditions are not imposed in Cattaneo, Jansson, and Newey (2018), they would not appear to be overly strong.

4.2 One-Way Model for Panel Data

For double-indexed data $(y_{(g,m)},{\mathit{x}}_{(g,m)})$, the fixed-effect model is$$y_{(g,m)}={\mathit{x}}_{(g,m)}^{\prime } \mathit{\beta }+{\alpha }_g+{\epsilon }_{(g,m)},\hspace{1em}g=1,\dots ,G_n,m=1\dots ,M,$$where α_g is a group-specific intercept and we assume that $\mathbb{E}({\epsilon }_{(g,m)}|{\mathit{x}}_{(g,1)},\dots ,{\mathit{x}}_{(g,M)})=0$. The regressors ${\mathit{x}}_{(g,m)}$ are assumed independent between groups but may be dependent within each group. The errors ${\epsilon }_{(g,m)}$ are independent between groups and (conditionally) uncorrelated within groups. The usual asymptotic embedding here has $G_n\to \infty$ with $M=O(1)$. The number of fixed effects grows at the same rate as the sample size; we have $n=G_n\times M$ and q_n = G_n so that $q_n/n=1/M$, which does not vanish. These conditions fit Assumption 1.

The fixed-effect estimator is the least-squares estimator of $y_{(g,m)}$ on ${\mathit{x}}_{(g,m)}$ and G_n dummy variables that indicate group membership. The estimated coefficients on these dummies are computed from M observations and are not consistent under our asymptotic approximation. In this example ${\mathit{M}}_n={\mathit{I}}_{G_n}\otimes {\mathit{T}}_M$, where ${({\mathit{T}}_M)}_{m,m\prime }:=\{m=m\prime \}-M^{-1}$ is the M × M matrix that transforms observations into deviations from their within-group mean. Consequently,$${\acute{\Sigma }}_n=\frac{M}{M-1}\sum _{g=1}^{G_n}\sum _{m=1}^M{\tilde{\mathit{x}}}_{(g,m)} {\tilde{\mathit{x}}}_{(g,m)}^{\prime }{y}_{(g,m)} ({\tilde{y}}_{(g,m)}-{\tilde{\mathit{x}}}_{(g,m)}^{\prime }{\widehat{\mathit{\beta }}}_n),$$where ${\tilde{y}}_{(g,m)}:={\sum }_{m\prime =1}^M{({\mathit{T}}_M)}_{m,m\prime } y_{(g,m\prime )}$, and ${\tilde{\mathit{x}}}_{(g,m)}$ and ${\tilde{\epsilon }}_{(g,m)}$ are defined in the same way. Further, using that ${\widehat{\mathit{\beta }}}_n=\mathit{\beta }+o_p(1)$,$$\begin{array}{c}{\acute{\Sigma }}_n=\frac{M}{M-1}\sum _{g=1}^{G_n}\sum _{m=1}^M{\tilde{\mathit{x}}}_{(g,m)} {\tilde{\mathit{x}}}_{(g,m)}^{\prime }\\ \hspace{1em}\left({\epsilon }_{(g,m)}{\tilde{\epsilon }}_{(g,m)}+\left({\mathit{x}}_{(g,m)}^{\prime }\mathit{\beta }+{\alpha }_g\right) {\tilde{\epsilon }}_{(g,m)}\right)+o_p(G_n).\end{array}$$

Because $\mathbb{E}({\epsilon }_{(g,m)}{\tilde{\epsilon }}_{(g,m)}|{\mathit{x}}_{(g,1)},\dots ,{\mathit{x}}_{(g,M)})={\sigma }_{(g,m)}^2 (M-1)/M$ the first term constitutes an unbiased estimator of ${\mathbf{\Sigma }}_n$. The second term on the right-hand side is mean zero because the errors are mean-independent of the regressors. Its variance, however, depends on the ${\alpha }_1,\dots ,{\alpha }_{G_n}$.

In a two-wave panel,$${\acute{\Sigma }}_n=\frac{1}{4}\sum _{g=1}^{G_n}\Delta {\mathit{x}}_g \Delta {\mathit{x}}_g^{\prime } (\Delta y_g-\Delta {\mathit{x}}_g^{\prime }{\widehat{\mathit{\beta }}}_n) \Delta y_g,$$where Δ denotes the first-difference operator; so, for example, $\Delta y_g:=y_{(g,2)}-y_{(g,1)}$. In this case, the variance estimator does not depend on the fixed effects. The factor $\frac{1}{4}$ appears because of the de-meaning and is inconsequential. The implied estimator of the covariance matrix ${\mathbf{\Omega }}_n$ is$$\begin{array}{c}{\left(\sum _{g=1}^{G_n}\Delta {\mathit{x}}_g \Delta {\mathit{x}}_g^{\prime }\right)}^{-1}\left(\sum _{g=1}^{G_n}\Delta {\mathit{x}}_g \Delta {\mathit{x}}_g^{\prime } (\Delta y_g-\Delta {\mathit{x}}_g^{\prime }{\widehat{\mathit{\beta }}}_n) \Delta y_g\right)\\ \hspace{1em}\times {\left(\sum _{g=1}^{G_n}\Delta {\mathit{x}}_g \Delta {\mathit{x}}_g^{\prime }\right)}^{-1}.\end{array}$$

The same estimator would be obtained if our procedure would be applied directly to the first-differenced model $\Delta y_g=\Delta {\mathit{x}}_g^{\prime }\mathit{\beta }+\Delta {\epsilon }_g.$ The standard estimator of Eicker (1963, 1967) and White (1980) applied to this model is$$\begin{array}{c}{\left(\sum _{g=1}^{G_n}\Delta {\mathit{x}}_g \Delta {\mathit{x}}_g^{\prime }\right)}^{-1}\left(\sum _{g=1}^{G_n}\Delta {\mathit{x}}_g \Delta {\mathit{x}}_g^{\prime } {(\Delta y_g-\Delta {\mathit{x}}_g^{\prime }{\widehat{\mathit{\beta }}}_n)}^2\right)\\ \hspace{1em}\times {\left(\sum _{g=1}^{G_n}\Delta {\mathit{x}}_g \Delta {\mathit{x}}_g^{\prime }\right)}^{-1},\end{array}$$and is known to be consistent here.

In the one-way panel model our Assumption 4 holds provided that$$\frac{\sum _{g=1}^{G_n}\left|\right|{\alpha }_g|{|}^{2+\theta }}{G_n}=O(1)$$and ${\max }_g{\max }_m\mathbb{E}(|\left|{\mathit{x}}_{(g,m)}\right|{|}^{2+\theta })=O(1)$ for some $\theta >0$. Given that Cattaneo, Jansson, and Newey (2018) imposed that ${\max }_g{\max }_m\mathbb{E}(|\left|{\mathit{x}}_{(g,m)}\right|{|}^4)=O(1)$ to validate Assumption 2, the condition on the fixed effects is the only additional requirement needed for our variance estimator to be consistent in this model.

4.3 Linear Model With Increasing Dimension

Finally, consider the regression model that takes (1) as the data generating process for independent and identically distributed observations $(y_{i,n},{\mathit{x}}_{i,n},{\mathit{w}}_{i,n})$, that is,$$y_{i,n}={\mathit{x}}_{i,n}^{\prime }\mathit{\beta }+{\mathit{w}}_{i,n}^{\prime }{\mathit{\gamma }}_n+u_{i,n},$$with$$n \mathbb{E}(||\mathbb{E}(u_{i,n}|{\mathit{x}}_{i,n},{\mathit{w}}_{i,n})|{|}^2)=n{\varrho }_n=o(1),$$ as in Cattaneo, Jansson, and Newey (2018). The generic nature of this model makes it difficult to specify a single set of simple primitive conditions for our Assumption 4 to hold. First consider the requirement that ${\sum }_{i=1}^n\left|\right|{\tilde{\mathit{Q}}}_{i,n}|{|}^4=O_p(n)$. A sufficient condition here is that$$n {\chi }_n=n \underset{\mathit{\delta }}{\min }\mathbb{E}\left(\left|\right|\mathbb{E}({\mathit{x}}_{i,n}|{\mathit{w}}_{i,n})-\mathit{\delta }\prime {\mathit{w}}_{i,n}|{|}^2\right)=O(1).$$

This will be the case, for example, when ${\mathit{w}}_{i,n}$ is discrete and a saturated regression model is used. It will also hold under smoothness conditions on the function $\mathbb{E}({\mathit{x}}_{i,n}|{\mathit{w}}_{i,n})$ when the ${\mathit{w}}_{i,n}$ are approximating functions, as discussed above. Alternatively, if ${\chi }_n=O(1)$, a sparsity condition on the projection matrix ${\mathit{M}}_n$ can be used. One such condition is ${\max }_i{n}_{i,n}=O_p(1)$, together with ${\sum }_{i=1}^n\mathbb{E}(|\left|{\mathit{Q}}_{i,n}\right|{|}^4)=O(n)$, where we have introduced $n_{i,n}:={\sum }_{j=1}^n\{({\mathit{M}}_n{)}_{i,j}\ne 0\}$. Cattaneo, Jansson, and Newey (2018) showed that the weaker condition ${\max }_i{n}_{i,n}=o_p(n^{1/3})$ can be used to support Assumption 3 under the same moment condition on the ${\mathit{Q}}_{i,n}$. Such a rate would appear difficult to obtain here. Finally, for ${\max }_i\left|\right|{\mu }_{i,n}\left|\right|/\sqrt{n}=o_p(1)$ to hold we again impose that $\mathbb{E}(|\left|{\mathit{x}}_{i,n}\right|{|}^{2+\theta })=O(1)$ for some $\theta >0$. When ${\mathit{w}}_{i,n}$ are approximating functions, a similar moment condition on the function that is being approximated will again suffice. In the case where ${\mathit{w}}_{i,n}$ are just many control variables included in the regression we can again use a sparsity condition. Let ${\kappa }_{i,n}$ denote the number of nuisance parameters on which $y_{i,n}$ depends. Then one alternative sufficient condition is that$$\underset{i}{\max }{\kappa }_{i,n}=O\left(n^{\frac{1}{2}\frac{\theta }{2+\theta }}\right),$$for some $\theta >0$, together with the assumption that the entries of ${\mathit{w}}_{i,n}$ have $2+\theta$ moments. A condition on ${\kappa }_{i,n}$ is different from a condition on $n_{i,n}$, as it only pertains to the regression of $y_{i,n}$ on ${\mathit{x}}_{i,n}$ and ${\mathit{w}}_{i,n}$ and does not restrict the auxiliary regression of ${\mathit{x}}_{i,n}$ on ${\mathit{w}}_{i,n}$. When the covariates are normally distributed or have bounded support, for example, we can allow for ${\max }_i{\kappa }_{i,n}=O\left(\sqrt{n}\right)$.

5 Numerical Illustrations

5.1 Simulations

We present numerical results for a setup taken from Cattaneo, Jansson, and Newey (2018). Data are generated as$$y_{i,n}=x_{i,n}\beta +{\mathit{w}}_{i,n}^{\prime }{\mathit{\gamma }}_n+{\epsilon }_{i,n},$$where $x_{i,n}\sim \text{iid} \mathit{N}(0,1),{\mathit{w}}_{i,n}$ contains a constant term and a collection of $q_n-1$ zero/one dummy variables, and ${\epsilon }_{i,n}\sim \text{iid} \mathit{N}(0,1)$. The dummy variables are drawn independently with success probability π and ${\mathit{\gamma }}_n=0$. The sample size was fixed to n = 700 throughout and we considered $q_n\in \{1,71,141,211,281,351,421,491,561,631\}$. All statistics reported below were computed over 10,000 Monte Carlo replications and all the variables were redrawn in each replication.

We consider three designs that vary in β and π. Design A is the design of Cattaneo, Jansson, and Newey (2018). It has β = 1 and $\pi =0.02$. ⁴ Each dummy variable takes on the value one for about 14 out of 700 observations, on average. Design B is a more sparse design where β = 1 is maintained and π is reduced to 0.01, leading to each dummy switching on for only 7 observations in each replication, on average. Design C, in turn, sets β = 2 and maintains $\pi =0.02$. This implies that the conditional variance of $y_{i,n} {\acute{u}}_{i,n}$ increases by a factor of four.

The results for the three design variations are presented in Tables 1–3. Empirical size of the two-sided t-test of the null that β = 1 and the average width of the corresponding confidence interval are given for all variance estimators discussed; HC0, its modifications HC1, HC2, and HC3, the bias-corrected estimator of Cattaneo, Jansson, and Newey (2018) (HCK) and the estimator presented here (HCA).

Table 1 Design A (β = 1 and $\pi =0.02$).

qn	1	71	141	211	281	351	421	491	561	631
$q_n/n$	0.0014	0.1014	0.2014	0.3014	0.4014	0.5014	0.6014	0.7014	0.8014	0.9014
Rejection frequency of 5%-level test
HC0	0.0502	0.0671	0.0771	0.1050	0.1287	0.1605	0.2054	0.2746	0.3691	0.5309
HC1	0.0495	0.0518	0.0485	0.0534	0.0531	0.0459	0.0426	0.0476	0.0438	0.0446
HC2	0.0498	0.0529	0.0517	0.0567	0.0574	0.0505	0.0495	0.0556	0.0530	0.0589
HC3	0.0495	0.0416	0.0280	0.0241	0.0131	0.0058	0.0022	0.0005	0.0001	0.0000
HCK	0.0498	0.0530	0.0520	0.0588	0.0593	0.0563	0.0598	0.0770	0.0976	0.2688
HCA	0.0519	0.0535	0.0529	0.0570	0.0566	0.0524	0.0521	0.0605	0.0567	0.0674
Average width of 95% confidence interval
HC0	0.1477	0.1484	0.1490	0.1494	0.1502	0.1510	0.1517	0.1525	0.1533	0.1541
HC1	0.1479	0.1567	0.1669	0.1789	0.1944	0.2141	0.2408	0.2797	0.3452	0.4944
HC2	0.1478	0.1560	0.1654	0.1766	0.1909	0.2093	0.2340	0.2698	0.3300	0.4646
HC3	0.1479	0.1645	0.1851	0.2113	0.2468	0.2964	0.3707	0.4937	0.7405	1.4803
HCK	0.1478	0.1559	0.1651	0.1758	0.1894	0.2065	0.2292	0.2606	0.3064	0.3482
HCA	0.1477	0.1557	0.1652	0.1763	0.1907	0.2087	0.2334	0.2688	0.3284	0.4611

Table 2 Design B (β = 1 and $\pi =0.01$).

qn	1	71	141	211	281	351	421	491	561	631
$q_n/n$	0.0014	0.1014	0.2014	0.3014	0.4014	0.5014	0.6014	0.7014	0.8014	0.9014
Rejection frequency of 5%-level test
HC0	0.0479	0.0606	0.0767	0.0963	0.1137	0.1457	0.1907	0.2450	0.3353	0.4639
HC1	0.0474	0.0462	0.0479	0.0471	0.0410	0.0432	0.0390	0.0340	0.0291	0.0256
HC2	0.0477	0.0477	0.0526	0.0523	0.0496	0.0543	0.0553	0.0544	0.0555	0.0584
HC3	0.0474	0.0376	0.0301	0.0227	0.0121	0.0050	0.0026	0.0005	0.0000	0.0000
HCK	0.0477	0.0582	0.0765	0.0963	0.1136	0.1457	0.1907	0.2450	0.3353	0.4639
HCA	0.0492	0.0512	0.0529	0.0540	0.0523	0.0566	0.0571	0.0553	0.0602	0.0705
Average width of 95% confidence interval
HC0	0.1479	0.1492	0.1504	0.1520	0.1537	0.1561	0.1588	0.1621	0.1671	0.1768
HC1	0.1481	0.1575	0.1685	0.1820	0.1989	0.2212	0.2518	0.2971	0.3757	0.5648
HC2	0.1480	0.1560	0.1653	0.1765	0.1905	0.2089	0.2335	0.2691	0.3291	0.4633
HC3	0.1481	0.1645	0.1848	0.2109	0.2459	0.2953	0.3688	0.4908	0.7352	1.4661
HCK	0.1480	0.1501	0.1507	0.1520	0.1537	0.1561	0.1588	0.1621	0.1671	0.1768
HCA	0.1479	0.1557	0.1649	0.1763	0.1903	0.2085	0.2333	0.2686	0.3278	0.4590

Table 3 Design C (β = 2 and $\pi =0.02$).

qn	1	71	141	211	281	351	421	491	561	631
$q_n/n$	0.0014	0.1014	0.2014	0.3014	0.4014	0.5014	0.6014	0.7014	0.8014	0.9014
Rejection frequency of 5%-level test
HC0	0.0505	0.0667	0.0800	0.1052	0.1273	0.1611	0.2045	0.2653	0.3708	0.5273
HC1	0.0501	0.0535	0.0495	0.0508	0.0490	0.0466	0.0464	0.0488	0.0462	0.0453
HC2	0.0504	0.0542	0.0519	0.0536	0.0529	0.0518	0.0521	0.0563	0.0545	0.0583
HC3	0.0501	0.0420	0.0310	0.0209	0.0117	0.0062	0.0019	0.0006	0.0000	0.0000
HCK	0.0504	0.0545	0.0525	0.0557	0.0553	0.0577	0.0631	0.0763	0.1026	0.2635
HCA	0.0555	0.0580	0.0584	0.0612	0.0588	0.0567	0.0624	0.0673	0.0710	0.0935
Average width of 95% confidence interval
HC0	0.1480	0.1483	0.1489	0.1496	0.1502	0.1508	0.1515	0.1523	0.1535	0.1543
HC1	0.1482	0.1566	0.1668	0.1792	0.1944	0.2139	0.2405	0.2794	0.3458	0.4950
HC2	0.1481	0.1559	0.1654	0.1769	0.1909	0.2091	0.2336	0.2696	0.3304	0.4652
HC3	0.1482	0.1644	0.1851	0.2116	0.2468	0.2961	0.3698	0.4933	0.7411	1.4817
HCK	0.1481	0.1558	0.1650	0.1761	0.1894	0.2064	0.2286	0.2601	0.3066	0.3491
HCA	0.1472	0.1552	0.1644	0.1761	0.1900	0.2080	0.2320	0.2673	0.3266	0.4535

As the setup features homoscedastic errors, inference based on HC0 will be liberal for large n when $q_n/n$ is not small (Chesher and Jewitt 1987; Cattaneo, Jansson, and Newey 2018). This is apparent from inspection of the tables. The degrees-of-freedom correction performed by HC1 alleviates most of this concern here. HC2, which is consistent under homoscedasticity, performs quite similarly to HC1. Both corrections do come with (on average) wider confidence intervals. The simulations also illustrate that HC3 yields conservative inference. The rejection frequency (under the null) approaches zero as $q_n/n$ increases. The confidence intervals are also wide. This implies that tests based on HC3 will suffer from low power. The relative inefficiency compared to HC1 and HC2 grows as $q_n/n$ increases.

In Design A, HCK gives close to correct size and confidence intervals of a comparable length as HC1 and HC2 for most of the values of $q_n/n$. As $q_n/n$ increases above $\frac{1}{2}$ this variance estimator does not always exist in each replication. When this happens in our Monte Carlo, HCK defaults to HC0. Consequently, for the larger values of $q_n/n$ in Table 1, the size of the HCK method increases above its nominal size and the average length of the confidence interval equally shrinks relative to HC1 and HC2. In the sparser Design B, the nonexistence of the HCK estimator is more frequent and also arises for smaller values of $q_n/n$. This explains the large overrejection rates observed in Table 2. The results for Design C in Table 3 are similar to those in Table 1 as the performance of HCK is invariant to the scale of the regression slopes.

The simulation results also confirm the theory behind our variance estimator. It yields close to correct inference in all three designs and for all values of $q_n/n$ considered, although some more overrejection is observed for the top-end values of this ratio. The average length of its confidence interval is comparable to those obtained for HC1, HC2, and HCK and is substantially smaller than those for HC3 for large values of $q_n/n$. These results show that this estimator is useful in problems with many covariates when HCK is not available, for example, due to sparseness of the regressor design or because of the presence of high-leverage observations more generally.

The supplemental appendix contains additional simulation results for a partially linear regression model and for a one-way panel model.

5.2 Empirical Example

We next use the different variance estimators available to infer the union membership premium. The data are a balanced panel on 545 working individuals and span 8 years (1980–1987), giving a total of 4360 observations. They are taken from Vella and Verbeek (1998) and are available from the data archive of the Journal of Applied Econometrics (http://qed.econ.queensu.ca/jae/1998-v13.2/vella-verbeek/).

We estimate the union premium as the coefficient from a least-squares regression of log wages on a dummy for union membership, after partialling-out a set of control variables. This set contains average hours worked per week, dummies for marital status and for poor health, a quadratic term in years of experience, and a large collection of dummies, as follows. First, as the data are a panel, both individual fixed effects and year fixed effects are included. Second, as the data contain information on the type of occupation (out of 9 categories) and the industry (out of a total of 12) in which the job is located we also control for these by including sector and industry dummies. Third, we allow for interaction effects between these categorical variables by including occupation-by-year, sector-by-year, and occupation-by-sector dummies as well as occupation-by-sector-by-year dummies. Baseline categories for the year-, sector-, and occupational dummies are chosen and their corresponding dummies are dropped as to avoid a dummy-variable trap. Certain interactions of the occupation and sector dummies never take on the value one and so are equally removed from the analysis. This gives a total of 1086 control variables relative to 4360 observations. Our point estimate of the union premium is 7.43%. This is in line with the literature (see, e.g., Jakubson 1991).

The standard errors on this point estimate, obtained by the various methods discussed, together with the implied 95% confidence intervals for the union premium are collected in Table 4. A total of 330 observations have leverage that exceeds the one-half threshold. 36 of these observations can be perfectly explained by the control variables because of the inclusion of the occupation-by-sector and the occupation-by-sector-by-year dummies and so contain no information on the union premium. The HCK variance estimator could not be computed here. (We note that the invertibility condition fails also after dropping the 36 non-informative observations.) This explains why an entry for this estimator is not available. HC0 yields the smallest standard error on our point estimate (1.72%). HC3 yields the largest (2.32%). The standard error of HCA (1.93%) is roughly in the middle of these two bounds. HC1 and HC2 give similar, slightly larger, standard errors as HCA (1.99% and 1.98%).

Table 4 Inference on the union premium.

	${\widehat{\beta }}_n=0.0743$
	SE	95% CI
HC0	0.0172	[0.0406, 0.1081]
HC1	0.0199	[0.0354, 0.1133]
HC2	0.0198	[0.0356, 0.1131]
HC3	0.0232	[0.0288, 0.1199]
HCA	0.0193	[0.0365, 0.1122]

6 Conclusion

This article has presented a heteroscedasticity-robust covariance-matrix estimator for linear regression models that is consistent under an asymptotic scheme where the number of control variables, q_n , grows at the same rate as the sample size, n. The estimator is similar to the proposal of Kline, Saggio, and Sølvsten (2019) but our consistency result covers more general settings. The estimator complements work by Cattaneo, Jansson, and Newey (2018), who derived inconsistency results for members of the HC-class of variance estimators, proved asymptotic-conservativeness of the HC3 estimator, and presented an alternative variance estimator (based on Hartley, Rao, and Kiefer 1969) that remains consistent when $\lim \underset{n}{\sup }{q}_n/n<\frac{1}{2}$. Under a set of additional high-level conditions, our estimator allows to weaken this restriction to $\lim \underset{n}{\sup }{q}_n/n<1$. Primitive conditions for these where given for partially linear models, fixed-effect panel data models, and generic regression models with increasing dimension. Simulations verify the theoretical properties. An empirical application to estimation of the union premium from panel data was also presented.

The idea underlying our variance estimator can be useful as a device to correct for bias more generally. Kline, Saggio, and Sølvsten (2019) utilized it to bias-correct quadratic forms in fixed-effect estimators. Chernozhukov et al. (2018) and Newey and Robins (2018) utilized related cross-fitting techniques to reduce bias in high-dimensional estimation problems that feature machine-learning estimators. Cattaneo, Jansson, and Ma (2019) characterized the bias in (nonlinear) two-step estimators when the first step features a high-dimensional linear regression.

Acknowledgments

An associate editor and three referees provided very constructive feedback on earlier versions of this article and pointed out a non sequitur in the derivation of a primitive condition. Matias Cattaneo generously shared and discussed his replication material. I am most grateful to each of them for their help.

Supplementary Materials

The supplementary material contains proofs, technical details, and additional simulation results.

Additional information

Funding

Financial support from the European Research Council through grant no. 715787 (MiMo) is gratefully acknowledged.

Notes

1 The leverage of an observation i is defined as the ith diagonal element of the hat matrix, that is, the matrix that transforms the observed outcomes into fitted values. It is bounded between zero and one and measures the influence of the observation on its own fitted value; a larger value reflects a higher influence (see, e.g., Hoaglin and Welsch 1978).

2 We follow Cattaneo, Jansson, and Newey (2018) and construct the HC-class variance estimators using the projection matrix ${\mathit{M}}_n$. The original proposals were made in a context where q_n is treated as fixed and did not differentiate between ${\mathit{x}}_{i,n}$ and ${\mathit{w}}_{i,n}$; they used the annihilator matrix that projects out both sets of variables. The difference is asymptotically negligible under our assumptions.

3 This is a slight strengthening of our requirement in Assumption 4 as we project out only ${\mathit{w}}_{i,n}$ to obtain ${\acute{u}}_{i,n}$ while Kline, Saggio, and Sølvsten (2019) projected out both the regressors ${\mathit{x}}_{i,n}$ and the covariates ${\mathit{w}}_{i,n}$ to obtain ${\check{u}}_{i,n}$.

4 The description of the simulation design in Cattaneo, Jansson, and Newey (2018, p. 1358) contains a typo that would imply that $\pi =0.0062$.