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.
107,968 characters · 20 sections · 89 citation commands
Fixed-$b$ Asymptotics for Panel Models with Two-Way Clustering
\setcounter{page}{1} \thispagestyle{empty} \pagestyle{plain}
When carrying out inference in a linear panel model, it is well known that failing to adjust the variance estimator of estimated parameters to allow for different dependence structures in the data can cause over-rejection/under-rejection problems under null hypotheses, which in turn can give misleading empirical findings (see Bertrand2004).
To study different dependence structures and robust variance estimators in panel settings, it is now common to use a component structure model $y_{it}=f(\alpha_{i},\gamma_{t},\varepsilon_{it})$ where the observable data, $y_{it}$, is a function of an individual component, $\alpha_{i}$, a time component, $\gamma_{t}$, and an idiosyncratic component, $\varepsilon_{it}$. See, for example, Davezies2021, mackinnon2021wild, Menzel2021, and CHS_Restat. As a concrete example, suppose $y_{it}=\alpha_{i}+\varepsilon_{it}$ for $i=1,...,N$ and $t=1,...,T$, where $\alpha_{i}$ and $\varepsilon_{it}$ are assumed to be i.i.d random variables. The existence of $\alpha_{i}$ generates serial correlation within group $i$, which is also known as the individual clustering effect. This dependence structure is well-captured by the cluster variance estimator proposed by Liang1986 and Arellano1987. One can also use the \textquotedblleft average of HACs" variance estimator that uses cross-section averages of the heterogeneity and autocorrelation (HAC) robust variance estimator proposed by Newey1987. On the other hand, suppose $y_{it}=\gamma_{t}+\varepsilon_{it}$ where $\gamma_{t}$ is assumed to be an i.i.d sequence of random variables. Cross-sectional/spatial dependence is generated in $y_{it}$ by $\gamma_{t}$ is through the time clustering effect. In this case one can use a variance estimator that clusters over time or use the spatial dependence robust variance estimator proposed by Driscoll1998. Furthermore, if both $\alpha_{i}$ and $\gamma_{t}$ are assumed to be present, e.g. $y_{it}=\alpha_{i}+\gamma_{t}+\varepsilon_{it}$, then the dependence of $\{y_{it}\}$ exists in both the time and cross-section dimensions, also as known as two-way clustering effects. Correspondingly, the two-way/multi-way robust variance estimator proposed by ColinCameron2011 is suitable for this case.
In macroeconomics, the time effects, $\gamma_{t}$, can be regarded as common shocks which are usually serially correlated. Allowing persistence in $\gamma_{t}$ up to a known lag structure, Thompson2011 proposed a truncated variance estimator that is robust to dependence in both the cross-section and time dimensions. Because of unsatisfying finite sample performance of this rectangular-truncated estimator, CHS_Restat propose a Bartlett kernel variant (CHS variance estimator, hereafter) and establish validity of tests based on this variance estimator using asymptotics with the cross-section sample size, $N$, and the time sample size, $T$, jointly going to infinity. The asymptotic results of the CHS variance estimator rely on the assumption that the bandwidth, $M$, goes to infinity as $T$ goes to infinity while the bandwidth sample size ratio, $b=\frac{M}{T}$ is of a small order. As pointed out by neave1970improved and Kiefer2005, the value of $b$ in a given application is a non-zero number that matters for the sampling distribution of the variance estimator. Treating $b$ as shrinking to zero in the asymptotics may miss some important features of finite sample behavior of the variance estimator and test statistics. As noted by Andrews1991, Kiefer2005, and many others, HAC robust tests tend to over-reject in finite samples when standard critical values are used. This is especially true when time dependence is persistent and large bandwidths are used. We document similar findings for tests based on CHS variance estimator in our simulations.
To improve the performance of tests based on CHS variance estimator, we derive fixed-$b$ asymptotic results (see Kiefer2005, sun2008optimal, Vogelsang2012, zhang2013fixed, sun2014let, bester2016fixed, and lls2021). Fixed-$b$ asymptotics captures some important effects of the bandwidth and kernel choices on the finite sample behavior of the variance estimator and tests and provides reference distributions that can be used to obtain critical values that depend on the bandwidth (and kernel). Our asymptotic results are obtained for $N$ and $T$ jointly going to infinity and leverage the joint asymptotic framework developed by Phillips1999. The limiting distribution of tests based on the CHS or BCCHS estimator are not asymptotically pivotal, so we propose a plug-in method of simulating fixed-$b$ critical values. One key finding is that the CHS variance has a multiplicative bias given by $1-b+\frac{1}{3} b^{2}\leq1$ resulting in a downward bias that becomes more pronounced as the bandwidth increases. By simply dividing the CHS variance estimator by $1-b+\frac{1}{3}b^{2}$ we obtain a simple bias-corrected variance estimator that improves the performance of tests based on the CHS variance estimator even without using plug-in fixed-$b$ critical values. We label this bias-corrected CHS variance estimator as BCCHS.
As a purely algebraic result, we show that the CHS variance estimator is the sum of the Arellano cluster and Driscoll-Kraay variance estimators minus the \textquotedblleft averages of HAC\textquotedblright\ variance estimator. We show that dropping the \textquotedblleft averages of HAC\textquotedblright \ component in conjunction with bias correcting the Driscoll-Kraay component removes the asymptotic bias in the CHS variance estimator and has the same fixed-$b$ limit as the BCCHS variance estimator. We label the resulting variance estimator of this second bias correction approach as the DKA (Driscoll-Kraay+Arellano) variance estimator. Similar ideas are also used by davezies2018asymptotic and mackinnon2021wild where they argue the removal of the negative and small order component in the variance estimator brings computational advantage in the sense that the variance estimates are ensured to be positive semi-definite. In our simulations we find that negative CHS variance estimates can occur up to 6.4% of the time. An advantage of the OKA variance estimator is guaranteed positive semi-definiteness. The DKA variance estimator also tends to deliver tests with better finite sample coverage probabilities although there are exceptions: when the data is independent and identically distributed (i.i.d.) in both the cross-section and time dimensions, we show the DKA estimator has a different fixed-$b$ limit and results in tests that are conservative including the case where the bandwidth is small\footnote{In the small bandwidth case we find that the limit of the DKA variance estimator is twice as big as the population variance\ for i.i.d.\ data - a finding similar to Theorem 2 in mackinnon2021wild in a multiway clustering setting. We thank a referee for pointing out the similarity between our results for the DKA variance estimator and the results in mackinnon2021wild for multiway cluster variance estimators when the data is i.i.d.}. The fixed-$b$ limit of the CHS variance estimator is also different in the i.i.d.\ case but tests based on it remain robust when the bandwidth is small.
In a finite sample simulation study, we compare sample coverage probabilities of confidence intervals based on CHS, BCCHS, and DKA variance estimators using critical values from both the standard normal distribution and the fixed-$b$ limits. The fixed-$b$ limits of the test statistics constructed by these three variance estimators are not pivotal, so we use a simulation method to obtain the critical values via a plug-in estimator approach to handle asymptotic nuisance parameters. While the plug-in fixed-$b$ critical values can substantially improve coverage rates relative to using standard critical values when using the CHS variance estimator, improvements from simply using the bias corrections are impressive. In the case of data-dependent bandwidths, the plug-in fixed-$b$ critical values provide further improvements in finite sample coverage probabilities when neither $T$ nor $N$ is large. Conversely, when both $N$ and $T$ are very small, bias correction alone can give more accurate finite sample coverage probabilities than bias correction with plug-in fixed-$b$ critical values. Similar results hold for tests based on the DKA variance estimator.
Overall, four different approaches for within- and across-cluster dependent robust tests are proposed: two are simple bias correction by BCCHS and DKA estimators and the other two are bias correction by BCCHS and DKA estimators with plug-in fixed-$b$ critical values. Even though tests based on BCCHS and DKA are asymptotically equivalent under the main assumptions, their finite sample performance is distinguishable. Based on theory and simulation results, we provide comprehensive empirical guidance hinged on the researcher's assessment of the data, model and priority of the test.
The rest of the paper is organized as follows. In Section 2 we sketch the algebra of the CHS estimator and rewrite it as a linear combination of three well-known variance estimators. In Section 3 we derive fixed-$b$ limiting distributions of CHS-based tests for pooled ordinary least squares (POLS) estimators in a simple location panel model and a linear panel regression model. In Section 4 we derive the fixed-$b$ asymptotic bias of the CHS estimator and propose two bias-corrected variance estimators. We also derive fixed-$b$ limits for tests based on the bias-corrected variance estimators. When the data is i.i.d.\ a key assumption for our asymptotic results no longer holds and we show that the asymptotic limits change in this case. Section 5 presents finite sample simulation results that illustrate the relative performance of $t$-tests based on the variance estimators. Some theoretical results for two-way-fixed-effects (TWFE) estimator are also discussed along with the simulation. In Section 6 we illustrate the practical implications of the bias corrections and use of fixed-$b$ critical values in an empirical example. Section 7 concludes the paper with guidance for empirical practice and a discussion on the limitations of proposed approaches.
We first motivate the estimator of the asymptotic variance of the pooled ordinary least squares (POLS) estimator under arbitrary dependence in both the time and cross-section dimensions. Consider the linear panel model
where $y_{it}$ is the dependent variable, $x_{it}$ is a $k\times1$ vector of covariates, $u_{it}$ is the error term, and $\beta$ is the coefficient vector. Let $\widehat{\beta}$ be the POLS estimator of $\beta$. For illustrative purposes the variance of $\widehat{\beta}$ can be approximated as \[ \text{Var}\left( \widehat{\beta}\right) \approx\widehat{Q}^{-1}\Omega _{NT}\widehat{Q}^{-1}, \] where $\widehat{Q}:= \frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}x_{it} x_{it}^{\prime}$ and $\Omega_{NT}:= \text{Var}\left( \frac{1}{NT}\sum_{i=1} ^{N}{\sum_{t=1}^{T}{v_{it}}}\right) $ with $v_{it}:= x_{it}u_{it}$.
Without imposing assumptions on the dependence structure of $v_{it}$, it has been shown, algebraically, that $\Omega_{NT}$ has the following form (see Thompson2011 and CHS_Restat):
Based on this decomposition of $\Omega_{NT}$, Thompson2011 and CHS_Restat each propose a truncation-type variance estimator. In particular, CHS_Restat replaces the Thompson2011 truncation scheme with a Bartlett kernel and establish the consistency result of their variance estimator while allowing two-way clustering effects with serially correlated stationary time effects.
As an asymptotic approximation, appealing to consistency of the estimated variance allows the asymptotic variance to be treated as known when generating asymptotic critical values for inference. While convenient, such a consistency result does not capture the impact of the choice of $M$ and kernel function on the finite sample behavior of the variance estimator and any resulting size distortions of test statistics. To capture some of the finite sample impacts of the choice of $M$ and kernel, we apply the fixed-$b$ approach of Kiefer2005.
Noticeably, the CHS variance estimator can be decomposed into three well-known variance estimators, which will be helpful when we apply the fixed-$b$ approximation. Using straightforward algebra, one can show that the CHS variance estimator defined in equation (2.12) of CHS_Restat can be rewritten as
where, with the Bartlett kernel defined as $k\left( \frac{m}{M}\right) =1-\frac{m}{M}$ and $M$ being the truncation parameter,
Notice that ((ref)) is the \textquotedblleft cluster by individuals\textquotedblright\ estimator proposed by Liang1986 and Arellano1987, ((ref)) is the \textquotedblleft HAC of cross-section averages" estimator proposed by Driscoll1998, and ((ref)) is the \textquotedblleft average of HACs\textquotedblright\ estimator (see Petersen2009 and Vogelsang2012). In other words, ${\widehat{\Omega}}_\text{CHS}$ is a linear combination of three well-known variance estimators that have been proposed to handle particular forms of dependence structure. While there are some existing asymptotic results for the components in ((ref)) that are potentially relevant (e.g. Hansen2007, Vogelsang2012, and CHS_Restat), these results are derived either under one-way dependence or are not sufficiently comprehensive to directly obtain a fixed-$b$ result for ${\widehat{\Omega} }_\text{CHS}$. Some new theoretical results are needed.
To set ideas and intuition we first focus on a simple panel mean model (panel location model) of a $k\times1$ random vector $y_{it}$ and then extend the analysis to the linear regression case. We use a large-$N$ and large-$T$ framework where $N/T\rightarrow c$ for some constant $c$ such that $0<c<\infty$. As a natural way to model panel data with two-way effects, we follow CHS_Restat and assume that $y_{it}$ is generated as follows.
The time component, ${\gamma}_{t}$, is allowed to have serial correlation given that panel data typically has serial correlation beyond that induced by individual effects. As pointed out by CHS_Restat at the beginning of their Section 3, the data-generating process in Assumption 1 is a strict generalization of the representation developed by hoover1979relations, aldous1981representations and kallenberg1989representation (the so-called AHK representation) precisely because ${\gamma}_{t}$ is allowed to have serial correlation. The AHK representation is not sufficient here because it was developed for data drawn from an infinite array of jointly exchangeable random variables in which case ${\gamma}_{t}$ would not have serial correlation.
Using the representation in Assumption 1, CHS_Restat develop the following decomposition of $y_{it}$. Denoting $a_{i}=\text{E}\left( y_{it}-\theta|\alpha_{i}\right) ,\ g_{t}=\text{E}\left( y_{it}-\theta |\gamma_{t}\right) $, and $e_{it}=(y_{it}-\theta)-a_{i}-g_{t}$, one can decompose $y_{it}-\theta$ as \[ y_{it}-\theta=a_{i}+g_{t}+e_{it}=: v_{it}. \] CHS_Restat show that the components are mean zero and that $e_{it}$ is also mean zero conditional on $a_{i}$ and conditional on $g_{t}$. The individual component, $a_{i}$, is i.i.d.\ across $i$, and the time component, $g_{t}$, is stationary. Conditional on $\gamma_{t}$, the $e_{it}$ component is independent across $i$. Finally, the three components are uncorrelated with each other with $a_{i}$ and $g_{t}$ being independent of each other. See Section 3.1 of CHS_Restat for details.
We can estimate $\theta$ using the pooled sample mean estimator given by $\widehat{\theta}={\left( NT\right) }^{-1}\sum_{i=1}^{N}{\sum_{t=1} ^{T}{y_{it}}}$. Rewriting the sample mean using the component structure representation for $y_{it}$ gives
The CHS_Restat variance estimator of $\widehat{\theta}$ is given by ((ref)) with ${\hat{v}} _{it}=y_{it}-\widehat{\theta}$ used in ((ref)) - ((ref)). To obtain fixed-$b$ results for $\widehat{\Omega}_\text{CHS}$ we rewrite the formula for $\widehat{\Omega}_\text{CHS}$ in terms of the following two partial sum processes of ${\hat{v}}_{it}$:
Note that the $a_{i}$ component drops from ((ref)) because $\sum _{i=1}^{N}{\left( a_{i}-\bar{a}\right) }=0$. The Arellano component ((ref)) of $\widehat{\Omega}_\text{CHS}$ is obviously a simple function of ((ref)) with $t=T$. The HAC components ((ref)) and ((ref)) can be written in terms of ((ref)) and ((ref)) using fixed-$b$ algebra (see Vogelsang2012). Therefore, the CHS_Restat variance estimator has the following equivalent formula:
Define three $k\times k$ matrices $\Lambda_{a}$, $\Lambda_{g}$, and $\Lambda_{e}$ such that: \[ \Lambda_{a}\Lambda_{a}^{\prime}=\text{E}(a_{i}a_{i}^{\prime})\text{, \ \ \ }\Lambda_{g}\Lambda_{g}^{\prime}=\sum_{\ell=-\infty}^{\infty} {\text{E}[g_{t}g_{t+\ell}^{\prime}]}\text{, \ \ \ }\Lambda_{e}\Lambda _{e}^{\prime}=\sum_{\ell=-\infty}^{\infty}{\text{E}[e_{it}e_{i,t+\ell} ^{\prime}].} \] The following assumption is used to obtain an asymptotic result for ((ref)) and a fixed-$b$ asymptotic result for $\widehat{\Omega} _\text{CHS}$. Through out the paper, let $\Vert.\Vert$ denote the Euclidean norm for matrices, and let $\lambda_{min}[.]$ denote the smallest eigenvalue of a square matrix.
Assumption 2(i) assumes the mean of $y_{it}$ exists and $y_{it}$ has finite fourth moments. Assumption 2(ii) assumes weak dependence of $\gamma_{t}$ using a mixing condition. Assumption 2 (i) - (ii) follow CHS_Restat. Assumption 2(iii) is a non-degeneracy restriction on the projected individual and time components. Clearly, when data is i.i.d over both the cross-section and time dimensions, this condition does not hold. Because the fixed-$b$ limits of $\widehat{\Omega}_\text{CHS}$ and its associated test statistics turn out to be different in the i.i.d case, we discuss it separately in Section 4. Assumption 2(iii) also rules out the pathological case described in Example 1.7 of Menzel2021: when $y_{it}=\alpha_{i}\gamma_{t}+\varepsilon_{it}$ with $\text{E}(\alpha_{i})=\text{E}(\gamma_{t})=0$, one can easily verify that $a_{i}=g_{t}=0$, in which case the limiting distribution of appropriately scaled $\widehat{\theta}$ is non-Gaussian. Assumption 2(iv) uses the fixed-$b$ asymptotic nesting for the bandwidth. The following theorem gives an asymptotic result for appropriately scaled $\widehat{\theta}$ and a fixed-$b$ asymptotic result for appropriately scaled $\widehat{\Omega}_\text{CHS}$.
The proof of Theorem 1 is given in the Appendix. The limit of $\sqrt{N}\left( \widehat{\theta}-\theta\right) $ was obtained by CHS_Restat. Because $z_{k}$ and $W_{k}(1)$ are vectors of independent standard normals that are independent of each other, $\Lambda_{a}z_{k}+\sqrt{c}\Lambda_{g}W_{k}(1)$ is a vector of normal random variables with variance-covariance matrix $\Lambda _{a}\Lambda_{a}^{\prime}+c\Lambda_{g}\Lambda_{g}^{\prime}$. The $\Lambda _{g}P\left( b,{\widetilde{W}}\left( r\right) \right) \Lambda_{g}^{\prime}$ component of ((ref)) is equivalent to the fixed-$b$ limit obtained by Kiefer2005 in stationary time series settings. Obviously, ((ref)) is different than the limit obtained by Kiefer2005 because of the $h(b)\Lambda_{a}\Lambda_{a}^{\prime}$ term. As the proof illustrates, this term is the limit of the \textquotedblleft cluster by individuals\textquotedblright\ ((ref)) and \textquotedblleft average of HACs\textquotedblright\ ((ref)) components whereas the $c\Lambda_{g}P\left( b,{\widetilde{W}}_{k}\left( r\right) \right) \Lambda_{g}^{\prime}$ term is the limit of the \textquotedblleft HAC of averages\textquotedblright\ ((ref)). Interestingly, Kiefer2005 showed that
where $I_{k}$ is a $k\times k$ identity matrix. The fact that both terms in the limit of $N\widehat{\Omega}_\text{CHS}$ are proportional to $h(b)$ suggests a simple bias correction that is discussed in Section 4.1. Because of the component structure of ((ref)), the fixed-$b$ limits of $t$ and Wald statistics based on $\widehat{\Omega}_\text{CHS}$ are not pivotal. We provide details on test statistics after extending our results to the case of a linear panel regression.
It is straightforward to extend our results to the case of a linear panel regression given by ((ref)). The POLS estimator of $\beta$ is
where $\widehat{Q}:=\frac{1}{NT}\sum_{i=1}^{N}{\sum_{t=1}^{T}{x_{it}x_{it}^{\prime}}}$ as in Section 2. Following CHS_Restat, we assume the components of the panel regression are generated from the component structure: \[ (y_{it},x_{it}^{\prime},u_{it})^{\prime}=f\left( {\alpha}_{i},{\gamma} _{t},{\varepsilon}_{it}\right) \] where $f$ is an unknown Borel-measurable function, the sequences $\left\{ {\alpha}_{i}\right\} $, $\left\{ {\gamma}_{t}\right\} $, and $\{{\varepsilon}_{it}\}$ are mutually independent, ${\alpha}_{i}$ is i.i.d across $i$, ${\varepsilon}_{it}$ is i.i.d across $i$ and $t$, and ${\gamma }_{t}$ is a strictly stationary serially correlated process. Define the vector $v_{it}=x_{it}u_{it}$. Similar to the simple mean model we can write $a_{i}=\text{E}\left( v_{it}|a_{i}\right) $, $g_{t}=\text{E}\left( v_{it}|\gamma_{t}\right) $, $e_{it}=v_{it}-a_{i}-g_{t}$, giving the decomposition \[ v_{it}=a_{i}+g_{t}+e_{it}. \] The CHS_Restat variance estimator of $\widehat{\beta}$ is given by ((ref)) with ${\hat{v}}_{it}$ in ((ref)) - ((ref)) now defined as ${\hat{v}}_{it}=x_{it}{\hat{u}}_{it}$ where ${\hat{u}}_{it}=y_{it}-x_{it}^{\prime}\widehat{\beta}$ are the POLS residuals.
The following assumption is used to obtain an asymptotic result for ((ref)) and a fixed-$b$ asymptotic result for $\widehat{\Omega}_\text{CHS}$ in the linear panel case.
Assumption 3 can be regarded as a counterpart of Assumptions 1 and 2 with Assumption 3(ii) being strengthened. It is very similar to its counterpart in CHS_Restat with a main difference the use of the fixed-$b$ asymptotic nesting for the bandwidth, $M$. For the same reason mentioned in the previous section, we discuss the case where $(x_{it},u_{it})$ are i.i.d separately in Section 4.
The next theorem presents the joint limit of the POLS estimator and the fixed-$b$ joint limit of CHS variance estimator.
The proof of Theorem 2 is given in the Appendix. We can see that the limiting random variable, $V_{k}(b,c)$, depends on the choice of truncation parameter, $M$, through $b$. The use of the Bartlett kernel is reflected in the functional form of $P\left( b,{\widetilde{W}}_{k}(r)\right) $ as well as the scaling term $h(b)$ on $\Lambda_{a}\Lambda_{a}^{\prime}$. Use of a different kernel would result in different functional forms for these limits. Because of ((ref)), it follows that
The scalar $h(b)$ can be viewed as a multiplicative bias term that depends on the bandwidth sample size ratio, $b=M/T$. We leverage this fact to implement a simple feasible bias correction for the CHS variance estimator that is explored below.
Using the theoretical results developed in this section, we next examine the properties of test statistics based on the POLS estimator and CHS variance estimator. We also analyze tests based on two variants of the CHS variance estimator. One is a bias-corrected estimator. The other is a variance estimator guaranteed to be positive semi-definite that is also bias-corrected.
In regression model ((ref)) we focus on tests of linear hypothesis of the form: \[ \text{H}_{0}:R\beta=r,\text{ \ \ \ }\text{H}_{1}:R\beta\neq r, \] where $R$ is a $q\times k$ matrix ($q\leq k$) with full rank equal to $q$, and $r$ is a $q\times1$ vector. Using $\widehat{\text{V}}_\text{CHS}(\widehat{\beta })$ as given by ((ref)), define a Wald statistic as
When $q=1$, we can define a $t$-statistic as \[ t_\text{CHS}=\frac{{R\widehat{\beta}-r}}{\sqrt{{R}\widehat{\text{V}} _\text{CHS}(\widehat{\beta})R^{\prime}}}. \] Appropriately scaling the numerators and denominators of the test statistics and applying Theorem 2, we obtain under $\text{H}_{0}$:
The limits of $W_\text{CHS}$ and $t_\text{CHS}$ are similar to the fixed-$b$ limits obtained by Kiefer2005 but have distinct differences. First, the form of $V_{k}(b,c)$ depends on two variance matrices rather than one. Second, the variance matrices do not scale out of the statistics. Therefore, the fixed-$b$ limits given by ((ref)) and ((ref)) are not pivotal. We propose a plug-in method for the simulation of critical values from these asymptotic random variables.
For the case where $b$ is small, the fixed-$b$ critical values are close to $\chi_{q}^{2}$ and $N(0,1)$ critical values respectively. This can be seen by computing the probability limits of the asymptotic distributions as $b\rightarrow0$. In particular, using the fact that $\text{plim} _{b\rightarrow0}P\left( b,{\widetilde{W}}_{k}(r)\right) =I_{k}$ (see Kiefer2005), it follows that
where $h(\cdot)$ and $P(\cdot)$ are defined in Theorem 1. Therefore, it follows that
and \[ p\lim_{b\rightarrow0}\left[ \frac{RQ^{-1}B_{k}(c)}{\sqrt{{RQ^{-1}} V_{k}(b,c){Q^{-1}R^{\prime}}}}\right] =\frac{RQ^{-1}B_{k}(c)}{\sqrt {{RQ^{-1}\text{Var}(}B_{k}(c){)Q^{-1}R^{\prime}}}}\sim N(0,1). \] In practice, there will not be a substantial difference between using $\chi_{q}^{2}$ and $N(0,1)$ critical values and fixed-$b$ critical values for small bandwidths. However, for larger bandwidths more reliable inference can be obtained with fixed-$b$ critical values.
We now leverage the form of the mean of the fixed-$b$ limit of the CHS variance estimator as given by ((ref)) to propose a biased corrected version of the CHS variance estimator. The idea is simple. We can scale out the $h(b)$ multiplicative term evaluated at $b=M/T$ to make the CHS variance estimator an asymptotically unbiased estimator of $\Lambda_{a}\Lambda _{a}^{\prime}+c\Lambda_{g}\Lambda_{g}^{\prime}$, the variance of $B_{k}(c)=\Lambda_{a}z_{k}+\sqrt{c}\Lambda_{g}W_{k}(1)$. Define the bias-corrected CHS variance estimators as
and the corresponding Wald and $t$-statistics under the null hypothesis $R\beta=r$ are defined as
Because ${\widehat{\Omega}}_\text{BCCHS}$ is a simple scalar multiple of ${\widehat{\Omega}}_\text{CHS}$, we easily obtain the fixed-$b$ limits
Notice that while the fixed-$b$ limits are different when using the bias-corrected CHS variance estimator, they are scalar multiples of the fixed-$b$ limits when using the original CHS variance estimator. Therefore, the fixed-$b$ critical values of the $W_\text{BCCHS}$ and $t_\text{BCCHS}$ are proportional to the fixed-$b$ critical values of $W_\text{CHS}$ and $t_\text{CHS}$. As long as fixed-$b$ critical values are used, there is no practical effect on inference from using the bias-corrected CHS variance estimator. Where the bias correction matters is when $\chi_{q}^{2}$ and $N(0,1)$ critical values are used. In this case, the bias-corrected CHS variance can provide more accurate finite sample inference. This will be illustrated by our finite sample simulations.
As noted by CHS_Restat, the CHS variance estimator does not ensure positive-definiteness, which is also the case for the clustered estimator proposed by ColinCameron2011. davezies2018asymptotic and mackinnon2021wild point out that the double-counting adjustment term in the estimator of ColinCameron2011 is of small order, and removing the adjustment term has the computational advantage of guaranteeing positive semi-definiteness. Analogously, we can think of ${\widehat{\Omega}}_\text{NW}$, as given by ((ref)), as a double-counting adjustment term. If we exclude this term, the variance estimator becomes the sum of two positive semi-definite terms and is guaranteed to be positive definite. Another motivation for dropping ((ref)) is that, under fixed-$b$ asymptotics, ((ref)) simply contributes downward bias in the estimation of the $\Lambda_{a}\Lambda _{a}^{\prime}$ term of $\text{Var}(B_{k}(c))$ through the $-b+\frac{1}{3} b^{2}$ part of $h(b)$ in the $h(b)\Lambda_{a}\Lambda_{a}^{\prime}$ portion of $V_{k}(b,c)$. Intuitively, the Arellano cluster estimator takes care of the serial correlation introduced by $a_{i}$, and the DK estimator takes care of the cross-section and time dependence introduced by $g_{t}$. From this perspective, ${\widehat{\Omega}}_\text{NW}$ is not needed.
Accordingly, we propose a variance estimator which is the sum of the Arellano variance estimator and the bias-corrected DK variance estimator (labeled as DKA hereafter) defined as \[ {\widehat{\Omega}}_\text{DKA}=:\widehat{\Omega}_\text{A}+h(b)^{-1}\widehat{\Omega} _\text{DK}, \] where $\widehat{\Omega}_\text{A}$ and $\widehat{\Omega}_\text{DK}$ are defined in ((ref)) and ((ref)). Notice that we bias correct the DK component so that the resulting variance estimator is asymptotically unbiased under fixed-$b$ asymptotics. This can improve inference should $\chi_{q}^{2}$ or $N(0,1)$ critical values be used in practice. The following theorem gives the fixed-$b$ limit of the scaled DKA variance estimator.
The proof of Theorem 3 can be found in the Appendix. Define the statistics $W_\text{DKA}$ and $t_\text{DKA}$ analogous to $W_\text{BCCHS}$ and $t_\text{BCCHS}$ using the variance estimator for $\widehat{\beta}$ given by $\widehat{\text{V}}_\text{DKA}\left( \widehat{\beta}\right) ={\widehat{Q}} ^{-1}{\widehat{\Omega}}_\text{DKA}{\widehat{Q}}^{-1}$. Applying Theorems 2 and 3, we obtain the fixed-$b$ limits of Wald/t test statistics associated with DKA variance estimator under the null: \[ W_\text{DKA}\Rightarrow W_\text{BCCHS}^{\infty},\text{ \ \ }t_\text{DKA}\Rightarrow t_\text{BCCHS}^{\infty}, \] which are the same as the limits of $W_\text{BCCHS}$ and $t_\text{BCCHS}$ given by ((ref)) and ((ref)).
While the DKA variance estimator is guaranteed to be positive semi-definite, this useful property comes with a potential cost. As is shown in Theorem 2 of mackinnon2021wild, if the score $x_{it}u_{it}$ is i.i.d. over $i$ and $t$, or if clusters are formed at the intersection between individuals and time\footnote{In our setting where the clustering only happens at individual and time levels, clustering at the intersection is the same as the independence across individuals and times.}, the probability limit of two-way cluster-robust variance estimators that drop the double-counting adjustment term, referred to as a two-term variance estimator, is twice the size of the true variance. In other words, if the researcher believes there is clustering when there is none, the use of a two-term estimator would overestimate the asymptotic variance. The associated Wald and $t$-statistics will be scaled down causing over-coverage (under-rejection) problems under the null hypothesis. The following assumption and theorem give fixed-$b$ results for the CHS, BCCHS and DKA statistics for the case of i.i.d. data.
Theorem 4 shows that the fixed-$b$ limits in the i.i.d. case are different for all three test statistics than the limits given by ((ref)), ((ref)) for CHS and ((ref)), ((ref)) for BCCHS and DKA.
Suppose tests are carried out using $\chi_{q}^{2}$ \ and $N(0,1)$ critical values. The limits in Theorem 4 can be used to compute asymptotic null rejection probabilities, or equivalently, asymptotic coverage probabilities for the case of i.i.d. data. For a two-tailed 5% $t$-test, the coverage probabilities are given by \[ P\left( \left\vert t_\text{CHS}^{\infty,iid}\right\vert \leq1.96\right) ,\text{ \ }P\left( \left\vert t_\text{BCCHS}^{\infty,iid}\right\vert \leq1.96\right) ,\text{ \ }P\left( \left\vert t_\text{DKA}^{\infty,iid}\right\vert \leq 1.96\right) . \] For small bandwidths, we can analytically compute these asymptotic coverage probabilities. As $b\rightarrow0$, the limits of $t_\text{CHS}^{\infty,iid}$ and $t_\text{BCCHS}^{\infty,iid}$ converge to $N(0,1)$ random variables giving asymptotic coverage of 95%. In contrast as $b\rightarrow0$, the limit of $t_\text{DKA}^{\infty,iid}$ is a $N(0,\frac{1}{2})$ random variable and the asymptotic coverage is 99.4%, and DKA over-covers and is conservative. This result for DKA tests is similar to Corollary 1 of mackinnon2021wild. For non-small bandwidths the limiting random variables are non-standard. We used simulation methods to compute these probabilities. We approximated the Wiener processes using scaled partial sums of 1,000 i.i.d. $N(0,1)$ random increments and used 50,000 replications to simulate the percentiles.
Table 4.1 reports 97.5% critical values for $t_\text{CHS}^{\infty,iid}$, $t_\text{BCCHS}^{\infty,iid}$, and $t_\text{DKA}^{\infty,iid}$ for a range of values of $b$ that will be used in our finite sample simulations. The critical values of $t_\text{CHS}^{\infty,iid}$ and $t_\text{BCCHS}^{\infty,iid}$ equal 1.96 when $b=0$ and increase as $b$ increase. This suggests that CHS and BCCHS tests will under-cover when the data is i.i.d. and bandwidths are not small. In contrast, the critical values of $t_\text{DKA}^{\infty,iid}$ are always smaller than 1.96 and remain smaller as $b$ increases. Thus, DKA tests over-cover regardless of the bandwidth.
Table 4.1\ also reports asymptotic coverage probabilities using the $N(0,1)$ critical value. We see that as $b$ goes from $0$ to $1.0$, coverage decreases from 95% to 66.7% for CHS, 95% to 87.2% for BCCHS, and is always close to 99% for DKA. These asymptotic calculations predict that CHS and BCCHS will over-reject (be liberal) when data is i.i.d. and non-small bandwidths are used. DKA is predicted to be conservative regardless of bandwidth. The table also reports some results for a random variable, $\hat{t}_\text{BCCHS} ^{\infty,iid}$, that is discussed in the next section.
As we have noted, the fixed-$b$ limits of the test statistics given by ((ref)), ((ref)) and ((ref)), ((ref)) are not pivotal due to the nuisance parameters $\Lambda_{a}$ and $\Lambda_{g}$. A feasible method for obtaining asymptotic critical values is to use simulation methods with unknown nuisance parameters replaced with estimators, i.e. use a plug-in simulation method.
To estimate $\Lambda_{a}$ and $\Lambda_{g}$ we use the estimators:
where $b_\text{dk}=\frac{M_\text{dk}}{T}$ and $M_\text{dk}$ is the truncation parameter for the Driscoll-Kraay variance estimator.\footnote{Note that, in principle, $b_\text{dk}$ can be different from the $b$ used for CHS variance estimator. For simulating asymptotic critical values we used the data dependent rule of Andrews (1991) to obtain $b_\text{dk}$.} The consistency of $\widehat{\Lambda_{a}\Lambda_{a}^{\prime}}$ is given by ((ref)) in the proof of Theorem 3:
And by ((ref)) in the proof of Theorem 2, we have,
Therefore, $\widehat{{\Lambda}_{a}{\Lambda}_{a}^{\prime}}$ is a consistent estimator for ${\Lambda}_{a}{\Lambda} _{a}^{\prime}$ and $\widehat{{\Lambda}_{g}{\Lambda} _{g}^{\prime}}$ is a bias-corrected estimator of ${\Lambda} _{g}{\Lambda}_{g}^{\prime}$ with the mean of the limit equal to ${\Lambda}_{g}{\Lambda}_{g}^{\prime}$ and the limit converges to ${\Lambda}_{g}{\Lambda}_{g}^{\prime}$ as $b_\text{dk}\rightarrow0$. The matrices ${\widehat{\Lambda}}_{a}$ and ${\widehat{\Lambda}}_{g}$ are matrix square roots of $\widehat {{\Lambda}_{a}{\Lambda}_{a}^{\prime}}$ and $\widehat {{\Lambda}_{g}{\Lambda}_{g}^{\prime}}$ respectively such ${\widehat{\Lambda}}_{a}{\widehat{\Lambda}}_{a}^{\prime }=\widehat{{\Lambda}_{a}{\Lambda}_{a}^{\prime}}$ and ${\widehat{\Lambda}}_{g}{\widehat{\Lambda}}_{g}^{\prime }=\widehat{{\Lambda}_{g}{\Lambda}_{g}^{\prime}}$
We propose the following plug-in method for simulating the asymptotic critical values of the fixed-$b$ limits. Details are given for a $t$-test with the modifications needed for a Wald test being obvious.
It is clear that under Assumption 3, as $N,T\rightarrow\infty$ and then as $b_\text{dk}\rightarrow0$, $\hat{t}_\text{CHS}$ converges weakly to the fixed-$b$ limit of $t_\text{CHS}$ in ((ref)) using results in ((ref)) and ((ref)); and so does $\hat{t}_\text{BCCHS}$ ($\hat{t}_\text{DKA}$). However, it is less clear what $\hat{t}_\text{CHS}$ and $\hat{t}_\text{BCCHS}$ ($\hat{t}_\text{DKA}$) are estimating when data is i.i.d.\ across both individual and time dimensions. In the i.i.d.\ case $\frac{1}{N} \widehat{\Lambda_{a}\Lambda_{a}^{\prime}}$ and $\frac{1}{T}\widehat {{\Lambda}_{g}{\Lambda}_{g}^{\prime}}$ each estimate $\Lambda_{xu}\Lambda_{xu}^{\prime}$ (the variance of $x_{it}u_{ut}$). Treating $\frac{1}{N}\widehat{\Lambda_{a}\Lambda_{a}^{\prime}}$ , $\frac{1}{T} \widehat{{\Lambda}_{g}{\Lambda}_{g}^{\prime}}$, and $\widehat{Q}$ as consistent plug-in estimators, it can be shown using arguments similar to the proof of Theorem 4 that $\hat{t}_\text{CHS}$ and $\hat{t}_\text{BCCHS}$ ($\hat{t}_\text{DKA}$) are simulating from the random variables \[ \hat{t}_\text{CHS}^{\infty,iid}=\frac{z_{1}+W_{1}(1)}{\sqrt{h(b)+P\left( b,\widetilde{W}_{1}(r)\right) }},\text{ \ \ }\hat{t}_\text{BCCHS}^{\infty ,iid}=\hat{t}_\text{DKA}^{\infty,iid}=h(b)^{1/2}\hat{t}_\text{CHS}^{\infty ,iid}. \] While these random variables do not depend on $c$ or nuisance parameters, they are clearly different than the limits given by Theorem 4. If we take the probability limit of these random variables as $b\rightarrow0$, it is easy to see\footnote{Obviously $\left( 1-b+\frac{1}{3}b^{2}\right) \rightarrow1$ as $b\rightarrow0$, and recall that $p\lim_{b\rightarrow0}P\left( b,\widetilde {W}_{q}(r)\right) =I_{q}$.} that both random variables converge to \[ \frac{z_{1}+W_{1}(1)}{\sqrt{2}}\sim N(0,1), \] because $\left( z_{1}+W_{1}(1)\right) \sim N(0,2)$. Recall from Theorem 4 that the fixed-$b$ limits of $t_\text{CHS}$ and $t_\text{BCCHS}$ are also approximately $N(0,1)$ when $b$ is small. Thus, the simulated critical values for $t_\text{CHS}$ and $t_\text{BCCHS}$ adapt to the i.i.d.\ case at least for small bandwidths. In contrast, the limit of $t_\text{DKA}$ in Theorem 4 is approximately $N(0,\frac {1}{2})$ when $b$ is small whereas $\hat{t}_\text{DKA}$ is simulating from a $N(0,1)$ random variable. Therefore, simulated critical values for $t_\text{DKA}$ do not adapt to i.i.d.\ data when $b$ is small and $t_\text{DKA}$ over-covers and is conservative.
When the plug-in critical values are used, we can make theoretical predictions for coverage probabilities in the i.i.d. case for bandwidths that are not small ($b>0$) by computing coverage probabilities of the limiting random variables $t_\text{BCCHS}^{\infty,iid}$ and $t_\text{DKA}^{\infty,iid}$ using critical values from the asymptotic random variable $\hat{t}_\text{BCCHS}^{\infty,iid}$ (same as $\hat{t}_\text{DKA}^{\infty,iid}$)\footnote{Coverage probabilities are the same for $t_\text{CHS}^{\infty,iid}$ and $t_\text{BCCHS}^{\infty,iid}$ using critical values from $\hat{t}_\text{CHS}^{\infty,iid}$ and $\hat{t} _\text{BCCHS}^{\infty,iid}$ given the common scaling factor $h(b)^{1/2}$.}. Results are given in Table 4.1 in the $\hat{t}_\text{BCCHS}^{\infty,iid}$ columns. As shown in the table, the critical values of $\hat{t}_\text{BCCHS}^{\infty,iid}$ increase with $b$ but slowly. This helps reduce the under-rejection problems of BCCHS but does not remove them as we see in the coverage probability column for BCCHS that uses critical values from $\hat{t}_\text{BCCHS}^{\infty,iid}$. When DKA uses critical values from $\hat{t}_\text{BCCHS}^{\infty,iid}$, coverage probabilities are similar to the $N(0,1)$, and the coverages do not vary much across $b$ because the critical values of $t_\text{DKA}^{\infty,iid}$ and $\hat{t}_\text{BCCHS}^{\infty,iid}$ roughly move together as $b$ increases.
The asymptotic calculations in Table 4.1 predict that CHS and BCCHS will tend to under-cover (liberal) when the data is i.i.d. with the coverage approaching the nominal level for small bandwidths. DKA is predicted to have over-coverage (conservative) when the data is i.i.d. regardless of the bandwidth.
To illustrate the finite sample performance of the various variance estimators and corresponding test statistics, we present a Monte Carlo simulation study with 10,000 replications in all cases. We focus on a simple linear panel model:
where the true parameters are $(\beta_{0},\beta_{1})=(1,1)$. To allow direct comparisons with Table 1 of CHS_working, we consider a data generating process (DGP) that is linear in the components:
where the latent components $\{{\alpha}_{i}^{x},{\alpha}_{i}^{u},{\varepsilon }_{it}^{x},{\varepsilon}_{it}^{u}\}$ are each i.i.d $N(0,1)$, and the error terms $\widetilde{\gamma}_{t}^{(j)}$ for the AR(1) processes are i.i.d $N(0,1-\rho_{\gamma}^{2})$ for $j=x,u$. The component weights $(\omega_\alpha,\omega_\gamma,\omega_\varepsilon)$ are used to adjust the relative importance of those components.
To further explore the role played by the component structure representation, we consider a second DGP where the latent components enter $x_{it}$ and $u_{it}$ in a non-linear way:
where $\Phi(\cdot)$ is the cumulative distribution function of a standard normal distribution and the latent components are generated in the same way as DGP(1).
Sample coverage probabilities of 95% confidence intervals for $\widehat {{\beta}}_{1}$, the OLS estimator of the slope parameter from ((ref)), are provided for the following variance estimators: Eicker-Huber-White (EHW), cluster-by-$i$ (Ci), cluster-by-$t$ (Ct), DK, CHS, BCCHS, and DKA. For the variance estimators that require a bandwidth choice (DK, CHS, BCCHS, and DKA) we report results using the Andrews1991 AR(1) plug-in data-dependent bandwidth, labeled as $\hat{M}$, designed to minimize the approximate mean square error of a variance estimator (same formula for all four variance estimators). In the case of a scalar $x_{it}$, the formula is given by\footnote{Using equation (6.4) from Andrews1991, we use 0 weight for constant regressor and the weights equal to the inverse of the squared innovation variances for other regressors. Because CHS_Restat parameterize the Bartlett kernel as $1-\frac{m}{M+1}$ whereas we use $1-\frac{m}{M}$, we add 1 to the data-dependent formula so that our Bartlett weights match those used by CHS_Restat.} \[ \hat{M}=1.8171\left( \frac{\widehat{\rho}^{2}}{\left( 1-\widehat{\rho }^{2}\right) ^{2}}\right) ^{1/3}T^{1/3} + 1, \] where $\widehat{\rho}$ is the OLS estimator from the regression $\bar {\hat{v}}_{t}=\rho\bar{\hat{v}}_{t-1}+\eta_{t}$ where $\bar{\hat{v}}_{t}=\frac{1}{N}\sum_{i=1}^{N}\hat{v}_{it}$, $\hat{v}_{it}=x_{it}\hat{u}_{it}$, and $\hat{u}_{it}$ are the OLS residuals from ((ref)). We label the ratio of $\hat{M}$ relative to the time sample size as $\hat{b}=\hat{M}/T$. In some cases $\hat{M}$ can exceed $T$ especially when the time dependence is strong relative to $T$. Therefore, we truncate $\hat{M}$ at $T$ whenever $\hat{M}>T$. We also report results for a grid of bandwidth choices. For tests based on CHS and DKA, we use both the standard normal critical values and the plug-in fixed-$b$ critical values. The simulated critical values use 1000 replications with the Wiener process approximated by scaled partial sums of 500 independent increments drawn from a standard normal distribution. While these are relatively small numbers of replications and increments for an asymptotic critical value simulation, it was necessitated by computational considerations given the need to run an asymptotic critical value simulation for each replication of the finite sample simulation.
We first focus on DGP(1) to make direct comparisons to the simulation results of Table 1 of CHS_working, a working paper version of CHS_Restat\footnote{The reason we refer to the 2022 working paper version of CHS_Restat is because their results for small sample sizes $(N,T) = (25,25)$ are not included in the published paper, CHS_Restat.}. Empirical null coverage probabilities of the confidence intervals for ${\widehat{\beta}}_{1}$ are presented in Table (ref). We start with both\ the cross-section and time sample sizes equal to 25. The weights on the latent components are ${\omega}_{\alpha} =0.25$, ${\omega}_{\gamma}=0.5$, ${\omega}_{\varepsilon}=0.25$. Because of the relatively large weight on the common time effect, $\gamma_{t}$, the cross-section dependence dominates the time dependence. We can see that the confidence intervals using EHW, Ci, and Ct \footnote{Finite sample adjustments are applied to these three variance estimators. $\text{HC}_{1}$ is used for EHW estimator. The \textquotedblleft cluster-by-$i$\textquotedblright, and \textquotedblleft cluster-by-$t$\textquotedblright\ are also adjusted by the usual degrees-of-freedom factor.} suffer from a severe under-coverage problems as they fail to capture both cross-section and time dependence.
With the time effect, $\gamma_{t}$, being mildly persistent ($\rho_{\gamma }=0.425$), the DK and CHS confidence intervals using the normal approximation undercover with small bandwidths with empirical rejection rates mostly below 0.85. The under-coverage problem becomes more severe as $M$ increases because of the well-known downward bias in kernel variance estimators that reflects the need to estimate ${\beta}_{0}$ and ${\beta}_{1}$. Coverages of DK and CHS using $\hat{M}$ are similar to the smaller bandwidth cases, e.g. $M=2$ or $3$, which makes sense given that the average $\hat{M}$ across replications is 2.6 (about 0.1 in terms of $\hat{b}$). However, as the note to the table indicates, large values of $\hat{M}$ can occur in which case $\hat{b}$ is not close to zero. Because they are bias corrected, the BCCHS and DKA variance estimators provide coverage that is less sensitive to the bandwidth. This is particularly true for DKA. If the plug-in fixed-$b$ critical values are used, coverages are closest to 0.95 and very stable across bandwidths with DKA having the best coverage. Because the CHS variance estimator is not guaranteed to be positive definite, we report the number of times that CHS/BCCHS estimates are negative out of the 10,000 replications. In Table 1 there were no cases where CHS/BCCHS estimates are negative.
Tables 5.2 - 5.5 give results for DGP(2) where the latent components enter in a non-linear way. Tables 5.2 - 5.4 have both sample sizes equal to 25 with weights across latent components being the same as DGP(1) ($\omega_{\alpha }=\omega_{\varepsilon}=0.25$, $\omega_{\gamma}=0.5$). Table 5.2 has mild persistence in $\gamma_{t}$ ($\rho_{\gamma}=0.25$). Table 5.3 has moderate persistence ($\rho_{\gamma}=0.5$) and Table 5.4 has strong persistence ($\rho_{\gamma}=0.75$). Tables 5.2-5.4 have similar patterns as Table (ref): confidence intervals with variance estimators non-robust to individual or time components under-cover with the under-coverage problem increasing with $\rho_{\gamma}$. With $\rho_{\gamma}=0.25$, CHS has reasonable coverage (about 0.86) with small bandwidths but under-covers severely with large bandwidths. BCCHS performs much better because of the bias correction and fixed-$b$ critical values provide some additional modest improvements. DKA has better coverage especially when fixed-$b$ critical values are used with large bandwidths. As $\rho_{\gamma}$ increases, all approaches have increasing under-coverage problems with DKA continuing to perform best. Table 5.5 has the same configuration as Table 5.4 but with both sample sizes increased to 50. Both BCCHS and DKA show some improvements in coverage. This illustrates the well-known trade-off between the sample size and magnitude of persistence for accuracy of asymptotic approximations with dependent data. Regarding bandwidth choice, the data-dependent bandwidth performs reasonably well for CHS, BCCHS, and DKA. Finally, the chances of CHS/BCCHS being negative are very small but not zero and chances decrease as both $N$ and $T$ increase.
To show that large values of $\hat{M}$ are not unusual in DGP(2), we report in Figure 1 the frequency of $\hat{b}$ among the 10,000 Monte Carlo replications used in Table 5.4. In this case, more than 21% of replications have $\hat{b}\geq0.2$. This explains why bias correction and fixed-$b$ critical values noticeably reduce the under-coverage problem when $\widehat {M}$ is used.
To show how the relative values of $N$ and $T$ can matter in practice, we provide additional results for the same DGP as Tables 5.4 and 5.5 for $N$ and $T$ over a range of values\footnote{We thank a referee for this suggestion,}. The results are given in Table 5.6. There are two main takeaways from the table: i) bias correction with and without fixed-$b$ critical values always improves coverage probabilities relative to the original CHS test, ii) bias correction alone does slightly better than bias correction with fixed-$b$ critical values when both $N$, $T$ are extremely small ($N=T=10$).
The middle row gives results for both $N$ and $T$ equal to 10. Here the under-coverage problem is substantial for the original CHS test. Bias correction helps especially if the DKA estimator is used. Interestingly, fixed-$b$ critical values help relative to CHS but less than bias correction alone. This is not surprising because the simulated critical values are functions of variance estimators based on small sample sizes. Going up the rows maintains $N=10$ with $T$ increasing to $160$. As expected, coverages approach 0.95 as $T$ increases. The top four rows have $T=160$ with $N$ increasing from $10$ to $80$. With $T$ fixed and $N$ increasing, the first three tests that fail to capture within-time/cross-sectional dependence have deteriorating coverage (undercoverage). In contrast the DK and CHS tests perform well in those cases; bias correction with fixed-$b$ critical values continues to provide further improvements.
Going down from the middle rows shows what happens as $N$ increases when $T$ is small ($T=10$). Coverage of the original CHS tests remains quite low as $N$ increases. Bias correction without fixed-$b$ critical values improves coverage rates but coverage does not improve as $N$ increases. Bias correction with fixed-$b$ critical values performs best and improves as $N$ increases. The results for DKA are interesting. As $N$ increases, undercoverage becomes more severe when normal critical values are used, whereas with fixed-$b$ critical values coverge is best and stable across $N$. The bottom four rows hold $N$ fixed at $160$ and show what happens as $T$ increases from $10$ to $80$. CHS and the bias-corrected versions show better coverage as $T$ increases. CHS and DKA with fixed-$b$ critical values perform best in these cases.
In Theorems 1 - 3, a non-degeneracy assumption on the components is imposed. A special case that violates this assumption is i.i.d.\ data in both individual and time dimensions (random sampling). As we showed in Theorem 4, the fixed-$b$ limit of the test statistics is different in the i.i.d.\ case. By setting ${\omega}_{\alpha}=0$, ${\omega}_{\gamma}=0$ in DGP(1), we present coverage probabilities for the i.i.d case in Table 5.7. There are some important differences between the coverage probabilities in Table 5.7 relative to previous tables. First, notice that the coverages using EHW, Ci, and Ct are close to the nominal level as one would expect. The patterns of coverage probabilities for CHS, BCCHS and DKA are as predicted by Theorem 4 and the asymptotic calculations given in Table 4.1. Coverages of CHS are close to 0.89 for small bandwidths and under-coverage problems occur with larger bandwidths. BCCHS is less prone to under-coverage as the bandwidth increases and plug-in fixed-$b$ critical values help to reduce, but do not eliminate, the under-coverage problem. In contrast DKA over-covers regardless of the bandwidth and whether or not fixed-$b$ critical values are used. As $N$,$T$ get larger, we would expect the coverages of CHS/BCCHS to approach 95% in the i.i.d. case (assuming a small bandwidth) but not for DKA where over-coverage would persist.
To gauge the extent to which the under-coverage of CHS/BCCHS and over-coverage of DKA is caused by the mis-match between the plug-in fixed-$b$ critical values and the i.i.d. fixed-$b$ limits, we report in Table 5.8 simulated coverage probabilities using fixed-$b$ critical values based on the limits in Theorem 4. We see that, regardless of the bandwidth, coverages are much closer to 95%. Therefore, a significant portion of the size distortions in Table 5.7 is because of the mis-match.
These results raise a practical question: how does a researcher know which case is being dealt with? In panel data settings, random sampling is almost never a reasonable assumption in the time dimension and clustering dependence often exists due to unobserved heterogeneity in both individual and time dimensions. Some concern may arise as it is common in practice for empirical researchers to include fixed-effect dummy variables, also as known as two-way fixed-effect estimator (TWFE, hereafter), to remove at least some of dependence generated by individual and time unobserved heterogeneity, and in some cases, such as DGP (1), all of the dependence structure would be removed and we are back to the i.i.d case. However, it is important to note that DGP (1) is a very special case. In general, fixed-effect approaches do not guarantee the resulting scores to be free from clustering dependence. Indeed, other data generating mechanisms exist where TWFE will not completely remove the dependence caused by individual and time components in the score as shown in CHS_Restat. We discuss an example and its implications in the next subsection.
One should also note that if we compare absolute size-distortions, DKA-based tests are not necessarily more concerning than BCCHS-based tests: while in opposite directions, the magnitudes of size distortions are mostly comparable and DKA-based tests do a better job as $M$ increases. Moreover, given that the DKA-based tests tend to be more conservative, a rejection using a DKA-based test delivers strong evidence against the null hypothesis. For a researcher that wants to avoid spurious null rejections (relative to the desired significance level), then DKA-based tests are preferred. On the other hand, if a rejection is not obtained with BCCHS tests, this is strong evidence that the null cannot be rejected. Suppose a rejection is obtained with BCCHS but not with DKA. In this case a researcher had to balance potential over-rejections from BCCHS with potential lower power of DKA which depend on the extent to which the researcher thinks two-way clustering is present in the model.
A popular alternative to the pooled OLS estimator is the additive TWFE estimator where individual and time period dummies are included in ((ref)). It is well known that individual and time dummies will project out any latent individual or time components that only linearly enter $x_{it}$ and $u_{it}$ individually (as would be the case in DGP(1)) leaving only variation from the idiosyncratic component $e_{it\text{.}}$. In this case, we would expect the sample coverages of CHS and DKA to be similar to the i.i.d case in Table 5.7. However, under the general component structure representation, the TWFE transformation may not fully remove the individual and time components if they enter in a nonlinear manner and we would expect results for CHS and DKA similar to Tables 5.1-5.6.
As an illustration, in Table 5.9 we report results for the TWFE estimator using the same configuration as Table 5.4 for DGP(2). The sample coverage probabilities are different from Table 5.4 but are very similar to the results in Table 5.7 for the i.i.d. case. Therefore, for DGP(2), the TWFE dummy variables remove the bulk of the variation from the individual and time components.
In contrast CHS_Restat provide an example where the TWFE dummy variables do not remove the component structure. Consider a third DGP given by
where the latent components $\{{\alpha}_{1i},{\alpha}_{2i},{\alpha} _{3i},\gamma_{1t},\gamma_{2t},\gamma_{3t},{\varepsilon}_{it}^{x},{\varepsilon }_{it}^{u}\}$ are $N(0,1)$ random variables that are independent across $i$ and $t$ and independent with each other. As CHS_Restat argue, there is no endogeneity between $x_{it}$ and $u_{it}$, and it is not difficult to show that $E(x_{it}|{\alpha}_{i})=E(x_{it}|\gamma_{t})=E(u_{it}|{\alpha} _{i})=E(u_{it}|\gamma_{t})=0$. While $x_{it}$ and $u_{it}$ do not have the component structure, the score, $x_{it}u_{it}$, does because $E(x_{it} u_{it}|{\alpha}_{i})={\alpha}_{2i}{\alpha}_{3i}$ and $E(x_{it}u_{it} |\gamma_{t})=\gamma_{t2}\gamma_{t3}$. Therefore, the TWFE dummy variables will not remove the component structure from $x_{it}u_{it}$.
Table 5.10 gives results for DGP(3) for TWFE with $N=T=25$. We see that tests based on variance estimators that are not robust to two-way cluster dependence have substantial under-coverage problems. The original CHS does a better job but tends to under-cover with large bandwidths. BCCHS works better and plug-in fixed-$b$ critical values provide additional improvements in coverage. DKA works quite well, with small improvements using plug-in fixed-$b$ critical values, and coverage probabilities are close to 95%.
The results for BCCHS and DKA in Table 5.10 suggest that fixed-$b$ limits given by ((ref)) and ((ref)) for POLS can continue to hold for tests based on the TWFE estimator of $\beta$. Let $\Ddot{x}_{it}$ and $\Ddot{u}_{it}$ denote the individual and time dummy demeaned versions of $x_{it}$ and $u_{it}$ respectively. Suppose that $\Ddot{x}_{it}\Ddot{u}_{it}$, has the individual and time component structure. Because CHS_Restat show that \[ \Ddot{x}_{it}\Ddot{u}_{it}=\Tilde{x}_{it}\Tilde{u}_{it}+o_{p}(1), \] where $\Tilde{x}_{it}=x_{it}-\text{E}[x_{it}|\alpha_{i}]-\text{E} [x_{it}|\gamma_{t}]+\text{E}[x_{it}]$ and $\Tilde{u}_{it}$ is similarly defined, equivalent versions Theorem 2, ((ref)), ((ref)) and Theorem 3 are easily established for the TWFE estimator provided the stronger exogeniety assumption, $\text{E}[\Tilde {x}_{it}u_{it}]=0$, holds\footnote{Strict exogeneity over time, $E(u_{it} |x_{i1},x_{i2},...,x_{iT})=0$, is sufficient for $E[\Tilde{x}_{it}u_{it}]=0$ to hold.}.
We illustrate how the choice of variance estimator affects $t$-tests and confidence intervals using an empirical example from Thompson2011. We test the predictive power of market concentration on the profitability of industries where the market concentration is measured by the Herfindahl-Hirschman Index (HHI, hereafter). This example features data where dependence exists in both cross-section and time dimensions with common shocks being correlated across time.
Specifically, consider the following linear regression model of profitability measured by $\text{ROA}_{m,t}$, the ratio of return on total assets for industry $m$ at time $t$: \[ \text{ROA}_{m,t}=\beta_{0}+\beta_{1}\text{ln}(\text{HHI}_{m,t-1})+\beta_{2}\text{PB}_{m,t-1}+\beta _{3}\text{DB}_{m.t-1}+\beta_{4}\bar{\text{ROA}}_{t-1}+u_{m,t} \] where $\text{PB}$ is the price-to-book ratio, $\text{DB}$ is the dividend-to-book ratio, and $\bar{\text{ROA}}$ is the market average $\text{ROA}$ ratio.
The data set used to estimate the model is composed of 234 industries in the US from 1972 to 2021. We obtain the annual level firm data from Compustat and aggregate it to industry level based on Standard Industry Classification (SIC) codes. The details of data construction can be found in Section 6 and Appendix B of Thompson2011.
In Table 6.1, we present the POLS estimates for the five parameters and $t$-statistics (with the null $\text{H}_{0}:\beta_{j}=0$ for each $j=1,2,...,5$) based on the various variance estimators. We use the data dependent bandwidth, $\hat{M}$, in all relevant cases. We can see the $t$-statistics vary non-trivially across different variance estimators. The estimated coefficient of $\text{ln}(\text{HHI}_{m,t-1})$ is significant at a 1% level based on two-sided t-tests using any standard errors among comparison, including the DKA standard error. As is discussed in Section 5.3, a rejection using DKA is strong evidence of market concentration being powerful in predicting the profitability of industries. On the other hand, the estimated coefficient of $\text{DIV}/\text{Book}$ is significant at the 5% significance level in a two-sided test when EHW, cluster-by-industry, cluster-by-time, and DK variances are used, while it is only marginally significant when CHS is used and marginally insignificant when BCCHS and DKA are used.
In Table 6.2 we present 95% confidence intervals. For CHS/BCCHS and DKA we give confidence interval using both normal and plug-in fixed-$b$ critical values. For the bias corrected variance estimators (BCCHS and DKA) the differences in confidence intervals between normal and fixed-$b$ critical values are not large consistent with our simulation results.
In Table 6.3, we include the results for TWFE estimator to see how the inclusion of firm level and time period dummies matter in practice. The presence of the dummies results in the intercept and $\bar{\text{ROA}}_{t-1}$ being dropped from the regression. Overall, test statistics based on CHS, BCCHS, and DKA agree with each other in magnitude and they are much smaller relative to EHW-based test statistics. As we have seen in Table 5.7, when the scores are independent in both the cross-section and time dimensions, test statistics based on those non-twoway robust standard errors tend to be smaller (higher coverage) on average except for DKA-based tests. Because the non-twoway test statistics in Table 6.3 are larger than the CHS/BCCHS statistics suggests TWFE does not fully remove two-way dependence and two-way cluster-robust standard errors are appropriate.
The 95% confidence intervals for TWFE case are presented in Table 6.4. Confidence intervals tend to be wider with fixed-$b$ critical values. This is expected given that fixed-$b$ critical values are larger in magnitude than standard normal critical values.
This paper investigates the fixed-$b$ asymptotic properties of the CHS variance estimator and tests. An important algebraic observation is that the CHS variance estimator can be expressed as a linear combination of the cluster variance estimator, \textquotedblleft HAC of averages" estimator, and \textquotedblleft average of HACs" estimator. Building upon this observation, we derive fixed-$b$ asymptotic results for the CHS variance estimator when both the sample sizes $N$ and $T$ tend to infinity. Our analysis reveals the presence of an asymptotic bias in the CHS variance estimator which depends on the ratio of the bandwidth parameter, $M$, to the time sample size, $T$. This bias is multiplicative and leads to simple feasible bias corrected version of the CHS variance estimator (BCCHS). We propose a second bias corrected variance estimator, DKA, by dropping the \textquotedblleft HAC of averages" that is guaranteed to be positive semi-definite. We show that the fixed-$b$ limiting distribution of tests based on CHS, BCCHS and DKA are not asymptotically pivotal, and we propose a straightforward plug-in method for simulating fixed-$b$ asymptotic critical values. Overall,\ we propose four test statistics that build on the CHS test: BCCHS and DKA tests using chi-square/standard normal critical values, and BCCHS and DKA tests using plug-in fixed-$b$ critical values\footnote{CHS tests that use simulated fixed-$b$ critical values are exactly equivalent to BCCHS tests based on simulated fixed-$b$ critical values because the fixed-$b$ limits explicitly capture the bias in the CHS variance estimator.}. Extensive simulations studies are reported that compare finite sample performance of the proposed approaches with existing approaches in terms of finite sample null coverage probabilities. The simple bias-correction approaches provide non-trivial improvements\ in coverage probabilities and bias-correction with plug-in fixed-$b$ critical values provide additional improvements except in the i.i.d. case and when both $N$ and $T$ are very small.
Our results clearly suggest that the bias corrected variance estimators, BCCHS and DKA, provide more reliable inference in practice with or without plug-in fixed-$b$ critical values. While plug-in fixed-$b$ critical values involve some computation cost in practice, we can generally recommend fixed-$b$ critical values be used in practice given that i) fixed-$b$ critical values improve finite sample coverage probabilities when large bandwidths are used, ii) data dependent bandwidths can be large, and iii) coverage probabilities with or without fixed-$b$ critical values are similar when bandwidths are small. However, there are important exceptions. When both the cross-section and time sample sizes are very small, then BCCHS and DKA based tests using plug-in fixed-$b$ critical values could yield slightly worse empirical null coverages than using chi-square/standard normal critical values because the plug-in estimators are noisy. Therefore, the choice between using fixed-$b$ or chi-square/standard normal critical values for BCCHS and DKA tests depends on the sample sizes in additional to any relevant computational costs.
The choice between tests based on BCCHS and DKA is nuanced. While DKA ensures positive definiteness and usually provides tests with better empirical null coverage probabilities, these benefits do not come without a cost. Although rare in panel settings, if the scores, $x_{it}u_{it}$, are i.i.d. over both individual and time dimensions, the DKA estimator has a different fixed-$b$ limiting distribution and tests based on the DKA estimator can be conservative. In contrast, while the BCCHS estimator also has a different fixed-$b$ limiting distribution in the i.i.d. case, it has correct asymptotic coverage probabilities when the bandwidth is small. However, if the bandwidth is not small, CHS and BCCHS tests under-cover in the i.i.d. case. Therefore, the practical choice between DKA and BCCHS depends on a researcher's assessment of the data, the model, and the priority of inference. If the data is thought to be independent in both dimensions, then one should not consider cluster-robust variances estimators in the first place. If the data is thought to have individual and serially correlated time cluster dependence and the researcher places higher priority on controlling over-rejections while having a conservative test (with the cost of lower power) should there not be cluster dependence, the DKA estimator is preferred. BCCHS would be preferred if the additional under-coverage relative to DKA is viewed as reasonable in order to have higher power should there not be cluster dependence.
It is important to acknowledge some limitations of our analysis and to highlight areas of future research. We found that finite sample coverage probabilities of all confidence intervals exhibit under-coverage problems when the autocorrelation of the time effects becomes strong relative to the time sample size. In such cases, potential improvements resulting from the fixed-$b$ adjustment is limited. Part of this limitation arises because the test statistics are not asymptotically pivotal, necessitating plug-in simulation of critical values. The estimation uncertainty in the plug-in estimators can introduce sampling errors to the simulated critical values that can be acute when persistence is strong. Finding a variance estimator that results in a pivotal fixed-$b$ limit would help address this problem although appears to be challenging.
An empirically relevant question is whether the component structure is a good approximation when the component representation in Assumption 1 is not exact. Ideally, inferential theory should be studied under a DGP where the dependence is generated not only through individual and time components but also through the idiosyncratic component. Obtaining fixed-$b$ results for this generalization appears challenging. Some unreported simulation results point to some theoretical conjectures but a formal analysis is beyond the scope of this paper and is left for future research.
A second empirically relevant case we do not address in this paper is the unbalanced panel data case. There are several challenges in establishing formal fixed-$b$ asymptotic results for unbalanced panels. Unbalanced panels have time sample sizes that are potentially different across individuals and this potentially complicates the choice of bandwidths for the individual-by-individual variance estimators in the average of HACs component of the variance. For the Driscoll-Kraay component, the averaging by time would have potentially different cross-section sample sizes for each period. Theoretically, obtaining fixed-$b$ results for unbalanced panels also depends on how the missing data is modeled. For example, one might conjecture that if missing observations in the panel occur randomly (missing at random), then extending the fixed-$b$ theory would be straightforward. While that is true in pure time series settings (see rho2019heteroskedasticity), the presence of the individual and time random components in the panel setting complicate things due to the fact that the asymptotic behavior of the components in the partial sums is very different from the balanced panel case. Obtaining useful results for the unbalanced panel case is challenging and is a focus of ongoing research.
We thank three anonymous referees, Antonio Galvao, Jeff Wooldridge, Bruce Hansen, Harold Chiang, and participants in AMES 2023 in Beijing and MEG 2022 in East Lansing for helpful comments and suggestions.