EconBase
← Back to paper

Variance Estimation with Dependence and Heterogeneous Means

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.

45,952 characters · 6 sections · 16 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Variance Estimation with Dependence and Heterogeneous Means

abstractThis paper considers the problem of estimating the variance of a sum of a triangular array of random vectors with heterogeneous means. When random vectors exhibit two-way cluster dependence or weak dependence, standard variance estimators designed under homogeneous means can underestimate the true variance, which results in subsequent tests being oversized. To restore validity, this paper proposes a simple conservative variance estimator robust to heterogeneous means and shows its asymptotic validity.

\vskip 0.2in

Introduction

Heterogeneous means arise naturally in a variety of statistical settings in which only aggregate moment restrictions are imposed. A prominent example is in design-based settings (e.g., abadie2020sampling, abadie2023should) where researchers condition on heterogeneous unit-specific potential outcomes. In such settings, unit-specific scores need not be mean zero, even though the sum over scores is mean zero. Further, time series analysis is often concerned with nonstationarity that includes heterogeneous means.

Consider a triangular array of random variables $\left\{ Y_{n,i}\right\} , n\geq1$, and a hypothesis test of the mean value i.e., $H_0: \frac{1}{n}\sum_{i=1}^{n}E\left[Y_{n,i}\right]=0$. Many statistical tests rely on consistently estimating the variance, which is usually done by calculating the standard error of the demeaned random variables. For instance, when observations are independent (though not necessarily identically distributed), a variance estimator for the scaled observed mean $\frac{1}{\sqrt{n}}\sum_{i=1}^{n}Y_{n,i}$ is $\frac{1}{n}\sum_{i=1}^{n}\left(Y_{n,i}-\frac{1}{n}\sum_{i=1}^{n}Y_{n,i}\right)^{2}$. Suppose we have heterogeneous means where $E\left[Y_{n,i}\right]$ may differ across observations indexed by $i$, but $\frac{1}{n}\sum_{i=1}^{n}E\left[Y_{n,i}\right]=0$. With the law of large numbers, we still have $\frac{1}{n}\sum_{i=1}^{n}Y_{n,i}\xrightarrow{p}0$, but the plug-in variance estimator targets an estimand of the form:

equation[equation omitted — 228 chars of source]

Independence across observations results in a conservative variance estimator that overestimates the true variance by $\frac{1}{\sqrt{n}}\sum_{i=1}^{n}E\left[Y_{n,i}\right]^{2}$ asymptotically, which still ensures that a standard hypothesis test based on a central limit theorem continues to control size.

However, this paper shows that once there is dependence, it is possible to construct adversarial data generating processes (DGP) such that this conservativeness no longer holds, and hence standard tests can no longer control size. A canonical setting primarily used for exposition in this paper is a panel dataset with time indexed by $t$ and cross sectional clusters indexed by $g$. We allow observations within a cross sectional cluster to be arbitrarily dependent: in a balanced panel where a researcher observes individuals over time, and observations across time are allowed to be arbitrarily correlated for a given individual, each cluster $g$ is defined to be an individual. We further allow serial correlation across cross sectional clusters over time, albeit at a decaying rate. The issue of heterogeneous means has not been explored in such a setting to the best of my knowledge, so this paper first presents results on anticonservativeness of standard plug-in variance estimators, then presents a novel variance estimator that restores validity.

The following (simpler) time series example illustrates this paper's key insight into the problem and a simple solution to restore validity.

exampleSuppose we have an $m$-dependent series, where $m=1$, and $3$ time periods. This example can be extended to the $T\rightarrow\infty$ case for a sequence of heterogeneous means. The variance for the sum of a scalar random variable $y_{t}$ is: \begin{equation} Var\left(\sum_{t=1}^{3}y_{t}\right) =\sum_{t=1}^{3}Var\left(y_{t}\right)+2\sum_{t=1}^{2}Cov\left(y_{t},y_{t+1}\right). \end{equation} When using a standard heteroskedasticity and autocorrelation robust (HAR) variance estimator, we target $\sum_{t=1}^{3}E\left[y_{t}^{2}\right]+2\sum_{t=1}^{2}E\left[y_{t}y_{t+1}\right]$, but \begin{align*} D_{1} & :=\sum_{t=1}^{3}E\left[y_{t}^{2}\right]+2\sum_{t=1}^{2}E\left[y_{t}y_{t+1}\right]-Var\left(\sum_{t=1}^{3}y_{t}\right)\\ & =\sum_{t=1}^{3}E\left[y_{t}\right]^{2}+2\sum_{t=1}^{2}E\left[y_{t}\right]E\left[y_{t+1}\right]. \tag{\stepcounter{equation}\theequation} \end{align*} The expression $D_{1}$ above can be negative while having $\sum_{t=1}^3 E\left[y_{t}\right]=0$ when $E[y_1] = 0.5, E[y_2] = -1, E[y_3]= 0.5$. Then, $\sum_{t=1}^{3}E\left[y_{t}\right]^{2}+2\sum_{t=1}^{2}E\left[y_{t}\right]E\left[y_{t+1}\right]=1.5-2=-0.5$. To restore validity, we may add (a scaled version of) $\sum_{t=1}^{3}y_{t}^{2}$ to the variance estimator, so that we target $2\sum_{t=1}^{3}E\left[y_{t}^{2}\right]+2\sum_{t=1}^{2}E\left[y_{t}y_{t+1}\right]$ instead. Using this expression, \begin{align*} D_{2} & :=2\sum_{t=1}^{3}E\left[y_{t}^{2}\right]+2\sum_{t=1}^{2}E\left[y_{t}y_{t+1}\right]-Var\left(\sum_{t=1}^{3}y_{t}\right)\\ & =2\sum_{t=1}^{3}E\left[y_{t}\right]^{2}+2\sum_{t=1}^{2}E\left[y_{t}\right]E\left[y_{t+1}\right]+\sum_{t=1}^{3}Var\left(y_{t}\right)\\ & \geq\sum_{t=1}^{2}E\left[y_{t}\right]^{2}+2\sum_{t=1}^{2}E\left[y_{t}\right]E\left[y_{t+1}\right]+\sum_{t=1}^{2}E\left[y_{t+1}\right]^{2}+\sum_{t=1}^{3}Var\left(y_{t}\right)\geq 0, \tag{\stepcounter{equation}\theequation} \end{align*} so the variance estimand is conservative and robust to mean heterogeneity. When the means are homogeneous, we would overestimate the true variance by $D_{2}=\sum_{t=1}^{3}Var\left(y_{t}\right)$.

Abstracting from the above example with $T=3$, for an autoregressive process with coefficient $\rho$ and constant means, $Var\left(\sum_{t=1}^{T}y_{t}\right)=\sum_{t=1}^{T}Var\left(y_{t}\right)/\left(1-\rho^{2}\right)$, so the proposed variance estimator would overestimate the true variance by $\left(\sum_{t=1}^{T}Var\left(y_{t}\right)\right)/Var\left(\sum_{t=1}^{T}y_{t}\right)$ $=1-\rho^{2}$. This conservativeness converges to zero for local-to-unity processes where $\rho^{2}\rightarrow1$, but with low serial correlation where $\rho\rightarrow0$, the estimated variance is close to double the true variance. Since the estimator overestimates the true variance by at most a factor of two, the resulting test is still consistent in that, as $T\rightarrow\infty$, the power of the test goes to one. Further, this exposition on conservativeness assumes that there are constant means. Instead, if means are heterogeneous, simulations in this paper show that the proposed variance estimator can numerically achieve rejection rates that are close to the nominal rate.

The general setting in this paper accommodates arbitrary dependence within clusters, and weak dependence across clusters arising from serial correlation or otherwise. I show that the intuition from the above example extends to more general settings, where dependence can result in a variance estimand that underestimates the true variance. By adding an analog of a unit-specific second moment, we can obtain a novel conservative variance estimator that asymptotically guarantees size control.

The paper contributes to three strands of the literature on variance estimation under dependence. First, this paper relates to work on cluster-robust and HAR variance estimation in dependent arrays. xu2024clustering show that, with heterogeneous means and two-way cluster dependence in a design-based setting, the standard two-way cluster robust variance estimator in cameron2011robust (CGM) is anticonservative, so this paper extends the result to cross-cluster weak dependence. Anticonservativeness of a plug-in variance in time series is also new to the best of my knowledge. chiang2024standard (CHS) proposed a variance estimator that accommodates cross-cluster serial dependence in panels, generalizing the two-way cluster variance estimator in CGM, and a standard HAR variance estimator, say from newey1987asimple. However, the CHS variance estimator is not robust to mean heterogeneity. Consequently, this paper exploits the insight from davezies2021empirical that allows CGM to be conservative in two-way clustering, and extends it to the more general setting considered in CHS.

Second, the paper complements recent work that permits distributional heterogeneity in time series. For example, chan2022optimal and casini2023theory allow time-varying means and second moments. chan2022optimal achieves consistency of the long-run variance by differencing, which does not have a natural analog in the present setting with cross-sectional heterogeneity. casini2023theory maintains that the mean function is sufficiently regular to be consistently estimated so that it can be removed prior to applying HAR methods. In contrast, the present paper imposes no regularity structure on the mean sequence. Instead of attempting to consistently estimate the long-run variance under nonstationarity, this paper modifies the variance estimand so as to guarantee asymptotic conservativeness under arbitrary mean heterogeneity.

Third, the paper contributes to the theory of dependent arrays beyond separately exchangeable structures. Much of the two-way clustering literature relies on separate exchangeability through an Aldous–Hoover representation, which restricts admissible forms of heterogeneity. CHS generalized the setting to allow serial correlation, but nonetheless required a representation to exploit projection-based proof techniques. The present paper shows that the CHS representation can be restrictive and rules out reasonable DGP. By leveraging the limit theory developed in kojevnikov2021limit (KMS), the results here do not require the representation used in CHS, thereby permitting more general DGPs. While KMS accommodates triangular arrays indexed by the population size, it maintains homogeneity of means across observations. The present work relaxes this restriction and develops variance theory for arrays with heterogeneous means and weak cross-cluster dependence. Compared with alternative non-exchangeable approaches in two-way clustering such as yap2025asymptotic, the framework here additionally allows serial dependence across clusters.

Setting

Consider a two-way cluster dependent setting, where one dimension exhibits cross-cluster weak dependence. To keep the exposition concrete, suppose we have a panel with both a cross sectional and a time series component. The time series component contains time indices $t\in\left\{ 1,2,\cdots,T\right\} $ for every observation, and observations within the same time period may be arbitrarily correlated. The cross sectional component can be partitioned on some dimension $G$, with cluster indices $g\in\left\{ 1,2,\cdots,G\right\} $ such that every observation within the same $g$ cluster may be arbitrarily correlated, even across time. Let $i\in\left\{ 1,2,\cdots,n\right\} =:N_{n}$ index observations, $g(i)$ denote its cross sectional cluster, and $t(i)$ denote its time period. For example, in a canonical balanced panel with each cross sectional unit observed in every time period, we have $i=\left(g,t\right)$. While the setting in this paper permits unbalanced panels, the balanced panel with $G \asymp T$ is used as a running example.

I consider a triangular array of multivariate random vectors $\left\{ Y_{n,i}\right\} ,n\geq1,Y_{n,i}\in\mathbb{R}^{v}$, where $n$ denotes the number of observations. We define the moment of this random vector as $\left\Vert Y_{n,i}\right\Vert _{p}=\left(E\left[\left\Vert Y_{n,i}\right\Vert ^{p}\right]\right)^{1/p}$. All observations $i$ may be serially correlated. The distance $d_{n}\left(i,j\right)$ between any two observations $i$ and $j$ is 0 if they share a cluster on any dimension, and the difference in their time index otherwise.

Let $N_{t}^{T}$ denote the set of observations in time $t$, $N_{g}^{G}$ denote the set of observations in cluster $g$, and $N_{t,g}^{T\cap G}:=N_{t}^{T}\cap N_{g}^{G}$ the set of observations in both time $t$ and cluster $g$. Define the following objects to measure how concentrated observations are across clusters and time:

align*[align* omitted — 818 chars of source]

These objects can be interpreted in the following manner. $\delta_{n}^{\partial}\left(s;k\right)$ tells us the average number of observations that are $s$ periods away, while $\Delta_{n}\left(s,m;k\right)$ tells us the number of observations between $m$ and $s$ periods away (inclusive), both raised to the power of $k$. Finally, $c_{n}\left(s,m;k\right)$ gives us an upper bound on the convolution of $\delta_{n}^{\partial}\left(s;k\right)$ and $\Delta_{n}\left(s,m;k\right)$ when H\"older's inequality is applied, which is interpretable as a measure of how quickly neighborhoods grow and overlap. In the running example of a balanced panel, $c_n (s,m;k) \asymp (m-s)^k G^k$.

Sufficient conditions for convergences are often stated in terms of $c_{n}\left(s,m;k\right)$. Variance estimators can be written as a sum of products of random vectors, so finding a variance bound on the variance estimators requires bounding a product of four random vectors. With a bound on covariances, it remains to count and bound the number of correlated terms, resulting in a convolution of $\delta_{n}^{\partial}\left(s;k\right)$ and $\Delta_{n}\left(s,m;k\right)$ objects.

The rest of this section carefully characterizes $\psi$-dependence, which is the sense in which observations are weakly dependent across clusters and across time. Let $\mathcal{L}_{v,a}$ denote the collection of bounded Lipschitz real functions on $\mathbb{R}^{v\times a}$:

equation[equation omitted — 151 chars of source]

where $Lip(f)$ is the Lipschitz constant of $f$ i.e., the smallest constant $C$ such that $\left|f(x)-f(y)\right|\leq Cd_{n}(x,y)$ for all $x,y\in\mathbb{R}^{v\times a}$, and $\left\Vert \cdot\right\Vert $ is the sup-norm $\left\Vert f\right\Vert _{\infty}=\sup_{x}\left|f(x)\right|$. For any positive integers $a,b,s$, consider two sets of observations (of size $a$ and $b$) with distance between each other of at least $s$. Let $\mathcal{P}_{n}\left(a,b;s\right)$ denote the collection of all such pairs:

equation[equation omitted — 173 chars of source]

where

equation[equation omitted — 86 chars of source]
definition($\psi$ dependence). A triangular array $\left\{ Y_{n,i}\right\} ,n\geq1,Y_{n,i}\in\mathbb{R}^{v}$ is conditionally $\psi$-dependent, if for each $n\in\mathbb{N}$, there exist a measurable sequence $\theta_{n}=\left\{ \theta_{n,s}\right\}_{s\geq0},\theta_{n,0}=1$, and a collection of nonrandom functions $\left(\psi_{a,b}\right)_{a,b\in\mathbb{N}},\psi_{a,b}:\mathcal{L}_{v,a}\times\mathcal{L}_{v,b}\rightarrow[0,\infty)$, such that for all $\left(A,B\right)\in\mathcal{P}_{n}\left(a,b;s\right)$ with $s>0$ and all $f\in\mathcal{L}_{v,a}$ and $g\in\mathcal{L}_{v,b}$, \begin{equation} \left|Cov\left(f\left(Y_{n,A}\right),g\left(Y_{n,B}\right)\right)\right|\leq\psi_{a,b}\left(f,g\right)\theta_{n,s}, \end{equation} where $Y_{n,A} := (Y_{n,i})_{i \in A}$. We call the sequence $\theta_{n}$ the dependence coefficients of $\left\{ Y_{n,i}\right\} $.

Due to Section 2.4.1 of KMS, this $\psi$ dependence setting generalizes the more familiar strong-mixing processes used in the two-way cluster dependence panel literature with serial correlation. In particular, $\psi$-dependence only requires the covariance decay for Lipschitz functions rather than the entire $\sigma$-field.

The setup here generalizes the CHS setting that requires a representation where $Y_{gt}=f\left(\alpha_{g},\gamma_{t},\varepsilon_{gt}\right)$, for iid $\alpha_{g},\varepsilon_{gt}$ and $\gamma_t$ strictly stationary, which generalizes a standard Aldous-Hoover representation that additionally requires $\gamma_t$ iid. The following example is an example of a DGP that satisfies the conditions of this paper, and has homogeneous means, but does not permit the above representation.

exampleSuppose every firm $g$ is observed for $T$ time periods, each indexed by $t$. Data is generated from the following process. There is an innovation sequence $\left(\nu_{gt}\right)_{t\in\mathbb{Z}}$ iid $N\left(0,1\right)$ over $g,t$. For some $|\rho|\in(0,1)$, $D_{gt}=\rho D_{gt-1}+\nu_{gt},t\in\mathbb{Z}$, with $D_{gt}$ stationary and processes are independent across $g$. Hence, within a firm, $\left\{ D_{gt}\right\} _{t}$ AR(1) Markov chain with parameter $\rho$, and across firms, $\left\{ D_{gt}\right\} _{g}$ are independent copies of this AR(1) process. Such a process does not permit the CHS representation. Suppose to a contradiction that $D_{gt}$ admits a representation where $D_{gt}=f\left(\alpha_{g},\gamma_{t},\varepsilon_{gt}\right)$ for iid $\alpha_{g},\varepsilon_{gt}$, and strictly stationary $\gamma_t$. In the DGP, $\gamma_{t}=0$, because if $\gamma_{t}$ is stochastic, then there must be dependence across $g$ units, which is ruled out in this DGP. With the resulting representation $D_{gt}=f\left(\alpha_{g},0,\varepsilon_{gt}\right)$, it follows that $Cov(D_{gt}, D_{gs_1}) = Cov(D_{gt}, D_{gs_2})$ for $t \ne s_1 \ne s_2 \ne t$. However, this is untrue in the DGP as $Cov(D_{g,t}, D_{g,t-1}) = \rho/(1-\rho^2) \ne \rho^2/(1-\rho^2) =Cov(D_{g,t}, D_{g,t-2})$.

Main Theoretical Results

This section first presents a central limit theorem (CLT) for a sum of dependent random vectors in this $\psi$-dependent panel setting. After showing that a plug-in estimator can be anticonservative, a conservative variance estimator is proposed and shown to converge to its target estimand, and this target estimand is shown to be asymptotically conservative.

The results are stated in a setting where $\sum_{i\in N_{n}}E\left[Y_{n,i}\right]=0$, and we are interested in constructing the variance of $\sum_{i\in N_{n}}Y_{n,i}$. In light of the central limit theorem, the results can be extended more generally to $\sum_{i\in N_{n}}E\left[Y_{n,i}\right]=\mu_{n}$ for some sequence of means $\mu_{n}$. Let $\lambda_{\min} (V)$ denote the minimum eigenvalue of matrix $V$. With $V_{true} := Var\left(\sum_{i\in N_{n}}Y_{n,i}\right)$, let $\lambda_{n} := \lambda_{\min} (V_{true})$ denote its smallest eigenvalue.

assumptionThe triangular array $\left\{ Y_{n,i}\right\} $ is $\psi$-dependent with the dependence coefficients $\left\{ \theta_{n}\right\} $ satisfying the following conditions. (a) For some constant $C>0$, $\psi_{a,b}\left(f,g\right)\leq C\times ab\left(\left\Vert f\right\Vert _{\infty}+Lip(f)\right)\left(\left\Vert g\right\Vert _{\infty}+Lip(g)\right)$. (b) $\sup_{n\geq1}\max_{s\geq1}\theta_{n,s}<\infty$. (c) For some $p>4$, $\sup_{n\geq1}\max_{i\in N_{n}}\left\Vert Y_{n,i}\right\Vert _{p}<\infty$ a.s.
assumptionThere exists a positive sequence $m_{n}\rightarrow\infty$ such that for $k=1,2$, and the $p$ in (ref), \begin{enumerate} [label=(\alph*)] • $\frac{n}{\lambda_{n}^{1+k/2}}\sum_{s\geq0}c_{n}\left(s,m_{n};k\right)\theta_{n,s}^{1-\frac{2+k}{p}} \to 0$; and • $\frac{n^{2}\theta_{n,m_{n}}^{1-(1/p)}}{\lambda_{n}^{1/2}} \to 0$. \end{enumerate}

These assumptions are used to show the CLT by balancing the dependence strength $\theta_{n,s}$, how fast neighborhoods grow and overlap $c_{n}\left(s,m_{n};k\right)$, and the variance growth. After truncating dependence beyond distance $m_{n}$, we want to ensure that the cumulative effect of all dependent interactions is negligible relative to the variance such that it does not distort the Gaussian limit. These conditions mimic the conditions in KMS Theorem 3.2 to apply their CLT.

In an AR(1) process with positive coefficient $\rho$, $\lambda_n \asymp T$, so (ref)(a) is implied by $T^{-k/2} \sum_{s\geq 0} (m_n-s) \rho^{1-(2+k)/p} \to 0$. By observing that $T^{-k/2} \sum_{s\geq 0} (m_n-s) \rho^{1-(2+k)/p} = T^{-k/2} m_n O(1) - T^{-k/2} O(1)$, the condition is satisfied for, say, $m_n =T^{1/3}$ as $T \rightarrow \infty $. (ref)(b) is implied by $T^{3/2} (\rho^{1-(1/p)})^{m_n} \to 0$, which similarly holds when $m_n$ diverges with $T$.

In the running example of a balanced panel with common random effects within clusters, and an AR(1) process in the common time effects with coefficient $\rho$, $\lambda_n$ is of order $TG^2 \asymp T^2 G$, so (ref)(a) is implied by $(TG^2)^{-k/2} \sum_{s\geq 0} (m_n-s) G^k \rho^{1-(2+k)/p} \to 0$, and (ref)(b) is implied by $(TG)^{2} (TG^2)^{-1/2} (\rho^{1-(1/p)})^{m_n} = GT^{3/2} (\rho^{1-(1/p)})^{m_n} \asymp T^{5/2} (\rho^{1-(1/p)})^{m_n} \to 0$. As before, both conditions may be satisfied with $m_n =T^{1/3}$ as $T \rightarrow \infty $.

theoremSuppose Assumptions (ref) and (ref) hold a.s. Let $\nu$ be a nonstochastic vector with $\dim\left(\nu\right)=\dim\left(Y_{n,i}\right)$ and $\left\Vert \nu\right\Vert =1$. Then, for $S_{n}:=\sum_{i\in N_{n}}\nu^{\prime}\left(Y_{n,i}-E\left[Y_{n,i}\right]\right)$ and $\sigma_{n}^{2}:=Var\left(\nu^{\prime}S_{n}\right)$, \begin{equation} \sup_{t\in\mathbb{R}}\left| \Pr\left\{ \frac{S_{n}}{\sigma_{n}}\leq t\right\} -\Phi(t)\right|\xrightarrow{a.s.}0, \end{equation} where $\Phi$ denotes the distribution function of $N(0,1)$.

This result suffices for a uniform multivariate central limit theorem due to the Cramer-Wold device.

The CHS plug-in variance estimator designed under homogeneous means that accounts for cluster dependence and serial correlation is anticonservative when means are heterogeneous. Using $y_{t}:=\sum_{i\in N_{t}^{T}}Y_{n,i}$, the CHS variance estimator is:

align*[align* omitted — 684 chars of source]

where $\omega\left(m,M\right)$ is a kernel function at distance $m$ with bandwidth $M$. The first line is exactly the CGM estimator, the second line accounts for serial correlation across clusters, and the final line adjusts for double counting. Notably, CHS omit the final line of the plug-in estimator in their implementation, because it is asymptotically negligible in their framework. In contrast, the true variance is:

align[align omitted — 629 chars of source]
exampleSuppose we have one observation in each $\left(g,t\right)$ intersection, and we have $m=1$ dependence. We can compare the oracle CHS variance estimand (where $\omega(m,M)=1$ for all $m\leq 1$ and 0 otherwise) with the true variance. As before, this example can be extended to $T \rightarrow \infty$ and $G \rightarrow \infty$. The scalar $E\left[Y_{n,i}\right]$ values are in (ref). \begin{table}[H] \caption{Values of scalar $E[Y_{n,i}]$} \begin{tabular}{|c|c|c|c|c|} \hline & $t=1$ & $t=2$ & $t=3$ & $t=4$ \tabularnewline \hline $g=1$ & $-1$ & $-1$ & $1$ & $1$ \tabularnewline \hline $g=2$ & $1$ & $-1$ & $-1$ & $1$ \tabularnewline \hline \end{tabular} \end{table} Then, observe that $\sum_{i \in N_n}\sum_{j\in N_{g(i)}^{G}}E\left[Y_{n,i}\right]E\left[Y_{n,j}\right]=0$, $\sum_{i \in N_n}\sum_{j\in N_{t(i)}^{T}}E\left[Y_{n,i}\right]E\left[Y_{n,j}\right]=4$, and $\sum_{i \in N_n}\sum_{j\in N_{t(i),g(i)}^{T\cap G}}E\left[Y_{n,i}\right]E\left[Y_{n,j}\right]=8$. In the serially correlated components, $\sum_{t=1}^{T-m}E\left[y_{t}\right]E\left[y_{t+m}^{\prime}\right]=-3$ and $\sum_{t=1}^{T-m}\sum_{g=1}^{G}\sum_{i\in N_{t,g}^{T\cap G}}\sum_{j\in N_{t+m,g}^{T\cap G}}E\left[Y_{n,i}\right]E\left[Y_{n,j}\right]=-3$. Then, on aggregate, using the same argument as Example 1, \begin{align*} & \sum_{i \in N_n}\sum_{j\in N_{g(i)}^{G}}E\left[Y_{n,i}\right]E\left[Y_{n,j}\right] + \sum_{i \in N_n}\sum_{j\in N_{t(i)}^{T}}E\left[Y_{n,i}\right]E\left[Y_{n,j}\right] - \sum_{i \in N_n}\sum_{j\in N_{t(i),g(i)}^{T\cap G}}E\left[Y_{n,i}\right]E\left[Y_{n,j}\right] \\ & +2\sum_{t=1}^{T-m}E\left[y_{t}\right]E\left[y_{t+m}\right] - 2\sum_{t=1}^{T-m}\sum_{g=1}^{G}\sum_{i\in N_{t,g}^{T\cap G}}\sum_{j\in N_{t+m,g}^{T\cap G}}E\left[Y_{n,i}\right]E\left[Y_{n,j}\right] = -4, \tag{\stepcounter{equation}\theequation} \end{align*} so the standard CHS estimand is anticonservative. If we further ignore the $\sum_{t=1}^{T-m}\sum_{g=1}^{G}\sum_{i\in N_{t,g}^{T\cap G}}\sum_{j\in N_{t+m,g}^{T\cap G}}Y_{n,i}Y_{n,j}^{\prime}$ component due to asymptotic negligibility in the covariance component in the CHS framework, the underestimation problem of the CHS estimand would be even more severe.

To restore validity, the proposed conservative variance estimator and estimand are as follows:

align*[align* omitted — 373 chars of source]
align*[align* omitted — 437 chars of source]

Since $\omega(.)$ is a kernel function, $\omega(.)$ satisfies $\left|\omega(.)\right|\leq1$. For implementation, the $M$ is selected using the bandwidth selection procedure of andrews1991heteroskedasticity and $\omega(.)$ is the triangular kernel, following CHS. Here, $\lambda_{n}$ is analogous to $n$ in the variance estimation framework of KMS, but $\lambda_{n}$ can diverge at a rate faster than $n$ depending on the correlation structure. The balanced panel is an example where $\lambda_n$ diverges faster than $n$.

assumption$\frac{n}{\lambda_{n}^{2}}\sum_{s\geq0}c_{n}\left(s,M;2\right)\theta_{n,s}^{1-\frac{4}{p}} \to 0$.

Observe that (ref) is different from (ref)(a), where $k=2$ and $m_{n}=M$ uses the chosen bandwidth, which should diverge with $n$. This assumption mimics Assumption 4.1(iii) in KMS. The sufficient conditions for time series and the balanced panel are the same as (ref)(a). We no longer require an analog of (ref)(b) because the consistency result already ignores information beyond the bandwidth $M$.

theoremUnder Assumptions (ref) and (ref), and $\lambda_{\min} (V_{con})/ \lambda_n \geq 1+ o(1)$, the variance estimator is consistent: \begin{equation} V_{con}^{-1}\hat{V}_{con}\xrightarrow{p}I_{K}. \end{equation}

Unlike the plug-in variance estimator, this variance estimator is conservative and hence controls size. Instead of using the full variance, I compare the kernel adjusted variance:

align*[align* omitted — 482 chars of source]

Similar to CHS, we neglect the final $\sum_{t=1}^{T-m}\sum_{g=1}^{G}\sum_{i\in N_{t,g}^{T\cap G}}\sum_{j\in N_{t+m,g}^{T\cap G}}Cov\left(Y_{n,i},Y_{n,j}\right)$ term as it is asymptotically negligible. The propositions that follow claim that $V_{con}$ is conservative relative to $V_{adj}$, and $V_{adj}$ is asymptotically the same as $V_{true}$, so $V_{con}$ is asymptotically conservative relative to $V_{true}$.

proposition$V_{con}-V_{adj}$ is positive semidefinite.
assumptionFor arbitrary positive and finite constants $\vartheta_1, \vartheta_2, \vartheta_3 $, the following hold almost surely: \begin{enumerate} [label=(\alph*)] • $\frac{\vartheta_{1}n}{\lambda_{n}}\sum_{m=1}^{M}\left|\omega\left(m,M\right)-1\right|\delta_{n}^{\partial}\left(m;1\right)\theta_{n,m}^{1-\frac{2}{p}}\rightarrow0$; • $\frac{\vartheta_{2}}{\lambda_{n}}\sum_{m\geq1}\sum_{t=1}^{T-m}\sum_{g=1}^{G}\left|N_{t,g}^{T\cap G}\right|\left|N_{t+m,g}^{T\cap G}\right|\theta_{n,m}^{1-\frac{2}{p}}\rightarrow0$; and • $\frac{\vartheta_{3}}{\lambda_{n}}\sum_{m=M+1}^{T-1}\sum_{t=1}^{T-m}\left|N_{t}^{T}\right|\left|N_{t+m}^{T}\right|\theta_{n,m}^{1-\frac{2}{p}}\rightarrow0$. \end{enumerate}

This assumption is used to control the bias introduced by kernel weighting and truncation in variance estimation. The first of these assumptions is analogous to KMS Assumption 4.1(ii), ensuring that the kernel weights are well-behaved. For instance, using a uniform kernel of $\omega(m,M) = 1$ trivially satisfies this assumption. (ref)(b) ensures that the terms in the final line of (ref) omitted by CHS are indeed asymptotically negligible. (ref)(c) ensures that the dependence beyond the bandwidth is negligible for $V_{adj}$ to accurately converge to $V_{true}$.

In the running AR(1) time series example with $\rho \in (0,1)$, when using the triangular kernel $\omega(m,M) = 1 - \frac{m}{M+1}$, (ref)(a) is implied by $\sum_{m=1}^{M}\left|\omega\left(m,M\right)-1\right|\delta_{n}^{\partial}\left(m;1\right) \allowbreak \left( \rho^{1-\frac{2}{p}} \right)^m \to 0$. With the triangular kernel, $\sum_{m=1}^{M}\left|\omega\left(m,M\right)-1\right|\delta_{n}^{\partial}\left(m;1\right) \left( \rho^{1-\frac{2}{p}} \right)^m \asymp \sum_{m=1}^{M}\frac{m}{M+1} \left( \rho^{1-\frac{2}{p}} \right)^m = \frac{1}{M+1} O(1) = o(1)$. (ref)(b) is not required, because the last line of (ref) no longer features in the variance expression as an adjustment term. (ref)(c) is implied by $T^{-1} \sum_{m=M+1}^{T-1}\sum_{t=1}^{T-m} \left(\rho^{1-\frac{2}{p}}\right)^m \rightarrow0$. Observe that $T^{-1} \sum_{m=M+1}^{T-1}\sum_{t=1}^{T-m} \left(\rho^{1-\frac{2}{p}}\right)^m = T^{-1} \sum_{m=M+1}^{T-1} (T-m) \left(\rho^{1-\frac{2}{p}}\right)^m \leq \sum_{m=M+1}^{T-1} \left(\rho^{1-\frac{2}{p}}\right)^m \leq \left(\rho^{1-\frac{2}{p}}\right)^{M+1} / \left( 1- \rho^{1-\frac{2}{p}} \right)$, which converges to zero as $M\rightarrow \infty$, so (ref)(c) is satisfied.

In the running example of a balanced panel, (ref)(a) is implied by $\frac{1}{G} \sum_{m=1}^{M}\left|\omega\left(m,M\right)-1\right|\delta_{n}^{\partial}\left(m;1\right) \left( \rho^{1-\frac{2}{p}} \right)^m \to 0$. Since $\delta_{n}^{\partial}\left(m;1\right) \asymp G$, the argument from the time series exposition applies. (ref)(b) is implied by $(G^2T)^{-1} \sum_{m \geq 1} \sum_{t=1}^{T-m} \allowbreak G \left( \rho^{1-\frac{2}{p}} \right)^m \to 0$. Then, observe that $(GT)^{-1} \sum_{m \geq 1} \sum_{t=1}^{T-m} \left( \rho^{1-\frac{2}{p}} \right)^m = G^{-1} O(1) = o(1) $. (ref)(c) is implied by $(G^2T)^{-1} \sum_{m=M+1}^{T-1}\sum_{t=1}^{T-m} G^2 \left(\rho^{1-\frac{2}{p}}\right)^m \rightarrow0$, so the same argument as the time series example can be applied for (ref)(c) to be satisfied.

In the case where observations are independent, Assumptions (ref), (ref), and (ref) are automatically satisfied as $\theta_{n,s} =0$ for all $s$, so the conclusions can be applied. Consequently, the variance estimators are still consistent and we do not necessarily have a slower convergence rate by allowing for dependence.

propositionUnder Assumptions (ref) and (ref), $V_{true}^{-1}V_{adj}\rightarrow I_{K}$.

Under the conditions of (ref), the condition that $\lambda_{\min} (V_{con})/ \lambda_n \geq 1+ o(1)$ in (ref) is automatically satisfied, as $\lambda_{\min} (V_{con}) / \lambda_n= \frac{\lambda_{\min} (V_{con})}{\lambda_{\min} (V_{adj})} \frac{\lambda_{\min} (V_{adj})}{ \lambda_n}$, with $\lambda_{\min} (V_{con}) / \lambda_{\min} (V_{adj}) \geq 1 $ due to (ref) and $\lambda_{\min} (V_{adj}) / \lambda_{n} = 1 + o(1)$.

exampleFollowing the extent of conservativeness in time series characterized in Example 1, we may also characterize how conservative the variance estimator is with two-way cluster dependence (without cross-cluster dependence). Suppose we consider the framework of davezies2021empirical, where $Y_{n,i}=f\left(\alpha_{g},\gamma_{t},\varepsilon_{gt}\right)$, for $\alpha_{i},\gamma_{t},\varepsilon_{gt}$ iid with variances $\sigma_{\alpha}^{2}$, $\sigma_{\gamma}^{2}$, and $\sigma_{\varepsilon}^{2}$ respectively, and there is only one observation in each $\left(g,t\right)$ intersection. Following xu2024clustering, CGM2 (which is a special case of the proposal here) would overestimate the true variance by: \begin{equation} \sum_{g,t}\sum_{i,j\in N_{t,g}^{T\cap G}}E\left[Y_{n,i}Y_{n,j}\right]=\sum_{g,t}\sum_{i\in N_{t,g}^{T\cap G}}E\left[Y_{n,i}^{2}\right]=\sum_{i\in N_{n}}Var\left(Y_{n,i}\right), \end{equation} even with constant means. This result is similar to the result that we have observed for time series. Then, we can evaluate the true variance $Var\left(\sum_{i \in N_n}Y_{n,i}\right)$ when $Y_{n,i}=\alpha_{g(i)}+\gamma_{t(i)}+\varepsilon_{i}$. Namely, \begin{align*} V_{true} & =\sum_{g=1}^{G}\sum_{i\in N_{g}^{G}}\sum_{j\in N_{g}^{G}}Cov\left(Y_{n,i},Y_{n,j}\right)+\sum_{t=1}^{T}\sum_{i\in N_{t}^{T}}\sum_{j\in N_{t}^{T}}Cov\left(Y_{n,i},Y_{n,j}\right) -\sum_{i\in N_{n}}Var\left(Y_{n,i}\right)\\ & =\sum_{g=1}^{G}\left|N_{g}^{G}\right|^{2}\sigma_{\alpha}^{2}+\sum_{t=1}^{T}\left|N_{t}^{T}\right|^{2}\sigma_{\gamma}^{2}-\sum_{i\in N_{n}}\left(\sigma_{\alpha}^{2}+\sigma_{\gamma}^{2}+\sigma_{\varepsilon}^{2}\right). \tag{\stepcounter{equation}\theequation} \end{align*} If we further assume $\sigma_{\alpha}^{2}=\sigma_{\gamma}^{2}=\sigma_{\varepsilon}^{2}=1$ and that there are $T$ clusters and $T$ periods such that $n=T^{2}$, then we can simplify the expression even further $\sigma_{n}^{2}=2T^{3}-3T^{2}$. Then, as $T\rightarrow\infty$, since $Var(Y_{n,i}) = O(1)$, we have $\left(\sum_{i\in N_{n}}Var\left(Y_{n,i}\right)\right)/\sigma_{n}^{2}\asymp1/T=o(1)$. Consequently, we observe that, just like the time series case, having stronger dependence would make the robust variance estimator less conservative.

Although the proposed variance estimator is conservative, it nonetheless shows that validity can be restored with heterogeneous means. The challenge for future research is then to shrink the conservative variance estimators towards the true variance. Some progress has already been made in particular situations (e.g., abadie2020sampling, abadie2023should, xu2024clustering) when there is more information or structure in the data, but there is certainly room for further innovation.

While the theoretical exposition focuses on means, the extension to the method of moments (which includes standard regressions) is straightforward. With a regression equation $Y_i = X_i^\prime \beta + u_i$, and array $\{ (X_i^\prime, u_i) \}$ satisfying $\psi$-dependence, as long as $\frac{1}{n} \sum_{i \in N_n} X_i X_i^\prime$ has a minimum eigenvalue bounded away from zero, and $\{ X_i u_i \}$ and $\{ vec(X_i X_i^\prime) \}$ satisfy the conditions of the (ref), the OLS estimator $\hat{\beta}$ is asymptotically normal and consistent. To calculate the variance of the coefficient estimator, we can use the predicted residuals to construct the conservative variance estimator of $\sum_{i \in N_n} X_i u_i$, so when the conditions for (ref) and (ref) hold for $X_i u_i$ in place of $Y_{n,i}$, exploiting the consistency of $\hat{\beta}$, the conservative variance estimator will converge to its estimand, which is conservative for $Var\left( \sum_{i \in N_n} X_i u_i\right)$. Standard OLS arguments then apply. In light of this argument, the next section provides numerical illustrations based on the regression coefficient.

Numerical Illustrations

The numerical illustrations follow CHS to show how the variance estimation method that is robust to heterogeneous means compares with existing methods.

Simulation

Data is generated based on the linear model

equation[equation omitted — 79 chars of source]

where $\beta^h_{gt}$ is a heterogeneity term, and the remaining variables $(X_{gt}, U_{gt})$ mimic the setup from CHS, where

align*[align* omitted — 259 chars of source]

I set $(\beta_0, \beta_1) = (0.1, 0.1)$, and $(w_\alpha, w_\gamma, w_\varepsilon) = (0.15, 0.20, 0.15)$. The latent components $(\alpha^x_g, \alpha^u_g, \varepsilon^x_{gt}, \varepsilon^u_{gt})$ are all mutually independent $N(0,1)$. The latent common time effects $(\gamma^x_t, \gamma^u_t)$ are dynamically generated according to an AR(1) process with coefficient $\rho$. Finally, I set $\beta^h_{gt}$ to be an alternating sequence of $.1$ and $-.1$.

In each simulation draw, I regress $Y_{gt}$ on $X_{gt}$ and a constant, and estimate $\hat{\beta_1}$. Various methods are used to calculate the standard error for $\hat{\beta}_1$. The corresponding confidence interval is then constructed and I calculate if each method rejects the null value of $\beta_1 = 0.1$ to evaluate its ability to control size. The results are reported in (ref).

table[table omitted — 1,212 chars of source]

From (ref), it is evident that over-rejection of the null can be very severe if we do not use the variance estimator that is robust to heterogeneous means when we introduce heterogeneity due to $\beta^h_{gt}$. Even though the variance estimators that account for dependence inflate the standard errors to lower over-rejection, they are still insufficient compared to the proposal in this paper. When the correlation over time is sufficiently high, the test achieves a rejection rate that is close to the nominal rate, showing how there are designs where this procedure is not too conservative.

Empirical Application

The new variance estimator is applied to a panel of industry portfolios across time. We may want to accommodate heterogeneous means in this context when portfolio labels and time labels are inherently meaningful and are not exchangeable. The Fama-French three-factor model can be written as:

equation[equation omitted — 111 chars of source]

where $R_{gt} - R_{ft}$ is the excess return as $R_{gt}$ denotes the total return of portfolio $g$ in month $t$, and $R_{ft}$ denotes the risk-free rate of return. $R_{Mt}$ is the total market portfolio return, $SMB_t$ is the size premium (small minus big), $HML_t$ is the value premium (high minus low book to market ratio), and $\beta$ are factor coefficients.

I use the industry portfolio dataset from CHS that contains 44 industry portfolios across 119 months. To mimic the analysis in CHS, I use a within-transformation for both $Y$ and $X = (R_{M} - R_{f}, SMB, HML)$ to remove additive fixed effects. Then, I estimate the $\beta$ coefficients and calculate the standard errors using various methods. With $\tilde{X}$ denoting the transformed random variable, inference on coefficients relies on a CLT applied to $\sum_{g,t} \tilde{X}_{gt} \tilde{e}_{gt}$. Allowing for mean heterogeneity is reasonable as the residuals of returns likely have heterogeneous means over portfolios and time.

table[table omitted — 820 chars of source]

From (ref), standard errors estimated using HM are higher than the standard errors estimated using other methods. However, the standard errors are not unreasonably large, as we can retain the statistical significance of the coefficient on HML as other methods do at the 5% level, while the statistical significance of SMB can now be called into question, as CHS already show. A further takeaway from the results is that cross-cluster serial correlation is empirically important: if the serial correlation is indeed negligible, we would expect the CHS standard errors to be fairly similar to CGM standard errors, which is not the case for two of the coefficients.