Consistency of Residual-Based Covariance Estimators

Why Residuals Can Replace Unobserved Errors

OLS covariance formulas involve the regression errors, but those errors are unobserved. In practice, covariance estimates therefore use residuals. The key asymptotic question is whether fitting the regression changes the sample covariance calculation by an amount that vanishes as the sample grows. With random covariates, the answer depends on moment conditions: the covariates, errors, and fitted residuals all vary across observations.

We study the fixed-dimensional regression model with independent and identically distributed observations $(y_i,x_i,u_i)$, where $x_i$ is a $k$-vector and

y_i=x_i'\beta+u_i,\qquad i=1,\ldots,n.

Here $\beta$ is a fixed $k$-vector, $E[x_i u_i]=0$, and $Q=E[x_i x_i']$ is finite and positive definite. Write $Q_n=n^{-1}\sum_{i=1}^n x_i x_i'$, let $\hat\beta$ be the OLS estimator when $Q_n$ is invertible, and define $\hat u_i=y_i-x_i'\hat\beta$. The results below concern fixed $k$, not a number of regressors that grows with $n$.

Definition: The heteroskedasticity-robust covariance estimator considered here is the HC0 sandwich estimator $$ \widehat V_{\mathrm{HC0}}=(X'X)^{-1}\left(\sum_{i=1}^n x_i x_i'\hat u_i^2\right)(X'X)^{-1}. $$ Its middle factor is called the estimated “meat.” The corresponding population matrix is $\Omega=E[x_i x_i'u_i^2]$. The estimator of the asymptotic covariance of $\sqrt n(\hat\beta-\beta)$ is $Q_n^{-1}M_nQ_n^{-1}$, where $M_n=n^{-1}\sum_{i=1}^n x_i x_i'\hat u_i^2$; thus $\widehat V_{\mathrm{HC0}}=n^{-1}Q_n^{-1}M_nQ_n^{-1}$.

The distinction between $\widehat V_{\mathrm{HC0}}$ and the covariance estimate for $\sqrt n(\hat\beta-\beta)$ is a factor of $n$. Keeping this scaling explicit helps prevent a common error when comparing a sample covariance estimator with its population limit.

Consistency of the Residual-Based Robust Estimator

A useful sufficient condition for estimating $\Omega$ is $E[\|x_i\|^4]<\infty$ and $E[u_i^2\|x_i\|^2]<\infty$. These assumptions make the relevant products integrable. In particular, they do not require $E[u_i^2]<\infty$, a point that matters when covariates can be small relative to errors.

Theorem: Suppose the observations are iid, $Q=E[x_i x_i']$ is finite and positive definite, $E[x_i u_i]=0$, $E[\|x_i\|^4]<\infty$, and $E[u_i^2\|x_i\|^2]<\infty$. Then OLS is consistent, $M_n\xrightarrow{p}\Omega$, and $$ n\widehat V_{\mathrm{HC0}}\xrightarrow{p}Q^{-1}\Omega Q^{-1}. $$

Proof. The moment conditions imply $E[\|x_i u_i\|]<\infty$, because

E[\|x_i u_i\|] =E[\|x_i\||u_i|] \leq \{E[u_i^2\|x_i\|^2]\}^{1/2}<\infty.

The weak law of large numbers gives $Q_n\xrightarrow{p}Q$ and $n^{-1}\sum_i x_i u_i\xrightarrow{p}0$. Since $Q$ is positive definite, $Q_n$ is invertible with probability approaching one, and $Q_n^{-1}\xrightarrow{p}Q^{-1}$. The OLS identity

\hat\beta-\beta=Q_n^{-1}\left(\frac{1}{n}\sum_{i=1}^n x_i u_i\right)

therefore gives $\hat\beta\xrightarrow{p}\beta$. Put $\delta_n=\hat\beta-\beta$. Then $\hat u_i=u_i-x_i'\delta_n$, so

M_n-\frac{1}{n}\sum_{i=1}^n x_i x_i'u_i^2 =\frac{1}{n}\sum_{i=1}^n x_i x_i'\{-2u_i x_i'\delta_n+(x_i'\delta_n)^2\}.

Using the operator norm and $|x_i'\delta_n|\leq\|x_i\|\|\delta_n\|$, the norm of this difference is at most

2\|\delta_n\|\left(\frac{1}{n}\sum_{i=1}^n\|x_i\|^3|u_i|\right) +\|\delta_n\|^2\left(\frac{1}{n}\sum_{i=1}^n\|x_i\|^4\right).

Both empirical averages converge in probability to finite expectations. For the first one, finiteness follows directly from the stated assumptions:

E[\|x_i\|^3|u_i|] \leq \{E[u_i^2\|x_i\|^2]\}^{1/2}\{E[\|x_i\|^4]\}^{1/2}<\infty.

Because $\delta_n\xrightarrow{p}0$, the displayed norm bound converges in probability to zero. Also, the weak law applies entry by entry to $x_i x_i'u_i^2$, since $E[\|x_i x_i'u_i^2\|]=E[\|x_i\|^2u_i^2]<\infty$. Hence $M_n\xrightarrow{p}\Omega$. Finally, continuity of matrix multiplication and inversion gives $Q_n^{-1}M_nQ_n^{-1}\xrightarrow{p}Q^{-1}\Omega Q^{-1}$. This proves the result. $\square$

The proof isolates the important replacement step: the residual error is controlled by $\|\hat\beta-\beta\|$ multiplied by empirical averages that remain bounded in probability. The moment assumptions are sufficient for those averages to be well behaved; OLS consistency alone would not establish this.

Worked Example: A Heteroskedastic Random-Design Sandwich

Let $x_i=(1,z_i)'$, where $P(z_i=0)=P(z_i=1)=1/2$. Suppose $u_i$, conditional on $z_i$, has mean zero and variance $1+z_i$, and observations are iid. All the moment conditions of the theorem hold because $x_i$ and the conditional error variances are bounded.

First calculate the population Gram matrix:

Q=E[x_i x_i'] =\frac12\begin{pmatrix}1&0\\0&0\end{pmatrix} +\frac12\begin{pmatrix}1&1\\1&1\end{pmatrix} =\begin{pmatrix}1&1/2\\1/2&1/2\end{pmatrix}.

Its determinant is $1/2-1/4=1/4$, so it is positive definite, and

Q^{-1}=\begin{pmatrix}2&-2\\-2&4\end{pmatrix}.

For the meat, the conditional variance is 1 when $z_i=0$ and 2 when $z_i=1$. Thus

\Omega=E[x_i x_i'u_i^2] =\frac12\begin{pmatrix}1&0\\0&0\end{pmatrix} +\frac12(2)\begin{pmatrix}1&1\\1&1\end{pmatrix} =\begin{pmatrix}3/2&1\\1&1\end{pmatrix}.

Multiplying gives $Q^{-1}\Omega=\begin{pmatrix}1&0\\1&2\end{pmatrix}$, and therefore

Q^{-1}\Omega Q^{-1} =\begin{pmatrix}1&0\\1&2\end{pmatrix} \begin{pmatrix}2&-2\\-2&4\end{pmatrix} =\begin{pmatrix}2&-2\\-2&6\end{pmatrix}.

The theorem says that $n\widehat V_{\mathrm{HC0}}$ converges in probability to this matrix. The second diagonal entry, 6, is the limiting variance for the slope component of $\sqrt n(\hat\beta-\beta)$; it is not the limiting variance of the unscaled slope estimator, whose variance is of order $1/n$.

Residual Variance Under Homoskedasticity

If the errors have a common conditional variance, the familiar homoskedastic covariance estimator uses the residual sum of squares in place of the robust meat. This requires a finite error variance itself, a stronger condition than the weighted second moment used above.

Theorem: Suppose the observations are iid, $Q=E[x_i x_i']$ is finite and positive definite, $E[x_i u_i]=0$, $E[\|x_i\|^2]<\infty$, $E[u_i^2]<\infty$, and $E[u_i^2\mid x_i]=\sigma^2<\infty$ almost surely. Then $$ \hat\sigma^2=\frac{1}{n}\sum_{i=1}^n\hat u_i^2\xrightarrow{p}\sigma^2. $$ Consequently, for fixed $k$, replacing $n$ by $n-k$ in the denominator has the same probability limit.

Proof. OLS consistency follows from the same Gram-matrix identity as above: $Q_n\xrightarrow{p}Q$, $n^{-1}\sum_i x_i u_i\xrightarrow{p}0$, and hence $\delta_n=\hat\beta-\beta\xrightarrow{p}0$. Expand the residual sum of squares:

\frac{1}{n}\sum_{i=1}^n\hat u_i^2 =\frac{1}{n}\sum_{i=1}^n u_i^2 -2\delta_n'\left(\frac{1}{n}\sum_{i=1}^n x_i u_i\right) +\delta_n'\left(\frac{1}{n}\sum_{i=1}^n x_i x_i'\right)\delta_n.

The first term converges in probability to $E[u_i^2]=\sigma^2$, by the weak law and the conditional-variance assumption. The middle term converges to zero: its score average converges to zero, and $\delta_n\xrightarrow{p}0$. The final term also converges to zero, because $Q_n\xrightarrow{p}Q$ and $\delta_n\xrightarrow{p}0$. This proves the stated limit. Since $n/(n-k)\to1$, the degrees-of-freedom version has the same limit. $\square$

Under these homoskedastic conditions, $\Omega=E[x_i x_i'E[u_i^2\mid x_i]]=\sigma^2Q$. The robust and homoskedastic asymptotic covariance matrices therefore agree:

Q^{-1}\Omega Q^{-1}=\sigma^2Q^{-1}.

The agreement depends on homoskedasticity. It is not a general identity for random-design regression.

Worked Example: Residual Variance in an Intercept-Only Model

Let $x_i=1$ and let $u_i$ be iid with $P(u_i=2)=1/2$ and $P(u_i=-2)=1/2$. The model is $y_i=\beta+u_i$. Here $E[u_i]=0$ and $E[u_i^2]=4$, so the conditional variance is the constant $\sigma^2=4$. OLS gives $\hat\beta=\bar y=\beta+\bar u$, and each residual is $\hat u_i=u_i-\bar u$.

The residual variance estimate satisfies the exact identity

\frac{1}{n}\sum_{i=1}^n\hat u_i^2 =\frac{1}{n}\sum_{i=1}^n(u_i-\bar u)^2 =\frac{1}{n}\sum_{i=1}^n u_i^2-\bar u^2.

Since every $u_i^2=4$, the first term is exactly 4. Also, $E[\bar u]=0$ and $\operatorname{Var}(\bar u)=4/n$, so $\bar u\xrightarrow{p}0$ by Chebyshev's inequality. Thus the estimate equals $4-\bar u^2$ and converges in probability to 4. The homoskedastic covariance estimate for $\hat\beta$ is $\hat\sigma^2/n$, which is consequently asymptotic to $4/n$.

Worked Example: Why Heteroskedasticity Changes the Limit

Consider a regression through the origin, $y_i=x_i\beta+u_i$, with $P(x_i=1)=P(x_i=2)=1/2$. Conditional on $x_i$, let the error have mean zero and variance 1 when $x_i=1$, and variance 9 when $x_i=2$. These bounded covariates and finite conditional variances satisfy the robust-estimator moment conditions.

The scalar Gram matrix is $Q=E[x_i^2]=\tfrac12(1)+\tfrac12(4)=5/2$. The robust meat limit is

\Omega=E[x_i^2u_i^2] =\frac12(1)(1)+\frac12(4)(9)=\frac{37}{2}.

Therefore the robust limit for $n\widehat V_{\mathrm{HC0}}$ is

Q^{-2}\Omega =\frac{37/2}{(5/2)^2} =\frac{74}{25}=2.96.

By contrast, the pooled error variance is $E[u_i^2]=\tfrac12(1)+\tfrac12(9)=5$. Using this pooled variance in a homoskedastic formula would give $E[u_i^2]/Q=5/(5/2)=2$, not 2.96. The discrepancy occurs because the larger error variance is paired with the larger squared covariate, which gives those observations more weight in the OLS score variance.

Moment Conditions and Common Pitfalls

A frequent mistake is to use a moment condition that is stronger than the proof needs—or to claim a stronger condition follows from assumptions that do not imply it. For the robust-meat proof, the cross term is controlled using $E[\|x_i\|^3|u_i|]$. It is finite by applying Cauchy–Schwarz to $|u_i|\|x_i\|$ and $\|x_i\|^2$, as shown in the proof. One must not instead assert that $E[u_i^2]$ is finite merely because $E[u_i^2\|x_i\|^2]$ is finite.

For example, let $z$ take values $2^{-j}$ with probability $2^{-j}$, for $j=1,2,\ldots$, and let $\varepsilon$ be an independent random sign, with probabilities $1/2$ each. Set $u=\varepsilon/z$ and $x=z$. Then $E[xu]=0$, while

E[u^2x^2]=1,\qquad E[x^4]=\sum_{j=1}^\infty 2^{-5j}<\infty,\qquad E[u^2]=\sum_{j=1}^\infty 2^{2j}2^{-j}=\sum_{j=1}^\infty 2^j=\infty.

Thus the robust theorem's weighted moment conditions can hold even when the unweighted error variance is infinite. In contrast, the residual-variance theorem explicitly assumes $E[u_i^2]<\infty$, as it must to use the weak law for $n^{-1}\sum_i u_i^2$.

A second pitfall is to conclude that residual substitution is valid simply because $\hat\beta$ is consistent. The expansion of the residual meat shows what else is needed: the empirical averages multiplying $\|\hat\beta-\beta\|$ and its square must be bounded in probability. Finally, do not describe the homoskedastic covariance estimate as robust to arbitrary conditional variance patterns. The two formulas have the same limit under the stated homoskedastic condition, but the heteroskedastic example shows that they can differ otherwise.

Key takeaway: Residual-based covariance estimation works when OLS coefficient error is small and the relevant covariate–error products have finite moments. For the robust sandwich, the crucial weighted condition is $E[u_i^2\|x_i\|^2]<\infty$; consistency of a pooled residual variance requires the separate condition $E[u_i^2]<\infty$.

Check Your Understanding

Answer each question using the assumptions and calculations in this tutorial.

Which empirical averages control the difference between the residual-based meat and the unobservable error-based meat?
Why does positive definiteness of $Q$, together with $Q_n\xrightarrow{p}Q$, ensure that $Q_n^{-1}\xrightarrow{p}Q^{-1}$ with probability approaching one?
Explain why $E[u_i^2\|x_i\|^2]<\infty$ does not, by itself, imply $E[u_i^2]<\infty$.
Under what conditional variance condition do the homoskedastic and robust asymptotic covariance matrices agree?
In the through-origin example with $x_i\in\{1,2\}$, which feature of the design and error variance produces the difference between 2.96 and 2?

Consistency of Residual-Based Covariance Estimators

What You'll Learn