Driscoll-Kraay Covariance Estimation

Why Aggregate the Scores Over Time?

Panel observations can be dependent in two directions: errors may be related across units at a given date, and score contributions may be serially related over time. A covariance estimator that treats every unit-time observation as independent can miss both sources of uncertainty. Driscoll–Kraay estimation handles them by first combining the units’ score contributions within each date, then estimating the long-run covariance of that time series.

The key point is that aggregation happens before the time-series HAC calculation. Each date’s aggregate score includes cross-products between different units, so contemporaneous cross-sectional dependence is not discarded. HAC products across dates then allow the aggregate scores to be serially correlated.

Definition: For the linear panel model $y_{it}=x_{it}'\beta+u_{it}$, let $x_{it}$ be a $k$-vector of regressors and let $\hat u_{it}$ be the fitted residual. Define the time-$t$ aggregate score as $$ \hat g_t=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}x_{it}\hat u_{it}. $$ Its long-run covariance is estimated with a HAC estimator over $t$. The Driscoll–Kraay covariance estimate for the OLS coefficient estimator is $$ \widehat{\operatorname{Var}}(\hat\beta) =\frac{1}{NT}\hat Q^{-1}\hat\Omega\hat Q^{-1}, \qquad \hat Q=\frac{1}{NT}\sum_{t=1}^{T}\sum_{i=1}^{N}x_{it}x_{it}'. $$

The formula is written for the regressors and residuals used in the estimator. For example, if the coefficient estimator is based on a within-transformed regression, its transformed regressors and corresponding residuals are the ones to use in the score. The estimator is a sandwich: $\hat Q^{-1}$ maps score variation into coefficient variation on each side, while $\hat\Omega$ estimates the long-run variation in the aggregate scores.

From the Panel Regression to the HAC Sandwich

Let $Q=E[(NT)^{-1}\sum_{it}x_{it}x_{it}']$, or its appropriate probability limit, and suppose it is nonsingular. The OLS normal equations give the exact decomposition

\sqrt{NT}(\hat\beta-\beta) = \hat Q^{-1}\frac{1}{\sqrt{T}}\sum_{t=1}^{T}g_t, \qquad g_t=\frac{1}{\sqrt{N}}\sum_{i=1}^{N}x_{it}u_{it}.

To see the scaling, $\hat\beta-\beta=\hat Q^{-1}(NT)^{-1}\sum_{it}x_{it}u_{it}$, and $\sum_i x_{it}u_{it}=\sqrt{N}g_t$. Thus multiplying by $\sqrt{NT}$ leaves $T^{-1/2}\sum_t g_t$. If the aggregate score process has long-run covariance $\Omega$, the covariance of its normalized sum approaches $\Omega$. Consequently, the leading covariance of $\hat\beta$ is $(NT)^{-1}Q^{-1}\Omega Q^{-1}$.

Estimate that long-run covariance using the fitted aggregate scores. For lag $\ell$, define

\hat\Gamma_{\ell} =\frac{1}{T}\sum_{t=\ell+1}^{T}\hat g_t\hat g_{t-\ell}'.

A standard Driscoll–Kraay HAC estimate uses Bartlett weights $w_{\ell}=1-\ell/(L+1)$ through truncation lag $L$:

\hat\Omega =\hat\Gamma_0+ \sum_{\ell=1}^{L}w_{\ell} \left(\hat\Gamma_{\ell}+\hat\Gamma_{\ell}'\right).

The two terms at each positive lag make the estimated matrix symmetric. The lag-zero term includes all contemporaneous cross-products among units: expanding $\hat g_t\hat g_t'$ produces products $x_{it}\hat u_{it}\hat u_{jt}x_{jt}'/N$ for every pair $i,j$, not just $i=j$. The lagged terms similarly include cross-products between units at different dates. These features are what permit broad cross-sectional dependence as well as serial dependence in the aggregate scores.

Asymptotic conditions: Driscoll–Kraay inference relies on a growing time dimension and regularity conditions for the aggregate score process. Typical conditions include a nonsingular limiting regressor second-moment matrix, suitable moment bounds, weak enough time dependence for a long-run covariance to exist, and no single unit or date dominating the aggregate. For HAC consistency, the truncation lag grows while remaining small relative to $T$: $L\to\infty$ and $L/T\to0$. Cross-sectional dependence may be extensive, but the resulting aggregate process must still satisfy the time-series regularity conditions.

The construction does not require the errors to have a common variance across units or dates. Nor does it require cross-sectional independence. It does require an appropriate coefficient moment condition, such as $E[x_{it}u_{it}]=0$ for the estimator under consideration; covariance correction does not repair an inconsistent coefficient estimator. The precise assumptions depend on the panel asymptotic regime and estimator.

Worked Calculations

Worked Example: A Scalar HAC Estimate With Serial Dependence

Consider a hypothetical balanced panel with $N=2$, $T=3$, and an intercept-only regression. For each date, both units have outcome $y_{it}=2+e_t/\sqrt{2}$, where $e_1=1$, $e_2=-2$, and $e_3=1$. The fitted intercept is $\hat\beta=2$; the residuals for both units at date $t$ are $e_t/\sqrt{2}$. Their total is zero, as required by the intercept normal equation.

Since $x_{it}=1$, each aggregate score is

\hat g_t=\frac{1}{\sqrt{2}}\sum_{i=1}^{2}\hat u_{it}=e_t, \qquad (\hat g_1,\hat g_2,\hat g_3)=(1,-2,1).

Use $L=1$, so the Bartlett weight is $w_1=1-1/2=1/2$. The lag-zero and lag-one estimates are

\hat\Gamma_0=\frac{1}{3}(1^2+(-2)^2+1^2)=2, \qquad \hat\Gamma_1=\frac{1}{3}\{(-2)(1)+(1)(-2)\}=-\frac{4}{3}.

Therefore,

\hat\Omega=2+\frac{1}{2}\left(-\frac{4}{3}-\frac{4}{3}\right)=\frac{2}{3}.

Here $\hat Q=1$, so the estimated variance of the intercept is $\hat\Omega/(NT)=(2/3)/6=1/9$, and its Driscoll–Kraay standard error is $1/3$, or about $0.3333$ outcome units. The scores sum to zero, $\sum_t\hat g_t=0$, in agreement with the fitted intercept’s normal equation.

Worked Example: Contemporaneous Dependence Across Units

Suppose a hypothetical intercept-only panel has $N=2$ and $T=2$. Both units have outcome $5+u_{it}$, with residuals $(u_{11},u_{21})=(1,1)$ at the first date and $(u_{12},u_{22})=(-1,-1)$ at the second. The fitted intercept is 5, and the residuals sum to zero over the panel.

The aggregate scores are $\hat g_1=\sqrt{2}$ and $\hat g_2=-\sqrt{2}$. Focusing on the lag-zero component (that is, using $L=0$ to isolate contemporaneous variation),

\hat\Gamma_0=\frac{1}{2}\{(\sqrt{2})^2+(-\sqrt{2})^2\}=2.

The cross-products between the two units contribute to this value: at either date, $(u_{1t}+u_{2t})^2/2=2$. If one instead retained only the unitwise squared residuals, the corresponding diagonal contribution would be $\{(1+1)/2+ (1+1)/2\}/2=1$. Thus ignoring the positive cross-unit products would halve this lag-zero estimate. With $\hat Q=1$, the lag-zero Driscoll–Kraay variance estimate is $2/(NT)=2/4=0.5$, in squared outcome units. This calculation isolates the contemporaneous component; it does not claim that a zero-lag estimate is adequate when serial dependence matters.

Worked Example: A Two-Coefficient Sandwich

Suppose a hypothetical fitted panel regression has $N=4$, $T=3$, and two regressors. Its normalized design matrix is $\hat Q=I_2$, and its aggregate fitted scores are $\hat g_1=(1,0)'$, $\hat g_2=(-2,1)'$, and $\hat g_3=(1,-1)'$. The scores sum to the zero vector, consistent with the regression normal equations. Use $L=0$ here to illustrate the matrix calculation, and suppose the coefficient estimate is $\hat\beta=(1.2,-0.4)'$.

The contemporaneous score covariance estimate is the average of the three outer products:

\hat\Omega =\frac{1}{3} \left[ \begin{pmatrix}1&0\\0&0\end{pmatrix} +\begin{pmatrix}4&-2\\-2&1\end{pmatrix} +\begin{pmatrix}1&-1\\-1&1\end{pmatrix} \right] = \begin{pmatrix}2&-1\\-1&2/3\end{pmatrix}.

Because $\hat Q=I_2$ and $NT=12$, the covariance estimate is

\widehat{\operatorname{Var}}(\hat\beta) =\frac{1}{12}\hat\Omega = \begin{pmatrix} 1/6&-1/12\\ -1/12&1/18 \end{pmatrix}.

The standard error for the second coefficient is $\sqrt{1/18}\approx0.2357$ coefficient units. Its asymptotic normal statistic for $H_0:\beta_2=0$ is $-0.4/0.2357\approx-1.697$, giving a two-sided normal-approximation $p$-value of about $0.0898$. At the 5% level, this calculation does not reject the null. The zero-lag truncation is used only to show the matrix sandwich; when serial dependence is present, lagged score products are needed as well.

What the Estimator Does—and Does Not—Protect Against

Driscoll–Kraay estimation is useful when dependence across units makes unitwise or cross-sectionally independent score calculations unreliable, while the time dimension is sufficiently informative to estimate aggregate-score persistence. It allows heteroskedasticity because the score products are not built from a common variance assumption. It allows broad cross-sectional dependence because each aggregate score is formed by summing over units before the outer products are calculated. And it accommodates serial correlation through the weighted lagged products.

These protections concern the covariance estimate, not the specification or identification of the regression. They do not make an invalid exogeneity assumption valid, remove omitted-variable bias, or guarantee accurate inference with a very short time series. The estimator also does not accommodate arbitrary time dependence without limit: the aggregate score process must have a stable enough long-run covariance for time-series HAC estimation to work.

Common Mistake

A common mistake is to compute separate unit-level variance contributions and then add them, reasoning that a large number of units makes cross-unit dependence negligible. That calculation omits the $i\ne j$ cross-products in $\hat g_t\hat g_t'$. If units share shocks, those terms can be substantial, as the second example shows. The correction is to aggregate the score contributions across units at each date first, and then apply the HAC calculation to that aggregate sequence.

Key takeaway: Driscoll–Kraay estimation is a time-series HAC estimator applied to cross-sectionally aggregated panel scores. Its sandwich scaling is $1/(NT)$; its aggregate outer products retain cross-unit dependence, and its lagged products account for serial dependence in the aggregate scores.

Check Your Understanding

Try these questions before reviewing the derivation.

Why are the score contributions aggregated across units before applying the HAC calculation?
Write the lag-$\ell$ aggregate-score covariance estimate and explain why both $\hat\Gamma_\ell$ and $\hat\Gamma_\ell'$ appear in the long-run covariance estimate.
In an intercept-only fitted panel regression, what normal-equation restriction applies to the sum of the aggregate residual scores over time?
What dependence is represented by the cross-products between different units inside a date’s aggregate score outer product?
State the basic bandwidth growth conditions used for HAC consistency, and explain why the covariance correction does not by itself solve endogeneity.