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.
75,990 characters · 15 sections · 42 citation commands
Standard Errors for Two-Way Clustering with Serially Correlated Time Effects
A standard panel data set has observations double-indexed over firms\footnote{The index $i$ can refer to any entity, such as firms, individuals, or households. For simplicity we will refer to these entities as “firms”.} $i$ and time $t$. A panel is said to have a two-way dependence structure if there is dependence across individuals at any given time, and across time for any given individual. A common model for two-way dependence is the components structure $U_{it} = f(\alpha_i,\gamma_t,\varepsilon_{it})$, where $\alpha_i$ is a firm effect, $\gamma_t$ is a time effect, and $\varepsilon_{it}$ is an idiosyncratic effect. It is typical to view the time effects $\gamma_t$ as omitted macroeconomic variables, such as the state of the business cycle. Therefore, they are unlikely to be serially independent. Consequently, it is reasonable to treat $\gamma_t$ as a serially correlated time-series process.
Serial correlation in the common time-effects, however, creates an extra layer of serial dependence beyond two-way dependence. It induces dependence among observations which do not share a common firm or time index. This fundamentally complicates the dependence structure, rendering existing theory and methods inappropriate. Most importantly for practice, existing two-way clustered inference methods do not allow time effects with arbitrary serial correlation. Moreover, no formal asymptotic theory has previously been developed under the two-way clustering setting with serially correlated time effects. This paper sheds new light on literature of two-way clustering by formally investigating asymptotic theory for serially correlated time effects in this context and proposing novel and theoretically supported method of inference under such settings.
The most popular inference method for two-way dependent panels is the two-way clustered standard errors of \citet*{CGM2011}, which we shall henceforth refer to as CGM. A related recently method is the bootstraps of menzel2021bootstrap; see also DDG2019 for a generalization to empirical processes. Both of these approaches allow for the components structure $U_{it} = f(\alpha_i,\gamma_t,\varepsilon_{it})$, but only under the additional strong condition that the time effects $\gamma_t$ are serially independent. The CGM standard errors explicitly calculate the variance allowing for traditional two-way dependence, not allowing for dependence induced by serially correlated time effects. Consequently, these methods exclude, by construction, the possibility that the common time component $\gamma_t$ is an unmodelled macroeconomic effect.
An important intuitive extension due to \citet*{Thompson2011} allows $\gamma_t$ to be serially correlated up to a known fixed number of lags and suggests to estimate the asymptotic variance by including unweighted, lagged autocovariance estimates to the CGM estimator. This relaxes the CGM assumptions by allowing serial correlation structures that are of $m$-dependence. In practice, however, it is difficult to implement since the serial dependence structure is not known a priori. Also, even under this $m$-dependence setting, no asymptotic distribution theory was provided. This is particularly troubling since serial correlated time effects induces a complicated dependence structure and thus was unclear whether this inference procedure is theoretically justified and under which conditions it is so. Indeed, our simulations unveil that this unweighted, fixed number of lags approach shows unsatisfactory finite sample performances under various DGPs; see Section (ref). In addition, based on his own simulations, \citet*{Thompson2011} recommends omitting the correction for serial correlation unless the time dimension is large. Thus in practice, the Thompson estimator actually implemented by most (if not all) empirical researchers reduces to the CGM two-way estimator.
Furthermore, an asymptotic distribution theory for regression with two-way clustering with general serial correlated time effects is missing. \citet*{CGM2011} assert an asymptotic theory for estimation, but do not examine the impact of two-way dependence, nor examine standard error estimation. As previously mentioned, \citet*{Thompson2011} does not provide a distribution theory even under the $m$-dependence setting. DDG2019, mackinnon2021wild, and menzel2021bootstrap do provide rigorous theory, yet only for settings without serial dependence.
Clustered inference can alternatively be based on unstructured one-way dependence (over either $i$ or $t$, but not both simultaneously) using the popular clustered variance estimator of \citet*{liang1986longitudinal} and \citet*{arellano1987computing}. These methods, however, cannot account for two-way dependence. An alternative framework is one-way-cluster dependence across $i$ with weak serial dependence across $t$ driscoll_kraay_1998. A yet alternative framework has been provided by vogelsang2012 and hidalgo2021inference, which study panels with cross-sectional and temporal dependence under a different set of conditions. They allow dependence across firms and time, but the dependence between observations within a time period, as well as the dependence of a cross-sectional unit observed over time, both decay as observations get further apart in time and space. Consequently, these alternative frameworks do not allow arbitrary two-way clustering.
Clustered standard errors have become ubiquitous in applied economic research, as evidenced by a perusal of current applied journals, and by the enormous citations to several of the above-mentioned papers. \citet*{petersen2009estimating} provides an excellent review of these popular methods and their use in empirical research through 2009. Our perusal of current applied journals reveals that nearly all applications use either Liang-Zeger-Arelleno one-way clustering or CGM two-way clustering. While \citet*{Thompson2011} is also highly cited, our review indicates that empirical applications do not employ his correction for correlated time effects, but rather use the simpler CGM two-way clustering.
In this article, we modify the CGM and Thompson two-way clustered standard error to accommodates time effects with arbitrary stationary serial dependence. Our approach allows for cluster dependence within individuals $i$, within time periods $t$, and allows the common time component $\gamma_t$ to be serially dependent of arbitrary order. Ours is the first approach which allows this complexity of two-way dependence. This is accomplished by a correction involving kernel smoothing over the autocorrelations, with the number of autocorrelation lags increasing with sample size. To select the lag truncation parameter, we propose a simple rule based on Andrews1991.
We provide an asymptotic theory of inference under weak regularity conditions, including the assumption that the time effects $\gamma_t$ are strictly stationary and mixing. We show that the least squares estimator is asymptotically normal, our proposed variance estimator is consistent, and t-ratios are asymptotically standard normal, permitting conventional inference. The proofs of these results are far from trivial. For example, our consistency proof for our proposed cluster-robust variance estimator is nonstandard. We show that the problem can be re-written into a claim of bounding a fourth-order sum of cross-moments of dependent time series, for which a mixing bound due to yoshihara1976limiting can be applied. Furthermore, the same proof uses novel projection arguments which simplify the derivations.
We explore the performance of our proposed method in a simple simulation experiment which compares the coverage probability of confidence intervals constructed with six different standard error methods. We find that our proposed method has the best performance relative to the competitors in each simulation design considered, and in some settings the difference is substantial.
We also illustrate the relevance of the method with an empirical application to estimation of the slope coefficients in a standard Fama-French three-factor regression using two panels of stock returns. We find that our proposed standard errors are different -- and larger -- than conventional standard errors, for four of six regression estimates examined.
The Stata command is available to install by ssc install xtregtwo.
The rest of this paper is organized as follows. Section (ref) discusses two-way dependence with correlated time effects and provides an informal overview of the method with a practical guide. Section (ref) presents the formal theoretical results. Section (ref) provides simulation evidence on the practical performance of our proposed method. Section (ref) presents an empirical application to a standard Fama-French regression. The appendix collects a mathematical proof of the main result, auxiliary lemmas and their proofs, and additional details omitted from the main text.
Let $\left(Y_{it},X_{it}'\right)$ be a panel of observations over $i=1,...,N$ and $t=1,...,T$, where $Y_{it}$ is real-valued and $X_{it}$ is a $k \times 1$ vector. The model is the linear regression equation
with
The standard estimator for $\beta$ is least squares
The least squares residuals are $\widehat U_{it}=Y_{it}-X_{it}'\widehat \beta $.
We are interested in the variance of $\widehat \beta$. It can be calculated explicitly under the auxiliary assumption that the regressors are fixed and the error is strictly exogenous.\footnote{The strict exogeneity condition is not required for our main theory and is used only for an illustration purpose in the current section for the exact variance calculations.} We only use this assumption to motivate our covariance matrix estimator, however, and will not be needed for our asymptotic distribution theory.
The variance of $\widehat \beta$ can be written as follows. Define the firm sums $R_i=\sum_{t=1}^T X_{it} U_{it}$, the time sums $S_t=\sum_{i=1}^N X_{it} U_{it}$, and the cross-sums $G_m=\sum_{t=1}^{T-m} S_{t} S_{t+m}'$ and $H_m =\sum_{i=1}^N \sum_{t=1}^{T-m} X_{it} U_{it}X_{i,t+m}' U_{i,t+m}$. With a little algebra we obtain the following decomposition.
where
and
The expression (ref)-(ref) decomposes the variance of the least squares estimator into four components: (ref) is the variance of the firm sums; (ref) is the variance of the time sums; (ref) is a correction for double-counting of the common variance in (ref) and (ref); and (ref) is the autocovariances of the time sums, corrected for double-counting.
Estimators of the variance matrix $V_{NT}$ take the general form
where $\widehat \Omega_{NT}$ is some estimator of $\Omega_{NT}$. Different estimators make distinct assumptions on the covariances in (ref)-(ref) which lead to distinct estimators for $\Omega_{NT}$ in (ref). The Liang-Zeger-Arellano one-way cluster estimator assumes that observations are independent across $i$, implying that (ref)+(ref)+(ref) equals zero. The “cluster within $t$” estimator assumes that observations are independent across $t$, implying that (ref)+(ref)+(ref) equals zero. The CGM two-way estimator assumes that observations $it$ and $js$ are independent if $i\ne j$ or $t\ne s$, implying that (ref) equals zero. The respective estimators take the same form as the assumed non-zero expressions in (ref)-(ref). For example, the CGM variance estimator of $\Omega_{NT}$ is
where $\widehat R_i=\sum_{t=1}^T X_{it}\widehat U_{it}$ and $\widehat S_t=\sum_{i=1}^N X_{it}\widehat U_{it}$.
\citet*{Thompson2011} assumes that the across-firm autocovariances are non-zero for small lags $m$, but zero for lags beyond a known constant $M$. This implies that the sum over $m$ in (ref) can be truncated above $m=M$. This motivates his estimator of $\Omega_{NT}$, which is (ref) plus
where $\widehat G_m=\sum_{t=1}^{T-m} \widehat S_{t} \widehat S_{t+m}'$ and $\widehat H_m =\sum_{t=1}^{T-m}\sum_{i=1}^N X_{it}\widehat U_{it}X_{i,t+m}'\widehat U_{i,t+m}$. \citet*{Thompson2011} does not discuss selection of $M$, other than to indicate that it is known a priori. For his simulations and empirical applications he sets $M=2$, which we take to be his default choice.
We illustrate the dependence patterns assumed by the different estimators in Figure (ref). Each panel shows an array with each entry depicting firm/time pairs $(i,t)$, with the star $\star$ marking the reference point $(i,t)=(1,1)$, and dependence structures indicated by the grey shading. Panel (A) illustrates the case of independent observations (which corresponds to the unclustered Eicker-Huber-White estimator) where the observation $(i,t)=(1,1)$ is uncorrelated with all other observations. Panel (B) illustrates the case of independence across firms (which corresponds to the Liang-Zeger-Arellano one-way cluster estimator) where the observation $(1,1)$ is correlated with $(1,t)$ for $t>1$, but is uncorrelated with all other observations. Panel (C) similarly illustrates the case of independence across time (the “cluster within $t$” estimator). Panel (D) illustrates the case where the observation $(1,1)$ is correlated with $(1,t)$ for $t>1$ and with $(i,1)$ for $i>1$ (corresponding to the CGM two-way clustered estimator). Panel (E) illustrates the case where two-way clustering is augmented to allow dependence between $(1,1)$ and $(i,t)$ for all $t\le 3$. This corresponds to Thompson's estimator. Finally, panel (F) illustrates the case where observation $(1,1)$ is correlated with all other observations. The dark-to-light shading is meant to imply that the correlation between $(1,1)$ and $(i,t)$ is expected to diminish for $t>1$.
To understand the source of cross-firm and cross-time dependence it is illuminating to consider the linear components model $Y_{it} = \alpha_i + \gamma_t +\varepsilon_{it}$ under the assumption of i.i.d. firm effects $\alpha_i$ and idiosyncratic effects $\varepsilon_{it}$. If the time effect $\gamma_t$ is also i.i.d. then observations $(i,t)$ and $(j,s)$ are independent if $i\ne j$ and $t\ne s$. However, if $\gamma_t$ is serially dependent, then observations $(i,t)$ and $(j,s)$ can be dependent for arbitrary indices.
In most applications the time effect $\gamma_t$ is a proxy for omitted macroeconomic factors, and is therefore unlikely to be i.i.d. or have truncated serial dependence. Most macroeconomic variables have untruncated autocorrelation functions.
We empirically illustrate the importance of this feature. Consider two variables involved in a standard market value equation: log Tobin's average Q (market value divided by the stock of non-R&D assets), and log of the R&D stock (relative to the stock of non-R&D assets). Panel regressions of the former on the latter have been the focus of griliches1981market and many subsequent papers \citep*[e.g.,][]{hall2005market,bloom2013identifying,arora2021knowledge}. Using a panel of 727 firms for the years 1981--2001 from \citet*{bloom2013identifying} (see Appendix (ref) for details) we estimated time effects for each series. Estimated time effects are plotted in the two panels of Figure (ref), with log R&D stock on the left and log Tobin's Q on the right. The graphs reveal considerable serial correlation. Their estimated first-order autocorrelations are 0.425 (with a standard error of 0.003), and 0.467 (with a standard error of 0.007), respectively, which are quite large. Furthermore, the autocorrelations are strong at multiple lags as illustrated in their autocorrelograms, which are displayed in Figure (ref). Together, this means that the time effects $\gamma_t$ for these series have substantial serial dependence, which is not well described by finite $M$-dependence. This implies that the dependence structures assumed by \citet*{CGM2011}, menzel2021bootstrap, and \citet*{Thompson2011} are incorrect, but rather need to be modified to allow for serial correlation of arbitrary order, as we propose in the next section.
As described in the previous section, serially correlated time effects $\gamma_t$ imply that the cross-firm autocorrelations $G_m$ in the variance decomposition (ref) are non-zero at potentially any lag $m$. However, we cannot estimate these correlations well at all lags $m$ for fixed $T$, for the same reasons as arise in time-series variance estimation. Under the assumption that the time effects $\gamma_t$ are strictly stationary and weakly dependent (meaning that the autocorrelation function decays to zero) then it is sufficient to focus on the small lags $m$, using a weighted average of the terms in (ref), with the number of terms increasing with sample size. This motivates the following estimator of $\Omega_{NT}$
where $w(m,M)$ is a weight function and $EVC(\cdot)$ is an eigenvalue correction.\footnote{The $EVC(\cdot)$ replaces any negative eigenvalue by zero to ensure positive semidefiniteness.} Inserted into (ref) we obtain our covariance estimator for $\widehat \beta$. This variance estimator is a generalization of the \citet*{CGM2011}, \citet*{Thompson2011}, and \citet*{newey1987simple} estimators. The covariance matrix estimator $\widehat V_{NT}$ can be multiplied by degree-of-freedom adjustments if desired, but we are unaware of any finite sample justification for a particular choice. The estimator simplifies to that of \citet*{CGM2011} when $M=0$, and to that of \citet*{Thompson2011} when $M=2$ and $w(m,M)=1$. Taking the square root of the diagonal elements of $\widehat V_{NT}$ yields standard errors for the elements of $\widehat\beta$.
The estimator (ref) depends on the choice of weights $w(m,M)$. Standard choices include the uniform (truncated) weights $w(m,M)=1$, and the triangular (Newey-West) weights $w(m,M)=1-m/(M+1)$, the latter popularized by \citet*{newey1987simple} for time-series data. We recommend the triangular weights. One advantage of this choice for time series applications as emphasized by \citet*{newey1987simple} is that this ensures a non-negative variance estimator. Unfortunately, the clustered estimator (ref) is not necessarily non-negative, even for $M=0$, as observed by \citet*{CGM2011}. However, the estimator (ref) is considerably less likely to be negative when the weights are triangular than uniform, which is an important practical advantage.
The variance estimator (ref) critically depends on the number $M$, which is often called the lag truncation parameter. It is useful to note that $M$ does not need to be integer-valued. The choice of $M$ leads to a bias/precision trade-off, with larger values of $M$ leading to less bias in the estimator $\widehat \Omega_{NT}$ of $\Omega_{NT}$, but less precision. In principle it is desirable to use a larger value of $M$ when the errors $U_{it}$ are highly serially correlated, and a smaller value of $M$ otherwise, but the extent of serial correlation is generally unknown, leading to the need for an empirical-based choice of $M$.
In the context of time-series variance estimation Andrews1991 proposed a data-driven choice of $M$ based on minimizing the asymptotic mean square error of the variance estimator, which is equivalent to the expression in (ref). We can therefore apply his method for selection of $M$, treating the time-sums of the regression scores as time-series observations. Andrews' formula depends on the specific choice of weight function; we assume triangular weights.
For $j=1,...,k$, let $X_{jit}$ be the $j$th element of $X_{it}$. Define the time-sums $S_{jt}=\sum_{i=1}^n X_{jit}\widehat U_{it}$ of the regression scores. Fit by least squares the AR(1) equations $S_{jt} = \widehat \rho _j S_{j,t-1}+\widehat e_{jt}$. The Andrews rule\footnote{This is calculated from Andrews' equation (6.4), setting his weights $w_a$ to equal the inverse squared variances of the estimated AR(1) processes, which is appropriate for least squares estimation.} for the lag truncation $M$ is
For the case of a scalar regressor this simplifies to
For an even simpler choice, stock_watson suggested the following rule-of-thumb. Setting $\rho = 0.25$ in the above formula (which occurs in a regression when both the regressor and regression error are AR(1) processes with AR(1) coefficients 0.5) the Andrews rule simplifies to $M = 0.75 \cdot T^{1/3}$. For example, for $T=50$, 100, and 200, respectively, the Stock-Watson rule is $M=2.7$, 3.5, and 4.4, respectively. The Stock-Watson can be used in place of the Andrews rule (ref) if desired.
Consider a panel $\{D_{it} : 1 \le i \le N ; 1 \le t \le T \}$ of random vectors consisting of observed and/or unobserved variables that are relevant to the data generating process of the researcher's interest. For instance, in the linear regression model presented in Section (ref), set $D_{it} = (Y_{it},X_{it}',U_{it})'$. With a Borel-measurable function $f$ (generally unknown to the researcher), we consider the framework of panel dependence in $D_{it}$ generated through the stationary process
where $\alpha_i$, $\gamma_t$, and $\varepsilon_{it}$ are random vectors of arbitrary dimension, with the sequences $\{\alpha_i\}$, $\{\gamma_t\}$ and $\{\varepsilon_{it}\}$ mutually independent, $\alpha_i$ is i.i.d. across $i$, and $\varepsilon_{it}$ is i.i.d. across $(i,t)$.\footnote{These i.i.d. conditions may be relaxed in some ways, but we leave it for future research. The existing literature on two-way clustering explicitly or implicitly assumes these i.i.d. conditions, and we continue to impose them in this paper. Our focus, therefore, is to relax the i.i.d. assumption on the $\gamma_t$ component. One way to relax the i.i.d. conditions on the other components is to impose an MDS-type condition, in which case our main results remain to hold. Another is to allow for spatial mixing provided a researcher has spatial information associated with panel data. Also, see Section (ref).} The sequence $\gamma_t$ is a strictly stationary serially correlated process.
The representation (ref) generalizes the Aldous-Hoover-Kallenberg (AHK) representation Kallenberg2006 which has been widely used for two-way clustering theory. See DDG2019,mackinnon2021wild,menzel2021bootstrap. Indeed, it has been argued that the AHK representation is a natural modelling framework for two-way clustered data \citep*{mackinnon2021wild}.\footnote{\citet*{mackinnon2021wild} state “[a] natural stochastic framework for the regression model with multiway clustered data is that of separately exchangeable random variables.” Since separately exchangeable random variables may be represented by the AHK Kallenberg2006, we make this assertion.} A limitation of the AHK representation is that the time effects $\gamma_t$ are mutually independent. We relax this assumption, by directly assuming that (ref) holds, allowing $\gamma_t$ to be serially dependent. This is a strict generalization of the AHK representation.
In this section we provide an asymptotic distribution theory for the least squares estimator, our proposed covariance matrix estimator, and associated test statistics. We start in Section (ref) by examining a multivariate mean, followed in Section (ref) with linear regression.
In this subsection, we focus on estimation of a multivariate mean. Let $X_{it}$ be an $m\times 1$ random vector satisfying equation (ref) for $D_{it}=X_{it}$. The standard estimator of the population mean $\theta =E[X_{it}]$ is the sample mean $\widehat\theta = (NT)^{-1} \sum_{i=1}^N \sum_{t=1}^T X_{it}$.
To obtain an asymptotic representation for $\widehat\theta$, we use a Hoeffding-type decomposition. For simplicity and without loss of generality assume that $E[X_{it}]=0$. Define the random vectors $a_i=E[ X_{it}\mid \alpha_i]$, $b_t=E[ X_{it}\mid \gamma_t]$, and $e_{it}= X_{it}-a_i-b_t$. This gives rise to the following construct:
This expresses $X_{it}$ as a linear function of a firm effect $a_i$, time effect $b_t$, and error $e_{it}$. However, as (ref) is a derived relationship, the error $e_{it}$ is not (in general) i.i.d. The decomposition has the following properties. The random vectors $a_i$ and $b_t$ are independent since they are each functions of the independent sequences $\{\alpha_i\}$ and $\{\gamma_t\}$. The sequence $\{a_i\}$ is i.i.d., and the sequence $\{b_t\}$ is strictly stationary. By iterated expectations we deduce the following: (1) $a_i$, $b_t$, and $e_{it}$ are each mean zero; (2) $E[e_{jt} \mid \alpha_i]=0$ and $E[e_{is} \mid \gamma_t]=0$ for any $i$, $j$, $t$, and $s$; (3) $E[ a_i e_{jt}']=0$ and $E[ b_t e_{is}']=0$ for any $i$, $j$, $t$, and $s$; (4) the sequences $\{a_i\}$, $\{b_t\}$, and $\{e_{it}\}$ are mutually uncorrelated; (5) conditional on $(\gamma_t,\gamma_s)$, $e_{it}$ and $e_{js}$ are independent for $j \ne i$.
Taking averages we find that
The uncorrelatedness of the sequences implies that the three sums are uncorrelated. Hence the variance of $\widehat\theta$ equals the sum of the variance matrices of the three components. The variance of the first component equals $N^{-1}$ times the variance matrix of $a_i$, the variance of the second component approximately equals $T^{-1}$ times the long-run variance matrix of $b_t$ (since the latter is serially correlated), and the variance of the third component approximately equals $(NT)^{-1}$ times the long-run variance matrix of $e_t$. The technical details are deferred to Appendix (ref).
We now present sufficient conditions for this decomposition to be valid.
Assumption (ref) (i) assumes that the observed random vectors are generated following the nonlinear, nonseparable, factor structure (ref). Assumption (ref) (ii) imposes moment conditions. Assumption (ref) (iii) imposes weak dependence on the time effects. The moment and mixing conditions here are standard in time-series theory, including that of newey1987simple and \citet*{hansen1992consistent}.
Define the variance matrices:
which are independent of $i$ and $t$. Given the decomposition (ref), we can write the variance of the sample mean as a weighted sum of the variance components (ref)-(ref).
A proof is provided in Appendix (ref). Theorem (ref) shows that the asymptotic variance of the sample mean depends on three components, inversely proportional to the number of firms $N$, time dimension $T$, and their product $NT$. When either $\Sigma_a>0$ or $\Sigma_b>0$ (which occurs when there is a non-degenerate firm or time effect) then the third term in the asymptotic variance is of lower stochastic order. However, in the special case where $X_{it}$ is i.i.d., then $\Sigma_a=0$ and $\Sigma_b=0$ so the first two terms equal zero, the long-run variance of $e_{it}$ simplifies to $\Sigma_e=var(X_{it})$, and the variance expression simplifies to $var( \widehat \theta )=var(X_{it})/NT$. Consequently, the rate of convergence of the sample mean depends on the cluster structure.
For our distribution theory we require the following additional condition.
Assumption (ref) (i) requires the presence of at least one-way clustering. This assumption is analogous to the positive definiteness condition in DDG2019 and equation (16) of \citet*{mackinnon2021wild}. Our results will continue to hold even if $N$ and $T$ diverge at different rates, but we make the homogeneous rate assumption for ease of exposition. On the other hand, our results will not hold under fixed $N$ or fixed $T$. Assumption (ref) (ii) is the contrary case of no clustering. As we show below, our results hold under either condition. While Assumption (ref) is sufficient for our results, it is probably stronger than necessary, but is used for its simplicity and tractability. Assumption (ref) does rule out possible scenarios, including cases which lead to non-Gaussian limit distributions menzel2021bootstrap -- see our discussion in Section (ref). The non-Gaussian cases can be characterized by daga generating processes that consist of degenerate additive factors of $i$, degenerate additive factors of $t$, and a small number of interactive factors between $i$ and $t$ chiang2022standard. With this said, it is legitimate to rule out non-Gaussian cases as our focus is on standard errors (as in the title of this article), which would not make sense under non-Gaussian limit distributions.
A proof is provided in Appendix (ref). Theorem (ref) shows that the self-normalized sample mean is asymptotically normal. Self-normalization is used to allow for differing rates of convergence due to clustering structure.
In this section, we focus on the setting where there is precisely one unit of observation per cluster intersection. In applications, there may be heterogeneous per-cluster numbers of observations, e.g., unbalanced panels. The above theory will straightforwardly extend to such cases -- see Appendix (ref).
We now revisit the linear panel regression model (ref)--(ref) and the OLS estimator (ref). In this subsection we set $D_{it} = (Y_{it},X_{it}',U_{it})'$ , so that the framework (ref) is
Set $a_i=E[X_{it}U_{it}\mid \alpha_i]$ and $b_t=E[X_{it}U_{it}\mid \gamma_t]$. Let $\Sigma_a$ and $\Sigma_b$ be the variance/long-run variance matrices of $a_i$ and $b_t$, respectively.
Assumptions (ref) (i)--(v) above are the counterparts of Assumptions (ref) and (ref), extended to the regression model. The moment and mixing conditions are standard in time series regression. Assumptions (ref) (v)--(vi) are needed for consistent variance estimation. The $\alpha$-mixing condition of Assumption (ref) has been strengthened to $\beta$-mixing in Assumption (ref). This is because our proof of consistent variance estimation relies on a deep fourth-order summability result due to yoshihara1976limiting which relies on $\beta$-mixing.
Under these assumptions, the asymptotic variance of $\widehat \beta$ is
where $\Omega_{NT}$ is defined in (ref)-(ref). Our proposed estimator of $V_{NT}$ is (ref) with (ref).
For our theory we focus on a vector-valued parameter $\theta=R'\beta$ for some $k \times m$ matrix $R$. This includes individual coefficients when $m=1$. The estimator of $\theta$ is $\widehat \theta=R'\widehat \beta$, its asymptotic variance is $\Sigma_{NT}=R'V_{NT}R$, with estimator $\widehat\Sigma_{NT}=R'\widehat V_{NT}R$. For the case $m=1$, a standard error for $\widehat \theta$ is $\widehat\sigma_{NT}=\sqrt{R'\widehat V_{NT}R}$. We now establish consistency of our variance estimator
A proof is provided in Appendix (ref). It relies on some technical lemmas in Appendix (ref) that are of potential independent interest. This result shows that our proposed variance estimator is consistent. Notice that we state consistency as a self-normalized matrix ratio. This allow for the differing rates of convergence covered by Assumption (ref).
Theorem (ref) is new. It is the first demonstration of consistent variance estimation under two-way clustering with serially dependent time effects.
The proof of Theorem (ref) includes some technical innovations. Of particular note is the use of the fourth-order summability condition of yoshihara1976limiting, combined with projections on the individual-specific and time-specific factors. The summability condition is needed to calculate the variance of the variance estimator, which is a fourth-order sum.
A proof is provided in Appendix (ref).
Theorem (ref) shows that the least squares estimator $\widehat \theta$ is asymptotically normal. Normality holds when the estimator is standardized by its asymptotic covariance matrix $\Sigma_{NT}$ or by its estimator $\widehat\Sigma_{NT}$. This latter result shows (for the case $m=1$) that t-ratios constructed with our standard errors are asymptotically $N(0,1)$. It also shows (for the case $m>1$) that Wald-type tests constructed with our covariance matrix estimator are asymptotically $\chi_m^2$. Hence conventional inference methods can be used with the least squares estimator $\widehat \beta$, our variance estimator $\widehat V_{NT}$, and our standard errors.
Theorem (ref) is new. It is the first result which rigorously demonstrates asymptotic normality of least squares estimators and t-ratios under two-way clustering with serially correlated time effects. The asymptotic normality presented here adapts to the unknown convergence rate. It is, however, worthy to note that this result is pointwise in DGP. In the absence of correlated time effects, menzel2021bootstrap discusses the issues of uniform inference. In a linear panel data context, lu2022uniform provides uniformly valid inference procedure. Although it remains unclear to us whether their approaches can be adapted to our framework, it is an interesting future research avenue to investigate the potential uniformity properties of the asymptotics under various sets of sequences of DGPs.
In contrast, test statistics constructed with the popular CGM variance estimator will not have conventional asymptotic distributions when the time effects $\gamma_t$ are serially correlated. The CGM variance estimator is inconsistent in this situation, so test statistics will have distorted asymptotic distributions.
For panel data, researchers often include one-way or two-way fixed effects. This section has two contents regarding fixed-effect models. First, Section (ref) argues that two-way clustering is still necessary in general even if a researcher includes two-way fixed effects. Second, Section (ref) extends our theory to two-way fixed-effect regressions.
Some empirical economists believe that it is unnecessary to cluster standard errors if fixed effects are included in estimation. In this section, we argue that fixed effects will not generally solve the problem of two-way cluster dependence.
Consider the two-way fixed-effect model:
Note that $X_{it}$ and $U_{it}$ are generated by the latent variables $\alpha_i=(\alpha_{i0},\alpha_{i1},\alpha_{i2},\alpha_{i3})$, $\gamma_t = (\gamma_{t0},\gamma_{t1},\gamma_{t2},\gamma_{t3})$ and $\varepsilon_{it} = (\varepsilon_{it0},\varepsilon_{it1})$. Here, $\alpha_{i0}$ and $\gamma_{t0}$ are additive fixed effects. Suppose that $\alpha_{i0}$, $\alpha_{i1}$, $\alpha_{i2}$, $\alpha_{i3}$, $\gamma_{t0}$, $\gamma_{t1}$, $\gamma_{t2}$, $\gamma_{t3}$, $\varepsilon_{it0}$, and $\varepsilon_{it1}$ are mutually independent with mean 0 and variance 1.
To abstract away from finite-sample issues, consider the population double differences\footnote{For a formal account of the discrepancy between the population and sample double differences, see Section (ref).}
Note that they reduce to
under the example (ref). Certainly, the double differencing removes the FEs, $\alpha_{i0}$ and $\gamma_{t0}$, but still leaves the strong two-way dependence through $(\alpha_{i1},\alpha_{i2},\alpha_{i3})$ and $(\gamma_{t1},\gamma_{t2},\gamma_{t3})$. Observe that there is no endogeneity, as $ E[\widetilde U_{it} | \widetilde X_{it}]=0. $ Furthermore, the score
entails the non-degenerate projections
with respect to $\alpha_i$ and $\gamma_t$, respectively.
Hence, in this example, the components of $\alpha_i$ and $\gamma_t$ are not eliminated by the two-way fixed effects regression of $Y$ on $X$, so the score $\ddot X_{it} \ddot U_{it}$ is still two-way dependent, and two-way clustering is necessary for calculation of the covariance matrix. Furthermore, the score is not degenerate so our theory to be presented in Section (ref) applies.
This section shows that our standard errors extend to two-way fixed-effect models. In different settings, the existing literature has considered inference for two-way fixed-effect models. verdier2020estimation studies linear regression with two-way fixed effects for fixed $T$ for sparsely matched data. juodis2021shock considers bootstrap-based inference for linear models with two-way fixed effects under large $N$ and $T$ with a different set of assumptions.
Consider the two-way fixed-effect model
where $\xi_i$ and $\eta_t$ are fixed effects and $E[U_{it}]=0$.\footnote{Note that (ref) implicitly requires that $\xi_i = f_2(\alpha_i)$ and $\eta_t = f_3(\gamma_t)$.} Define the within-transformed outcome and within-transformed regressors by
respectively. The within transformations induce complex dependence structure between transformed variables. Using idempotency of the within transformation matrices (see Ch. 17.8 of hansen2021econometrics), the two-way within estimator $\widehat \beta$ for $\beta$ is defined as the OLS estimator of $\ddot Y_{it}$ on $\ddot X_{it}$ and thus satisfies
Also, define the variance estimator $\widehat \Sigma_{NT}$ as in (ref) with $(\ddot Y_{it},\ddot X_{it}')$ in place of $(Y_{it},X_{it}')$ and with the two-way within estimator $\widehat \beta$ defined in this section. Denote the population counterpart of the within-transformed regressor by
Note that $(\widetilde X_{it}',U_{it})$ only depends on $\alpha_i,\gamma_t$, and $\varepsilon_{it}$ and satisfies $E[\widetilde X_{it}]=0$.
A proof is provided in Appendix (ref). It is worthy noting that the proof is not a mere extension of the previous theorems, but requires nontrivial technicalities involving maximal inequalities in the time-series framework.
In this section, we use simulated data to examine the performance of our proposed robust variance estimator in comparison with six existing alternatives.
We generate data based on the linear model
where the right-hand side variables $(X_{it},U_{it})'$ are generated through the panel dependence structure
We set $(\beta_0,\beta_1)= (1,1)$ throughout. For the weight parameters, we use $(w_\alpha,w_\gamma,w_\varepsilon) = (0,0,1)$ to generate i.i.d. data and also use $(w_\alpha,w_\gamma,w_\varepsilon) = (0.25,0.50,0.25)$ to generate dependent data. The latent components $(\alpha^x_i,\alpha^u_i,\varepsilon^x_{it},\varepsilon^u_{it})$ are all mutually independent $N(0,1)$.
The latent common time effects $(\gamma^x_t,\gamma^u_t)$ are dynamically generated according to the AR(1) design:
The initial values are drawn from $N(0,1)$. We vary the AR coefficient $\rho \in \{0.25,0.50,0.75\}$ across sets of simulations.
For each realization of observed data $\{(Y_{it},X_{it}) : 1 \le i \le N, 1 \le t \le T\}$ constructed according to the data generating process described above, we estimate $(\beta_0,\beta_1)$ by OLS. Our objective is to evaluate the performance of our proposed robust variance estimator $\widehat V_{NT} = \widehat Q^{-1} \widehat\Omega_{NT} \widehat Q^{-1} $, where $\widehat Q$ and $ \widehat\Omega_{NT} $ are given in (ref) and (ref), and the tuning parameter $\widehat M$ is chosen according to the rule (ref). Through simulation studies, we examine the performance of this robust variance estimator (hereafter referred to as CHS; Chiang-Hansen-Sasaki) in comparison with six existing alternative variance estimators which are in popular use for panel data analysis. They include the heteroskedasticity robust estimator (EHW; Eicker-Huber-White; also known as HC0), the cluster robust estimator within $i$ (CR$i$; which corresponds to the Liang-Zeger-Arellano estimator), the cluster robust estimator within $t$ (CR$t$), the two-way cluster robust estimator \citep*[CGM;][]{CGM2011}, the wild bootstrap estimator \citep*[MNW;][]{mackinnon2021wild} for CGM, the bootstrap estimator \citep*[M;][]{menzel2021bootstrap},\footnote{Implementation of M requires some tuning parameters. We set the number of bootstrap iterations and the model selection tuning parameters, $\kappa_a$ and $\kappa_g$ following the simulation code for regressions by menzel2021bootstrap. In addition to the default method of M, we also ran M without its model selection feature to find its results the same as those of the default method. Hence, we only report results by the default method of M.} and the two-way cluster robust estimator with 2-dependence \citep*[T;][]{Thompson2011}.
Table (ref) reports simulation results. Reported values are the coverage frequencies for the slope parameter $\beta_1$ for the nominal probability of 95% based on 10,000 Monte Carlo iterations. The top and bottom panels show coverage probability results under the i.i.d. design and the dependence design, respectively. In each group of three consecutive rows, the panel sample sizes $(N,T)$ vary by rows. Cells are shaded based on the proximity of the simulated coverage probability to the nominal probability of 0.95; the darker shades indicate more correct coverage.
We observe the following four points in these results. First, under the i.i.d. design (rows (I)--(III)), EHW, CR$i$, CR$t$, CGM, MNW and CHS produce accurate coverage probabilities. On the other hand, M yields over-coverage consistently across different sample sizes, and T yields under-coverage especially under small $T$. Second, EHW and CR$i$ tend to behave poorly in general when there are two ways of cluster dependence as in rows (IV)--(XII). Third, CR$t$ also tends to behave poorly under larger extents of serial dependence as in rows (X)--(XII). Likewise, CGM, MNW, and M perform less preferably as the serial dependence becomes even stronger, as in rows (X)--(XII). These results are consistent with the fact that these methods do not account for serially correlated common time effects. Fourth, in contrast, T and CHS behave more robustly under strong serial dependence especially when $T$ is large, as in rows (X) and (XI). Whenever $T$ is small, as in rows (III), (VI), (IX) and (XII), however, T incurs severe under-coverage and hence CHS outperforms T in general. The last observation that T performs poorly for small sample sizes is consistent with the similar observations made by Thompson2011 in his Monte Carlo simulation studies. In summary, we demonstrate that confidence intervals constructed with our proposed standard errors lead to robustly superior coverage performance relative to the existing methods.
We ran additional simulations beyond those presented in this section. Their results are found in Appendix (ref). Specifically, Appendix (ref) illustrates power analyses, and Appendix (ref) presents simulations for the two-way fixed-effect estimator.
In this section, we highlight differences across the alternative standard error estimates for estimates of a simple asset pricing model. Consider the Fama-French three-factor model
where $R_{it}$ is the total return of portfolio/stock $i$ in month $t$, $R_{ft}$ is the risk-free rate of return in month $t$, $R_{Mt}$ is the total market portfolio return in month $t$, $SMB_t$ is the size premium (small$-$big), $HML_t$ is the value premium (high$-$low), and $\beta = (\beta_1,\beta_2,\beta_3)'$ are the factor coefficients.
We use two data sets of portfolio/stock returns. They are (A) 44 industry portfolios excluding four financial sectors (banking, insurance, real estate, and trading) and (B) individual stocks. For each of these data sets, (A) and (B), we use the monthly panel of length $120$ from January 2000 to December 2009. For the individual stock data set (B), we use the balanced portion of the panel data, consisting of $N=779$ stocks. The risk-free rate is based on the monthly 30-day T-bill beginning-of-month yield. See Appendix (ref) for the source of data.
Let $\ddot Y_{it}$ and $\ddot X_{it}$ denote the within-transformations\footnote{ Technically, our theory does not allow for these demeaned variables as they are $(N,T)$-dependent. But we expect that there should be no change to the theory if we allow for array data. Formalizing this would be a useful future research direction. } of $R_{it}-R_{ft}$ and $(R_{Mt}-R_{ft},SMB_t,HML_t)'$, respectively. ($\ddot X_{it}$ is homogeneous in the cross section.) We estimate $\beta$ by the within-estimator
to remove any additive portfolio/stock fixed effects. Thus, our proposed standard errors are computed based on $\widehat V_{NT}=\widehat Q^{-1}\widehat \Omega_{NT}\widehat Q^{-1}$ with $\widehat Q$ and $\widehat \Omega_{NT}$ from (ref) and (ref) with $Y_{it}$, $X_{it}$ and $\widehat U_{it}$ replaced by $\ddot Y_{it}$, $\ddot X_{it}$ and $\ddot Y_{it} - \ddot X_{it} \widehat\beta$, respectively. Table (ref) summarizes estimates $\widehat\beta$ of $\beta$ along with alternative standard error estimates of them for each of the two data sets described above.
For each row of Table (ref), the standard error due to the Eicker-Huber-White (EHW) estimator is smaller than any other standard error. On the other hand, the two-way cluster robust estimators of Cameron-Gelbach-Miller (CGM) and Thompson (T) and our proposed estimator (CHS) tend to yield the largest standard errors in each row. The remaining three standard errors stay in the middle between these two groups. The estimator by Menzel (M) behaves similarly to EHW in panel (A) while it behaves similarly to CGM in panel (B). This puzzling outcome arises from the model selection feature of M. In fact, M would also behave similarly to CGM in panel (A) as well if the tuning parameter of M were chosen to take a much smaller value.\footnote{As in the simulation section, we set the number of bootstrap iterations and the model selection tuning parameters for M following the simulation code for regressions by menzel2021bootstrap.} Because of these idiosyncratic behaviors of M that depend on discrete outcomes of model selection which in turn depend on tuning parameters, we will hereafter focus on the other estimators in comparing the results.
Observe in panel (A) that the statistical significance of the coefficient $\beta_2$ of SMB meaningfully diminishes as the standard error estimator becomes more robust (again, except for M). Specifically, it is significant at the 5% level with EHW, CR$i$ and M, but it becomes insignificant at this level with CR$t$, CGM, T, or CHS. Furthermore, while it is significant at the 10% level with EHW, CR$i$, CR$t$, and CGM, it becomes insignificant at this level with T and CHS. This part of the estimation results shows a case where accounting for serial correlation in common time effects may even overturn conclusions from statistical inference based on the other standard errors. Accounting for arbitrary untruncated time correlation, however, CHS yields a slightly smaller standard error estimate than T for this case.
In this paper, we propose new robust standard error estimators for panel data. The new estimators account for the cluster dependence within $i$, the cluster dependence within $t$, and serial dependence in the common time effects. In particular, all the existing robust standard error estimators fail to accommodate untruncated serial dependence in the common time effects, while this feature is relevant to empirical data used in economics and finance. Simulation studies show that the new standard errors produce robustly superior coverage performance than existing alternatives, including the heteroskedasticity robust estimator (Eicker-Huber-White; also known as HC0), the cluster robust estimator within $i$, the cluster robust estimator within $t$, the two-way cluster robust estimator \citep*{CGM2011}, the bootstrap estimator \citep*{menzel2021bootstrap}, and the two-way cluster robust estimator with 2-dependence \citep*{Thompson2011}.
In the rest of this section, we discuss limitations of our method and potentials of future research in relation to the existing literature. Since the seminal work by \citet*{CGM2011} and \citet*{Thompson2011}, a few important papers have proposed methods of robust inference in two way cluster dependence.
DDG2019 derive Donsker results under multiway cluster dependence. In the current paper, we only derive limit distributions for finite-dimensional parameters that are relevant to many empirical applications in economics and finance. In other words, DDG2019 provide a generalization of existing results by allowing for empirical processes, we on the other hand provide a generalization in a different direction by allowing for serial dependence in common time effects. Combining these two directions of generalization is left for future research.
menzel2021bootstrap proposes a method of (conservative) inference with uniform validity over a large class of distributions including the case of non-Gaussian degeneracy under two-way cluster dependence. In the current paper, for the purpose of providing a simple method of inference via analytic standard error formulas, we focus on the case of Gaussian degeneracy as well as non-degenerate cases. In other words, menzel2021bootstrap provides a generalization of existing results by allowing for uniformity over a large class, we on the other hand provide a generalization in a different direction by allowing for serial dependence in common time effects. Combining these two directions of generalization is also left for future research.
chiang2020inference derive a high-dimensional central limit theorem under multiway cluster dependence. In the current paper, we only consider finite-dimensional parameters that are relevant to many empirical applications in economics and finance. In other words, chiang2020inference provide a generalization of existing results by allowing for high dimensionality, we on the other hand provide a generalization in a different direction by allowing for serial dependence in common time effects. Again, combining these two directions of generalization is left for future research.
The independence conditions in Assumption (ref) (i) may be relaxed in a couple of directions. One way is relax the i.i.d. assumption on the $\varepsilon_{it}$ factor in (ref). With this said, the existing literature explicitly or implicitly makes this assumption, and we continue to focus on i.i.d. $\varepsilon_{it}$. Another way is to relax the i.i.d. assumption on the $\alpha_{i}$ factor in (ref). For instance, if a researcher obtains spatial information associated with panel data, then it may be a possibility to allow for spatial $\alpha$-mixing jenish2009central. We leave these extensions for future research.
Our standard errors allowing for serial correlation in time effects effectively use the Newey-West-type long-run variance estimation, but it is well known that in some cases a relatively long time series may be required for such estimators to perform well lazarus2018har. Inference based on moving-block bootstrap gonccalves2011moving may improve the finite-sample performance, and we suggest it as another direction for future research. Another promising recent proposal by chen2023fixed is to use fixed-b asymptotic theory to derive bias corrections and improve the distributional approximation.