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.
59,488 characters · 12 sections · 43 citation commands
Subsampling Under Two-way Clustering with Serial Correlation
Inference on clustering data has always been an important subject and has been extensively studied in the past. In the realm of multiway clustering, cameron_robust_2011 is probably one of the most cited papers, which proposes multiway-clustering-robust standard errors. petersen_estimating_2009 explains when firm clustering, time clustering, or two-way clustering is appropriate. mackinnon_wild_2021 studies the wild bootstrap procedures and theoretical conditions for valid inference with two-way clustered data; menzel_bootstrap_2021 also proposes using bootstrap method for inference under multiway clustering and gives the corresponding asymptotic theory. However, none of the above studies allows for correlation within at least one dimension of the clusters, which may accompany the multiway clustering problem up to some extent. The multiway clustering is quite common in the panel data, where one dimension is usually firm or individual and the other one is time. Yet, the common time effect is typically serially correlated. Hence, this problem poses a new challenge to the inference method, which requires the researchers to account for both multiway clustering and the potential serial correlation at the same time, as bertrand_how_2004 has shown the necessity of considering dependence in practical study. chiang_standard_2024 is probably one of the first papers that takes such problem into consideration. They give a new variance estimator and the corresponding asymptotic theory under a two-way clustering framework, with allowing the serial correlation up to a certain threshold. chen_fixed-b_2024 obtains a new asymptotic result of the chiang_standard_2024 estimator (CHS estimator hereafter) under a fixed ratio between the bandwidth and the sample size, which they call it fixed-b, and proposes two bias-corrected variants of the CHS estimator along with the critical values under the fixed-b asymptotics. The most recent study on this issue is probably hounyo_reliable_2024, which proposes two multiway wild cluster bootstrap methods based on CHS asymptotic-based approach. They show the wild bootstrap method works both theoretically and practically under a two-way clustering with serial correlation framework.\\ Interestingly, as opposed to bootstrap, using subsampling for inference under two-way clustering with serial correlation remains almost unexplored. The use of subsampling as an inference method can be dated back to mahalanobis_sample_1946, which suggests to employ subsamples to estimate variance in crop yields. carlstein_use_1986 uses non-overlapping subsamples to estimate variance in time series. The leave-many-out jackknife, which is very similar to subsampling, is studied as an approach to estimate variance by wu_jackknife_1986 and shao_general_1989. politis_subsampling_1999 presents a very complete theory on subsampling, including basic properties of subsampling distribution under i.i.d. case, stationary and non-stationary time series, and various other settings. As subsampling works under both independent data and time series, then it seems promising that it would work under two-way clustering with serial correlation. menzel_bootstrap_2021 considers using subsampling for inference under two-way clustering where both individual and time effects are i.i.d. draws. Other recent usage of subsampling as an inference method include chernozhukov_inference_2011 and kurisu_subsampling_2025, both of which consider extremal conditional quantile regression. Nevertheless, there is no similar result or theorem that works for the two-way clustering with serial correlation, and this paper aims to fill this gap. Using subsampling for this problem is intuitive, as each proper subsample is in fact a smaller sample and would reflect the true underlying distribution. Then, repeatedly drawing subsamples from the original sample is similar to drawing many samples from the true underlying distribution, which could not be done in reality. Therefore, we can recover some statistics of the true underlying distribution via this procedure.\\ The main contribution of this paper is that it shows the validity of using subsampling methods for inference under two-way clustering with serial correlation. We prove the subsampling distribution converges to the true underlying distribution in probability pointwise under mild conditions; henceforth, its quantiles also converge to their counterparts in the true underlying distribution and can be used to construct confidence intervals with desired levels. Furthermore, such method does not require knowing the form of the limiting distribution but only its existence. The existing literature on this problem is either CLT-based (chiang_standard_2024; chen_fixed-b_2024) or depends on the normality of t-statistics (hounyo_reliable_2024). However, due to a lack of subsample size selection algorithm, using quantiles of the subsample distribution for inference is limited at this stage. With an additional uniform integrability condition, we also develop a variance estimator based on the subsampling distribution, which is consistent to the variance of the true underlying distribution. Such variance estimator can be used to construct confidence intervals or test statistics as well. We further suggest a data-driven method to select subsample size by connecting this subsampling variance estimator to the well-studied spectral density estimator, when the parameter of interest is the population average. A bias-corrected version of this variance estimator that reduces the order of bias is also proposed. The simulation results agree with the theoretical argument and show our methods are comparable, in terms of the inference quality, to the existing leading method, except when the serial correlation is extremely strong.\\ The paper is organized as the following. Section (ref) talks about the subsampling distribution, and Section (ref) is about the variance estimator. The simulation result is presented in Section (ref). Appendix (ref) and (ref) give proofs to the main theorems, whereas Appendix (ref) contains proofs of technical lemmas. The derivations of some results and equalities are in Appendix (ref), along with other technical details.
We consider a panel data $\{X_{nt}:1\leq n\leq N,1\leq t\leq T\}$, where $N$ is the sample size of individuals and $T$ the sample size of time periods. We assume each $X_{nt}$ is generated via some Borel-measurable function $f$
with $\{\alpha_{n}\},\{\gamma_{t}\}$, and $\{\varepsilon_{nt}\}$ being mutually independent, $\{\alpha_{n}\}$ being i.i.d. across $n$, $\{\gamma_{t}\}$ being strictly stationary and serially correlated, and $\{\varepsilon_{nt}\}$ being i.i.d. across $(n,t)$.
Let $\theta$ be the parameter of interest. In the main body of this paper, we assume $\theta$ is real-valued, while all the theorems and results can be easily extended to the multivariate case. Define $\hat{\theta}_{NT}=\hat{\theta}_{NT}(\{X_{nt}:1\leq n\leq N,1\leq t\leq T\})$, an estimator of $\theta$ using the full sample. Furthermore, let $\tau_{NT}$ be some normalizing constant and $J_{NT}(x):=\mathbb{P}(\tau_{NT}(\hat{\theta}_{NT}-\theta)\leq x)$ to be the cumulative distribution function of the root $\tau_{NT}(\hat{\theta}_{NT}-\theta)$. Since our goal is to do inference properly on the estimator $\hat{\theta}_{NT}$, then we also need the limiting distribution of the root to exist, which is described by the following assumption.
Subsampling in panel data is not merely a combination of subsampling in i.i.d. data and the one in time series data. In i.i.d. data, we usually pick $b$ indices from $\{1,\dots,N\}$ without replacement and allow overlapping subsamples, and we are able to bound tail probabilities by Hoeffding-Serfling inequality (serfling_u-statistics_1980, Theorem A, p.201) due to independence. However, under two-way clustering with serial correlation, none of any pair of observations are independent. Yet, any two observations that do not share the same individual index and far away in time are close to independent due to the strong mixing condition on the time component $\{\gamma_{t}\}$. Therefore, to construct a subsample, we randomly partition $\{1,\dots,N\}$ into $\lceil\frac{N}{b}\rceil$ disjoint subsets with all but one having a size of $b$ and pick $l$ consecutive indices from $\{1,\dots,T\}$ to preserve the serial correlation. For instance, $\{X_{nt}:n\in I_{i},k\leq t\leq k+l-1\}$ is a subsample of individuals from $I_{i}\subset\{1,\dots,N\}$, with $|I_{i}|=b$, and starting from time $k$. In total, there are $\lceil\frac{N}{b}\rceil\cdot(T-l+1)$ possible subsamples. The corresponding subsample estimator $\hat{\theta}_{b,l,i,k}$ is defined as $\hat{\theta}_{b,l,i,k}=\hat{\theta}_{bl}(\{X_{nt}:n\in I_{i},k\leq t\leq k+l-1\})$. The merit of the subsampling is a “tautology": a subsample is in fact a smaller sample. The distribution of $\tau_{bl}(\hat{\theta}_{bl}-\theta)$ is exactly $J_{bl}$, and it has the same limiting distribution as the root $\tau_{NT}(\hat{\theta}_{NT}-\theta)$ as long as we think $b$ and $l$ as the smaller versions of $N$ and $T$ respectively. Not only both $b$ and $l$ must go to infinity, but also their relation has to be the same as the relation between $N$ and $T$ in the limit. For example, if $N$ goes to infinity faster than $T$, then $b\ll l$ would result in $J_{bl}$ and $J_{NT}$ having different limits. Another way to think this is to consider two sequences $N({m}):\mathbb{N}\rightarrow\mathbb{N}$ and $T({m}):\mathbb{N}\rightarrow\mathbb{N}$. $N$ and $b$ are both from the sequence $(N(m))$, and $T$ and $l$ are both from the sequence $(T(m))$. But for the $m_{1}$ and $m_{2}$ such that $N=N(m_{1})$, $T=T(m_{1})$ and $b=N(m_{2})$, $l=T(m_{2})$, $m_{1}>m_{2}$. With these, repeatedly drawing subsamples from the distribution $J_{bl}$ and evaluating the empirical distribution should somehow reflect the true underlying distribution $J_{NT}$ which eventually goes to $J$, the limit associated to our parameter of interest $\theta$. These lead to our first main theorem.\\ Define the subsampling distribution as
where $N_{b}=\lceil\frac{N}{b}\rceil$ and $q=T-l+1$.
The proof is in the Appendix (ref). Part (a) provides a basic pointwise convergence in probability result which is used to derive the convergence of quantile in part (c). In part (b), the pointwise convergence in probability is strengthened to an uniform convergence in probability if the limiting distribution is continuous. Part (c) provides a valid inference method using subsampling: we can use the quantile from the empirical distribution $L_{N,T,b,l}(x)$ to construct a confidence interval, of which coverage probability goes to the nominal level. The one-sided confidence interval $(\hat{\theta}_{NT}-\frac{c_{b,l}^{L}(1-\alpha)}{\tau_{NT}},\infty)$ has an asymptotic coverage probability being the nominal level $1-\alpha$. To construct a two-sided equal-tail confidence interval with a correct asymptotic coverage probability, one should use $(\hat{\theta}_{NT}-\frac{c_{b,l}^{L}(1-\frac{\alpha}{2})}{\tau_{NT}},\hat{\theta}_{NT}-\frac{c_{b,l}^{L}(\frac{\alpha}{2})}{\tau_{NT}})$, where left and right endpoints are interchanged like the bootstrap percentile-t interval. This subsampling quantile method does not require any specific form of the limiting distribution, whereas the existing literature on two-way clustering with serial correlation are built on the Gaussianity of the asymptotic $J$. This allows us to apply subsampling to some non-standard case.
Although the knowledge of the exact form of $J$ is not required, knowing the normalizing constant $\tau_{NT}$ is still necessary to construct valid confidence interval. That is, such method works in the degenerate case when it is known to be degenerate.
In some cases, we can also do inference on $\hat{\theta}_{NT}$ properly by estimating its asymptotic variance. For example, if the limiting distribution is Gaussian and has an unknown variance $V$, then $\hat{\theta}_{NT}$ approximately follows the distribution of $N(\theta,\frac{V}{\tau_{NT}^{2}})$. Furthermore, whenever we want to evaluate the efficiency or derive the mean squared-error of the estimator $\hat{\theta}_{NT}$, a variance estimator is usually required. We here provide a subsampling method to estimate the variance consistently by using the variance associate to the empirical distribution $L_{N,T,b,l}$ defined in ((ref)).
The proof is in Appendix (ref). Note that $\hat{\sigma}_{N,T,b,l}^{2}$ estimates $\operatorname*{Var}(\tau_{NT}\hat{\theta}_{NT})$. To estimate the variance of $\hat{\theta}_{NT}$, we simply divide $\hat{\sigma}_{N,T,b,l}^{2}$ by $\tau_{NT}^{2}$. The condition $\frac{\tau_{bl}}{\tau_{NT}}$, which lies in Assumption (ref), is in fact redundant. Throughout the proof, the only two requirements put on the normalizing sequence $\{\tau_{NT}$ are $\operatorname*{Var}(\tau_{NT}\hat{\theta}_{NT})$ converging to some finite $V$ and the fourth order uniform integrability on $\{\tau_{NT}(\hat{\theta}_{NT}-\mathbb{E}[\hat{\theta}_{NT}])\}$. Hence, if $\hat{\theta}_{NT}$ is consistent to some deterministic $\theta$, then $\hat{\sigma}_{N,T,b,l}^{2}\rightarrow_{p}0$. Although Theorem (ref) does not explicitly rely on Assumption (ref), knowing which family of distribution $\lim\tau_{NT}(\hat{\theta}_{NT}-\theta)$ belongs to is required if we want to construct confidence intervals using this variance estimator.
Up to this point, asymptotic theory only requires $b=o(N)$ and $l=o(T)$ but does not tell us which exact number to pick. Selecting the subsample size is quite important as it would largely affect the result of our variance estimator just as the bandwidth selection in many nonparametric estimation problems. Traditionally, researchers used some sort of mean squared-error to determine the optimal bandwidth choice. In this paper, we aim to minimize the asymptotic mean squared-error (AMSE) of $\hat{\sigma}_{N,T,b,l}^{2}$, which is defined as
where $Bias_{\infty}(\hat{\sigma}_{N,T,b,l}^{2},V)$ and $\operatorname*{Var}_{\infty}(\hat{\sigma}_{N,T,b,l}^{2})$ are defined as the leading terms of $Bias(\hat{\sigma}_{N,T,b,l}^{2},V)$ and $\operatorname*{Var}(\hat{\sigma}_{N,T,b,l}^{2})$ respectively. However, either $Bias(\hat{\sigma}_{N,T,b,l}^{2},V)$ or $\operatorname*{Var}(\hat{\sigma}_{N,T,b,l}^{2})$ depends on the specific form of $\hat{\theta}_{NT}$. Therefore, we focus on mean estimation throughout the rest of this section (and the next section).\\ Under the data generating process of ((ref)) and Assumption (ref), we are interested in $\theta=\mathbb{E}[X_{nt}]$, and we assume $\theta=0$ without loss of generality. Let $\hat{\theta}_{NT}=\frac{1}{NT}\sum_{n=1}^{N}\sum_{t=1}^{T}X_{nt}$ be an estimator to $\theta$. Furthermore, we use a projection technique to transform $X_{nt}$ into a linear combination of some random variables. Define $a_{n}=\mathbb{E}[X_{nt}|\alpha_{n}]$, $b_{t}=\mathbb{E}[X_{nt}|\gamma_{t}]$, and $e_{nt}=X_{nt}-a_{n}-b_{t}$. Then
The random variables $\{a_{n}\}$, $\{b_{t}\}$, and $\{e_{nt}\}$ satisfy the following: (i) $\{a_{n}\}$ is i.i.d., and $\{b_{t}\}$ is strictly stationary; (ii) $\{e_{nt}\}$ is identically distributed across $(n,t)$; (iii) $\mathbb{E}[a_{n}]=\mathbb{E}[b_{t}]=\mathbb{E}[e_{nt}]=0$; (iv) $\{a_{n}\}$ and $\{b_{t}\}$ are independent, and $\{a_{n}\}$, $\{b_{t}\}$, and $\{e_{nt}\}$ are mutually uncorrelated; (v) $e_{nt}$ and $e_{mp}$ are independent conditional on $(\gamma_{t},\gamma_{p})$ for any $n\neq m$ and any $t,p\in\mathbb{N}$. Furthermore, define the follow variance and autocovariances
In chiang_standard_2024, the authors have established the asymptotic theory for the mean estimator $\hat{\theta}_{NT}$ under the framework ((ref)) and the assumptions below.
Then, under Assumption (ref), (ref), and (ref),
where $V=V_{a}+c\sum_{k=-\infty}^{\infty}R_{b}(k)<\infty$ (see the derivation of $V$ in Appendix (ref)). Now, we are ready to introduce a data-driven approach for selecting the subsample size in the context of mean estimation, which has a close relation to the bandwidth choice in the classical spectral density estimation. To see this relationship, note that the limiting variance $V$ can also be written as
where $SP_{b}(\lambda)\equiv\frac{1}{2\pi}\sum_{k=-\infty}^{\infty}R_{b}(k)e^{-ik\lambda}$ is the spectral density function. Meanwhile, as the subsampling variance estimator $\hat{\sigma}_{N,T,b,l}^{2}$ is consistent and has the following expectation
then we have the approximation
Assumption (ref) says in order for $\hat{\sigma}_{N,T,b,l}^{2}$ to be a consistent estimator of V, $\frac{N}{T}-\frac{b}{l}\rightarrow0$. Then, we may assume the ratio between $b$ and $l$ is fixed and $b=\frac{N}{T}\cdot l$. Hence, there is only one parameter to choose and ((ref)) becomes
This expression connects the subsampling variance estimator to the spectrum estimation at frequency zero using the Bartlett window. Not surprisingly, the subsample size $l$ is closely related to the Bartlett window's width; in particular, $l$ is exactly the inverse of the Bartlett window's width if we also treat $\frac{1}{l}\sum_{k=-l+1}^{l-1}(1-\frac{|k|}{l})R_{e}(k)$ negligible. At this stage, we treat this part negligible, while it would be valuable in the future if an algorithm that does not neglect such part is devised. To find the optimal subsample size, We will employ a two-step procedure to choose $l$. First, we try to identify and estimate $V_{b}$, $V_{e}$, $R_{b}$, and $R_{e}$ by some consistent estimators $\hat{R}_{b}$ and $\hat{R}_{e}$ respectively. Then, we plug in those consistent autocovariance estimators and compute $l_{opt}$ minimizing the asymptotic mean squared error. It is important to note that $\hat{R}_{b}$ and $\hat{R}_{e}$, as well as $\hat{V}_{a}$, can also be used to estimate $V$. However, we here only use them to compute the optimal subsample size. \\ If we treat $\frac{\hat{V}_{e}}{l}+\frac{2}{l}\sum_{k=1}^{l-1}(1-\frac{k}{l})\hat{R}_{e}(k)$ negligible, then
and the problem boils down to the optimal bandwidth choice of the Bartlett kernel, which has been widely studied in the past. For instance, andrews_heteroskedasticity_1991 obtained the optimal bandwidths for a class of kernels by using an asymptotic truncated mean squared error criterion; newey_automatic_1994 proposed a plug-in procedure for selecting the bandwidth. In this paper, we adopt an iterative plug-in method developed by buhlmann_locally_1996, which was also used by buhlmann_block_1999 to select the block length of the moving blocks bootstrap in time series. To see how it connects to the optimal bandwidth choice of the Bartlett kernel, recall that $V=V_{a}+c\cdot2\pi SP_{b}(0)$. $V_{a}$ here can be thought as a constant, as a typical estimator $\hat{V}_{a}$ usually does not depend on the bandwidth parameter. Given the limiting ratio $c=\lim\frac{N}{T}$ cannot be inferred anyway, our goal becomes to select an $l$ that minimizes $AMSE(\sum_{k=-l+1}^{l-1}(1-\frac{|k|}{l})\hat{R}_{b}(k),\sum_{k=-\infty}^{\infty}R_{b}(k))$, or equivalently to select a $w$ that minimizes
where $W_{B}(kw)=\max\{0,1-|kw|\}$ is the Bartlett window, with the bandwidth $w$ being the inverse of $l$. The optimal $w$ that minimizes the problem ((ref)) is
To estimate $w^{*}$, one usually needs to use other types of kernels to estimate $\sum_{k=-\infty}^{\infty}R_{b}(k)$ and $\sum_{k=-\infty}^{\infty}|k|R_{b}(k)$.\\ Our two-step procedure will then be estimating $R_{b}(k)$ by $\hat{R}_{b}(k)$ and applying B\"uhlmann's iterative scheme to find the optimal subsample size $l_{opt}$. First, to estimate $R_{b}(k)=\mathbb{E}[b_{t}b_{t+k}]$, by the properties of $\{a_{n}\}$, $\{b_{t}\}$, and $\{e_{nt}\}$, for any $n\neq m$ and $k\in\{0\}\cup\mathbb{N}$, we have
Then, the estimator $\hat{R}_{b}(k)=\frac{1}{N(N-1)T}\sum_{n\neq m}\sum_{t=1}^{T-k}X_{nt}X_{m,t+k}$ consistently estimates $\mathbb{E}[X_{nt}X_{m,t+k}]$. We divide the time sum by $T$ instead of $T-k$ in order to decrease the volatility. When $k$ is close to $T$, there are fewer available data points for $X_{nt}X_{m,t+k}$. Hence, the quantity is very sensitive to the sample realization. By dividing $T$ instead of $T-k$, the influence of $X_{nt}X_{m,t+k}$ is decreased when $k$ is large. The iterative steps are as follows
where
The number $L$ is the number of iterations of the “global steps". According to buhlmann_locally_1996, $L$ needs to be no less than 4 in order for the optimal bandwidth $w_{opt}$ to have the correct asymptotic order. In the later simulation section, we set $L=20$ to guarantee a convergence.\\ The optimality of such algorithm depends on a set of assumptions. First, we need ((ref)) to be a good approximation with estimating $R_{b}(k)$ by $\hat{R}_{b}(k)=\frac{1}{N(N-1)T}\sum_{n\neq m}\sum_{t=1}^{T-k}X_{nt}X_{m,t+k}$, namely,
as $N,T\rightarrow\infty$. The spectral density at frequency 0 also has to be positive, $SP_{b}(0)>0$. Furthermore, we need $\sum_{k=0}^{\infty}(k+1)^{4}|R_{b}(k)|<\infty$, and the cumulants of $\{b_{t}\}$ with order $h\leq8$ are summable,
With these assumptions, $l_{opt}=l^{*}(1+O(T^{-\frac{2}{7}}))$ as $N,T\rightarrow\infty$, where $l^{*}=\arg\min_{l\in\mathbb{N}} AMSE(\hat{\sigma}_{N,T,b,l}^{2},V)$.
In the last section, we have pinned down the optimal subsample sizes. However, such choice does not eliminate the bias in general due to bias-variance trade-off. In this section, we will introduce a bias-corrected subsampling variance estimator that reduces the order of bias under the framework ((ref)) and Assumption (ref)-(ref).\\ Recall that
then the bias is
The detail is in Appendix (ref). The term $(\frac{b}{l}-c)\sum_{k=-\infty}^{\infty}R_{b}(k)$ involves the unknown constant $c=\lim\frac{N}{T}$, then for simplicity of exposition, we assume $\frac{N}{T}$ is close enough to $c$ such that $(\frac{N}{T}-c)\sum_{k=-\infty}^{\infty}R_{b}(k)$ is in a small order. With $l=l_{opt}$ defined in ((ref)) and $b=\frac{N}{T}\cdot l$, $b,l=O(T^{\frac{1}{3}})=O(N^{\frac{1}{3}})$, and ((ref)) simplifies to
Due to the finiteness of $\sum_{k=-\infty}^{\infty}|k|R_{b}(k)$ (see Appendix (ref)), $Bias(\hat{\sigma}_{N,T,b,l}^{2},V)=O(\frac{1}{l})$. This inspires the following bias-corrected variance estimator
for some constant $D$, where $\Tilde{b}\ll b$, $\Tilde{l}\ll l$ and $\Tilde{b},\Tilde{l}\rightarrow\infty$ as $N,T\rightarrow\infty$. The idea is that we use $\hat{\sigma}_{N,T,\Tilde{b},\Tilde{l}}^{2}$, which is the same type as $\hat{\sigma}_{N,T,b,l}^{2}$ but with even fewer subsamples, to estimate $\mathbb{E}[\hat{\sigma}_{N,T,b,l}^{2}]$. Therefore $\hat{\sigma}_{N,T,\Tilde{b},\Tilde{l}}^{2}-\hat{\sigma}_{N,T,b,l}^{2}$ estimates $Bias(\hat{\sigma}_{N,T,b,l}^{2},V)$, as $\hat{\sigma}_{N,T,b,l}^{2}$ is already an estimator to $V$. And by choosing D wisely, we can eliminate the leading terms in ((ref)) so that $Bias(\hat{\sigma}_{N,T,b,l}^{2,BC},V)$ has a smaller order compared to $Bias(\hat{\sigma}_{N,T,b,l}^{2},V)$. If we choose $D=\frac{\Tilde{l}}{l-\Tilde{l}}$, then we have
We consider the following linear regression
Define $a_{n}=\mathbb{E}[X_{nt}U_{nt}|\alpha_{n}]$, $b_{t}=\mathbb{E}[X_{nt}U_{nt}|\gamma_{t}]$, and $e_{nt}=X_{nt}U_{nt}-a_{n}-b_{t}$. Suppose we are interested in the OLS estimator $\hat{\beta}_{OLS}=(\frac{1}{NT }\sum_{n=1}^{N}\sum_{t=1}^{T}X_{nt}X_{nt}')^{-1}(\frac{1}{NT}\sum_{n=1}^{N}\sum_{t=1}^{T}X_{nt}Y_{nt})$ and want to estimate its asymptotic variance. Instead of subsampling the whole $\hat{\beta}_{OLS}$ and computing the variance of its subsample distribution, we can subsample the score $\frac{1}{NT}\sum_{n=1}^{N}\sum_{t=1}^{T}X_{nt}U_{nt}$ and compute its subsample distribution variance $\hat{\Sigma}_{N,T,b,l}$ first, and then estimate the asymptotic variance of $\hat{\beta}_{OLS}$ by $\hat{\phi}_{NT}^{-1}\hat{\Sigma}_{N,T,b,l}\hat{\phi}_{NT}^{-1}$, where $\hat{\phi}_{NT}=\frac{1}{NT }\sum_{n=1}^{N}\sum_{t=1}^{T}X_{nt}X_{nt}'$. By doing so, we transform such problem into a mean-estimation problem and allow ourselves to apply the subsample size selection algorithm and bias-correction technique described previously.\\ The problem of this approach is that $U_{nt}$ is unobservable; hence we might need to replace $X_{nt}U_{nt}$ by $X_{nt}\hat{U}_{nt}$ where $\hat{U}_{nt}=Y_{nt}-X_{nt}'\hat{\beta}_{OLS}$. Let $\check{\Sigma}_{N,T,b,l}$ be the variance of the subsample distribution associate to the practical score $\frac{1}{NT}\sum_{n=1}^{N}\sum_{t=1}^{T}X_{nt}\hat{U}_{nt}$. It turns out that the consistency of $\check{\Sigma}_{N,T,b,l}$ does not require any additional assumption other than the ones securing the asymptotic normality of $\hat{\beta}_{OLS}$ and the consistency of $\hat{\Sigma}_{N,T,b,l}$.
The proof is in Appendix (ref). (i)-(iv) of Assumption (ref) are used to established the asymptotic normality of $\sqrt{N}(\hat{\beta}-\beta)$, which is obtained by chiang_standard_2024. (v)-(vi) of Assumption (ref), together with the asymptotic normality of $\sqrt{N}(\hat{\beta}-\beta)$, establish the consistency of $\hat{\Sigma}_{N,T,b,l}$ and $\check{\Sigma}_{N,T,b,l}$ and the asymptotic normality of the corresponding t-statistics.
In this section, we evaluate the validity of the subsampling methods and compare them with other existing methods by simulations. We first consider a model with Gaussian limit in which both subsampling methods and other two-way clustering with serial correlation methods work. We also consider a degenerate and non-Gaussian limit model, in which only the subsampling quantile method should be able to provide the correct asymptotic coverage.
For $\{(Y_{nt},X_{nt},U_{nt}):1\leq n\leq N,1\leq t\leq T\}$, consider a linear regression
where $\beta_{0}=\beta_{1}=1$ and $Z_{nt}=[1\quad X_{nt}]'$. The regressor $X_{nt}$ and error $U_{nt}$ are generated as
with $(\alpha_{n}^{x},\varepsilon_{nt}^{x},\alpha_{n}^{u},\varepsilon_{nt}^{u})$ being mutually independent standard normal random variables. The time effects $(\gamma_{t}^{x},\gamma_{t}^{u})$ are independently generated via AR(1) process as
where the innovation terms $v_{t}^{x}$ and $v_{t}^{u}$ follow the distribution $N(0,1-\rho^{2})$ and are independent to $\gamma_{t+1}^{x}$ and $\gamma_{t+1}^{u}$ respectively. We choose $\rho\in\{0,0.25,0.5,0.75\}$ to see how our methods perform under different levels of serial correlation. Note that $\rho=0$ represents the case of two-way clustering without serial correlation. And our goal is to test the 95% confidence interval coverage probability of the OLS estimator $\hat{\beta}_{1}$ using the quantiles of the subsampling distribution and the subsampling variance estimator.\\ For the quantile method, we draw $N_{b}\cdot q$ many subsamples and estimate $\beta_{1}$ using each subsample by OLS. Then, we find the 2.5% and 97.5% quantile of the subsample OLS estimators $\{\hat{\beta}_{1,b,l,i,k}:1\leq i\leq N_{b},1\leq k\leq q\}$ and denote them as $c_{b,l}^{L}(0.025)$ and $c_{b,l}^{L}(0.975)$ respectively. The 95% confidence interval is constructed as $[\hat{\beta}_{1}-\frac{c_{b,l}^{L}(0.975)}{\sqrt{N}},\quad\hat{\beta}_{1}-\frac{c_{b,l}^{L}(0.025)}{\sqrt{N}}]$, where $\hat{\beta}_{1}$ is the OLS estimator of $\beta_{1}$ using the full sample. Table (ref) shows the coverage probabilities for $\beta_{1}$ with a nominal probability of $95\%$ and 1000 Monte Carlo repetitions under different dependence levels. We vary $(N,T)$ from 50 to 400 for different choices of $\rho$. The coverage probability goes to the nominal level for $\rho\in\{0,0.25,0.5\}$, while there is an under-coverage for $\rho=0.75$. Such phenomenon can also be found in some HAC robust tests, for example, andrews_heteroskedasticity_1991 and andrews_improved_1992. kiefer_new_2005 has pointed out that HAC robust test tends to over-reject in finite samples if the traditional method of selecting bandwidth (where the bandwidth goes to infinity slower than the sample size) is used, especially when the serial correlation is strong. As the serial correlation gets larger, it is expected to pick a larger $l$ that captures more of the serial correlation. This suggests a data-driven method of selecting subsample size is needed for the quantile method as well.\\
For the subsampling variance estimator, instead of subsampling the whole OLS estimators and then computing the subsample variance, we use the two estimators $\hat{V}_{N,T,b,l}^{2}$ and $\check{V}_{N,T,b,l}^{2}$, which are defined in Theorem (ref), as well as their bias-corrected versions. The confidence intervals are constructed using the critical value of a standard normal distribution. The data-driven method discussed in Section (ref) is used to select the subsample size. The block length $l$ is set to be $\max\{4,l_{opt}\}$, where $l_{opt}$ is defined in (ref)), and $b$ is set to be equal to $\frac{N}{T}\cdot l$. The reason that we bound the subsample size from below is to leave room for the smaller subsample used in the bias-correction. $\Tilde{b}$ and $\Tilde{l}$ used for bias-correction are set to be $\lfloor\sqrt{b}\rfloor$ and $\lfloor\sqrt{l}\rfloor$ respectively. Table (ref) reports the coverage probabilities for $\beta_{1}$ of our subsampling variance estimators together with the CHS estimator (chiang_standard_2024) and the Chen-Vogelsang estimator (CV; chen_fixed-b_2024). The nominal probability is $95\%$, and we run 1000 Monte Carlo repetitions for each $\rho\in\{0,0.25,0.5,0.75\}$. The sample size $(N,T)$ vary from 50 to 200 for each choice of $\rho$. \\ The results are as follow. First, both bias-corrected variance estimators, $\hat{V}_{N,T,b,l}^{2,BC}$ and $\check{V}_{N,T,b,l}^{2,BC}$, provide the correct nominal coverage probabilities with large sample size except for $\check{V}_{N,T,b,l}^{2,BC}$ under $\rho=0.75$. When the sample size is small, they have lower coverage due to the usage of a even smaller subsample. $\hat{V}_{N,T,b,l}^{2}$ and $\check{V}_{N,T,b,l}^{2}$ tend to be conservative when the serial correlation is weak. This is largely due to the subsample size selection. When the level of dependence is low, the algorithm usually returns a small optimal block length $l_{opt}$. However, by manually bounding $l$ from below, the effective block length $l$ tends to be longer and leads to a more conservative variance estimation. Meanwhile, CHS and CV are less robust to serial correlation than the subsampling variance estimators. These results are consistent with the asymptotic theory and show the effectiveness of the bias-correction under a moderate sample size.\\
Furthermore, we compare our method with the CHS variance estimator and the wild bootstrap method designed by hounyo_reliable_2024. We compare the 95% confidence interval coverage probabilities for $\beta_{1}$ using $\hat{V}_{N,T,b,l}^{2}$, $\check{V}_{N,T,b,l}^{2}$, their bias-corection counterparts, two kinds of CHS variance estimator, and two kinds of wild bootstrap methods under 1000 Monte Carlo simulations. The sample size $N$ and $T$ are both set to be 100, and $\rho$ takes value between 0 to 0.9. The results are shown in Figure (ref). The blue lines represents $\hat{V}_{N,T,b,l}^{2}$ and $\hat{V}_{N,T,b,l}^{2,BC}$, while $\check{V}_{N,T,b,l}^{2}$ and $\check{V}_{N,T,b,l}^{2,BC}$ are colored in aquamarine. The CHS variance estimator is proposed by chiang_standard_2024, and CHS_V is a variant of CHS variance estimator using a different weight function proposed by hounyo_reliable_2024. MWCB_equal and MWCB_sym are two wild bootstrap methods using different P values. MWCB_equal uses equal-tail bootstrap P values, while MWCB_sym uses symmetric bootstrap P values. As we can see, the wild bootstrap methods remain robust and delivers the correct coverage level across different values of $\rho$. The two CHS variance estimators are relatively sensitive to the serial correlation, as they only have the correct coverage probabilities when $\rho<0.4$. Meanwhile, our theoretical subsampling variance estimators deliver the nominal coverage level and are comparable to the wild bootstrap method except when $\rho=0.9$. The practical subsampling variance estimators behave slightly worse but are still more robust to the serial correlation compared to two CHS estimators. Under the case when $\rho=0.9$ and the time effects are highly correlated, the subsampling variance estimators exhibit significant under-coverage. However, $\hat{V}_{N,T,b,l}^{2,BC}$ still has a coverage probability around 90%, while its non-corrected counterpart only has about 85%. This shows the effectiveness of the bias-correction. Overall, the subsampling method is comparable to the leading one of the existing methods in terms of the inference quality, except when the level of serial correlation is extremely high.\\
In the part, we consider a non-separable heterogeneity model (Example 1.7 in menzel_bootstrap_2021).
$\alpha_{n},\gamma_{t}$, and $\varepsilon_{nt}$ are mutually independent. $\alpha_{n}$ and $\varepsilon_{nt}$ are two standard normal random variables, where $\gamma_{t}$ is an AR(1) process generated as
where the innovation terms $v_{t}$ follows the distribution $N(0,1-\rho^{2})$ and is independent to $\gamma_{t+1}$. Since both $\alpha_{n}$ and $\gamma_{t}$ have an expectation of 0, then $\frac{1}{\sqrt{NT}}\sum_{n=1}^{N}\sum_{t=1}^{T}X_{nt}\rightarrow_{d}15\cdot Z_{1}Z_{2}+0.1\cdot Z_{3}$ as described in Section (ref). The coefficients 15 and 0.1 are chosen to amplify the non-Gaussianity. Our goal is to test the 95% confidence interval coverage probability of the sample average $\frac{1}{\sqrt{NT}}\sum_{n=1}^{N}\sum_{t=1}^{T}X_{nt}$. For the subsampling quantile method, we again use a two-sided equal-tailed confidence interval. We vary $(N,T)$ from 50 to 200 for $\rho\in\{0,0.25,0.5,0.75\}$, and the subsample size $b$ and $l$ are chosen accordingly under each pair of $(N,T)$. As shown in Table (ref), with the same subsample size choice across different $\rho$'s, the subsampling quantile method is slightly conservative when the serial correlation is weaker but remains overall valid. This little drawback can be mitigated by a data-driven bandwidth selection algorithm. However, as far as we know, there is no theoretical justified algorithm designed for such method.
Furthermore, we also compare the subsampling quantile method with the CHS variance estimator and the wild bootstrap method of hounyo_reliable_2024. Under this degenerate non-Gaussian limit setting, CHS variance estimator and the wild bootstrap method are known to be theoretically invalid. The sample size $N$ and $T$ are both set to be 100, and $\rho$ takes value between 0 to 0.9. The results are shown in Figure (ref). With the non-Gaussian limit, the CLT based CHS variance estimator and the t-statistics based wild bootstrap method are too conservative. Especially, the wild bootstrap method consistently has a coverage probability around 99%. Meanwhile the subsampling quantile method can generally provide the correct coverage probabilities except when the serial correlation is extremely strong. When $\rho<0.8$, the subsampling quantile method offers a coverage between 95% and 96.5%. This contrast highlights the advantage of our subsampling quantile method. It works under both Gaussian and non-Gaussian limits, as well as both non-degenerate and degenerate cases if we know whether it is degenerate or not.
Additional simulation results are in Appendix (ref).
In this paper, under a two-way clustering with serial correlation structure, we show the subsampling distribution converges to the limiting distribution in probability pointwise as long as the limit is well-defined. The subsampling distribution is defined to be the empirical distribution using all possible and proper subsets of the original sample. The pointwise convergence can be strengthened to uniform convergence if the limiting distribution is continuous. Hence, the quantiles of the subsampling distribution converge to their counterparts in the limiting distribution, suggesting the confidence intervals constructed from the subsampling distribution have the correct asymptotic coverage probabilities. Most importantly, such convergence does not depend on any specific form of the limiting distribution nor any specific convergence rate. We conduct simulation studies under different limiting distributions and convergence rates to confirm the validity of this quantile approach. Especially, the result suggests that the subsampling quantile method still provides the correct asymptotic coverage probability under a degenerate non-Gaussian limit, where all existing methods on two-way clustering with serial correlation fail. We further propose a consistent variance estimator of the limiting distribution, which is the variance of the subsampling distribution. Besides, a data-driven method to select subsample size for this variance estimator is introduced by connecting this subsampling variance estimator to the well-studied spectral density estimator, when the parameter of interest is the population average. We adopt the iterative method introduced by buhlmann_locally_1996, and the optimal subsample sizes $b$ and $l$ are in the order of $N^{\frac{1}{3}}$ and $T^{\frac{1}{3}}$ respectively. Given this rate, we design a bias-corrected variance estimator which reduces the bias from $O(T^{-\frac{1}{3}})$ to $o(T^{-\frac{1}{3}})$. The bias-reduction is done by using another subsampling variance estimator with smaller subsample sizes to estimate the mean of the variance estimator. The simulation studies suggest our variance estimators are comparable, in terms of the inference quality, to the existing leading method, except when the serial correlation is extremely strong. When the serial correlation is extremely strong, our subsampling variance estimators rely on a large subsample size to deliver the desired coverage level. Nevertheless, when the serial correlation is moderate, our methods can deliver the correct coverage level under a reasonable sample size and are easy to implement.\\ There are at least two improvements in the paper that could be made in the future. First, we have not found a proper subsample size selection algorithm for the quantile method nor the variance estimator with a general parameter of interest. Yet, one possible approach is to use the minimum volatility method. When the subsample size is too large, there are not enough variations between different $\hat{\theta}_{b,l,i,k}$ and they are all close to $\hat{\theta}_{NT}$. Then, the subsampling distribution will be too concentrated, and the subsampling confidence intervals tend to under-cover. When the subsample size is too small, the intervals can either over-cover or under-cover. The under-coverage comes from the possibility that $\frac{\tau_{bl}}{\tau_{NT}}$ being too small and resulting in a very tight confidence interval. When the subsample size is in the right range, we should expect the subsampling confidence intervals having coverage probabilities close to the nominal level. This says that for a certain range of the subsample size, the endpoints of confidence intervals should be smooth as a function of subsample size. Therefore, a heuristic method is to pick subsample size from a sequence such that the endpoints are least volatile at the pick. For the variance estimator, we pick the number that yields the least varying variance estimator. Unfortunately, there is no theorem that supports the optimality of such pick in subsampling. bickel_choice_2008 proposes a similar method that works for m-out-of-n bootstrap, which has a lot of similarity to subsampling, to select the bootstrap sample size $m$ and claims optimality. However, their theorem is built on the important fact that the bootstrap distribution converges anyway even if $\frac{m}{n}\rightarrow1$, while the convergence of subsampling distribution heavily relies on the subsample size in a smaller order of the original sample size. Therefore, a theorem of optimality is still needed under the subsampling framework. Secondly, both inference methods rely on a prior knowledge of the convergence rate, and it appears to be a more urging problem to the quantile method. politis_subsampling_1999 has shown that it is possible to apply subsampling methods to a problem with an unknown convergence rate (see Chapter 8). Yet, their method and result are for cross-sectional or time-series data, not panel data. Although a slight modification is probably enough to make their method work in the panel setting, details await to be shown.\\ It might be possible to extend the theorems for subsampling under two-way clustering to a more sophisticated setting. For example, if the individuals are not i.i.d any more but spatially correlated, subsampling still has a chance to work. lahiri_prediction_1999 has shown that their proposed subsampling method provides accurate approximations to the sampling distributions of various functionals of the spatial cumulative distribution function predictor, especially the quantiles. Then, with an appropriate choice of the time indices, subsampling method could potentially work under two-way clustering with serial correlation and spatial dependence.