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.
56,301 characters · 10 sections · 46 citation commands
Standard Errors for Panel Data Models with Unknown Clusters
\thispagestyle{empty}
\onehalfspacing
\setcounter{page}{1} \pagenumbering{arabic}
Consider a linear panel regression with fixed-effects:
where $\alpha_{i}$ and $\mu_{t}$ are individual fixed-effects and time fixed effects; $x_{it}$ is a $k \times1$ vector of explanatory variables; $u_{it}$ is an unobservable error component. The outcome variable $y_{it}$ and fixed effects are scalars, and $\beta$ is a $k\times 1$ vector.
This paper is about the standard error of the fixed-effect ordinary least squares (OLS). One of the commonly used standard errors for OLS in empirical research is the white1980heteroskedasticity heteroskedasticity robust standard error in the cross-sectional setting. In the presence of serial and cross-sectional correlations, the conventional panel standard errors may be biased. newey1986simple introduced heteroskedasticity and autocorrelation consistent (HAC) covariance matrix estimator for time series, which allows serial correlations (also see andrews1991, newey1994). The cluster standard errors suggested by arellano1987 are often reported in studies of the panel model. This estimator is robust to heteroskedasticity in the cross-section and also arbitrary serial correlation, but it focuses on the large-$N$ small-$T$ scenario. The case of large-$N$ large-$T$ is then studied by ahn2014, hansen2007, among many others, while either cross-sectional or serial independence is required. hansen2007 examined the covariance estimator when the time series dependence is left unrestricted. In addition, vogelsang2012 studied the asymptotic theory that is robust to heteroskedasticity, autocorrelation, and spatial correlation, which extended and generalized the asymptotic results of hansen2007 for the conventional cluster standard errors including time fixed effects. stock2008 suggested a bias-adjusted heteroskedasticity-robust variance matrix estimator that handles serial correlations under any sequences of $N$ or $T$. Also, see petersen2009 who used a simulation study to examine different types of standard errors, including the clustered, Fama-MacBeth, and the modified version of Newey-West standard errors for panel data. In general, on the other hand, the conventional cluster standard errors assume that individuals across clusters are independent. Also, the cluster structure should be known such as schools, villages, industries, or states. See arellano2003, cameron2015 and greene2003. However, the knowledge of clusters is not available in many applications.
In a recent interesting paper, abadie2017should argue that clustering is an issue more of sampling design or experimental design. Clustered standard errors are not always necessary and researchers should be more thoughtful when applying them. One reason is that clustering may result in an unnecessarily wider confidence interval. Clustered standard errors are derived from the modelling perspective (model implied variance matrix) and are widely practiced, see, for example, angrist2008mostly, cameron2005microeconometrics, and wooldridge2003cluster, wooldridge2010econometric. In this paper, we continue to take the modeling perspective. Because of our use of thresholding method, the resulting confidence interval is not necessarily much wider, even if all cross-sectional units are allowed to be correlated. Furthermore, the proposed approach is also applicable when the knowledge of clustering is not available.
We provide a robust standard error that allows both serial and cross-sectional correlations. We do not impose parametric structures on the serial or cross-sectional correlations. We assume these correlations are weak and apply nonparametric methods to estimate the standard errors. To control for the autocorrelation in time series, we employ the Newey-West truncation. To control for the cross-sectional correlation, we assume sparsity for cross-section $(i,j)$ pairs, potentially resulting from the presence of cross-sectional clusters, but the knowledge on clustering (the number of clusters and the size of each cluster) is not assumed. We then estimate them by applying the thresholding approach of bickel2008a. We also show how to make use of information on clustering when available. In passing we point out that instead of robust standard errors, in a separate study, baichoiliao2019 proposed a feasible GLS (FGLS) method to take into account heteroskedasticity and both serial and cross-sectional correlations. The FGLS is more efficient than OLS.
The methods we employ in this paper, banding and thresholding, are regularization methods, and have been used extensively in the recent machine learning literature for estimating high-dimensional parameters. Nonparametric machine learning techniques have been proved to be useful tools in econometric studies.
The rest of the paper is organized as follows. In Section (ref), we describe the models and standard errors as well as the asymptotic results of OLS. Monte Carlo studies evaluating the finite sample performance of the estimators are presented in Section (ref). Section (ref) illustrates our methods in an application of US divorce law reform effects. Conclusions are provided in Section (ref) and all proofs are given in Appendix (ref).
Throughout this paper, $\nu_{min}(A)$ and $\nu_{max}(A)$ denote the minimum and maximum eigenvalues of matrix $A$. We use $\|A\| = \sqrt{\nu_{max}(A'A)}$, $\|A\|_{1} = max_{i}\sum_{j}|A_{ij}|$ and $\|A\|_{F} = \sqrt{tr(A'A)}$ as the operator norm, the $\ell_1$-norm and the Frobenius norm of a matrix A, respectively. Note that if $A$ is a vector, $\|A\|$ is the Euclidean norm, and $|a|$ is the Absolute-value norm of a scalar $a$.
We consider the following model:
where $\beta$ is a $k \times 1$ vector of unknown coefficients, $x_{it}$ is a $k \times 1$ vector of regressors, and $u_{it}$ represents the error term, often known as the idiosyncratic component. This formulation incorporates the standard fixed effects models as in hansen2007. For example, $x_{it}, y_{it}$ and $u_{it}$ can be interpreted as variables resulting from removing the nuisance parameters from the equation, such as first-differencing to remove the fixed effects. Indeed, it is straightforward to allow additive fixed effects by using the usual demean procedure.
For a fixed $t$, model ((ref)) can be written as:
where $y_{t} = (y_{1t}, ..., y_{Nt})'$ $(N \times 1)$, $x_{t} = (x_{1t}, ..., x_{Nt})'$ $(N \times k),$ and $u_{t} = (u_{1t}, ..., u_{Nt})'$ $(N \times 1)$. To economize notation, we define $y_{i} = (y_{i1}, ..., y_{iT})'$ $(T \times 1)$, $x_{i} = (x_{i1}, ..., x_{iT})'$ $(T \times k),$ and $u_{i} = (u_{i1}, ..., u_{iT})'$ $(T \times 1)$. So when the vector $y$ is indexed by $t$, it refers to an $N\times 1$ vector, and when $y$ is indexed by $i$ it refers to a $T\times 1$ vector. Similar meaning is applied to $x$ and $u$. There is no confusion when context is clear.
The (pooled) ordinary least square (OLS) estimator of $\beta$ from equations ((ref)) and ((ref)) may be defined as
The variance of $\widehat{\beta}$ depends on both $V_{X} \equiv \cfrac{1}{NT}\sum\limits_{i=1}^{N}\sum\limits_{t=1}^{T}x_{it}x_{it}'$, and particularly,
The goal of this paper is to consistently estimate $V$ in the presence of both serial and cross-sectional correlations in $\{u_{it}\}$.
There are two types of clustered standard errors suggested by arellano1987. The original individual clustered version is
with $\widehat{u}_{i} = y_{i}-x_{i}\widehat{\beta} $ are the OLS residuals, and this estimator allows for arbitrary serial dependence and heteroskedasticity within individuals. In addition, $\widehat{V}_{CX}$ assumes no cross-section correlation.
The time-clustered version, which allows for heteroskedasticity and arbitrary cross-sectional correlation, is
with $\widehat{u}_{t} = y_{t}-x_{t}\widehat{\beta}$. Here $\widehat{V}_{CT}$ assumes no serial correlation.
The above clustered standard errors are robust to either arbitrary serial correlation or arbitrary cross-sectional correlation, respectively. In practice, however, since the dependence assumption is unknown, an over-rejection problem may occur. Specifically, if there exist both serial and cross-sectional correlations, these estimators are not robust anymore, as our numerical evidence shows in Section (ref) (e.g., Tables (ref) and (ref)).
To control for the serial correlation, a simple modification of $\widehat{V}_{CT}$ using newey1986simple is
where $\omega(\cdot)$ is the kernel function and $L$ is the bandwidth. This estimator is suggested by driscoll1998. When $N$ is large, however, ((ref)) accumulates a large number of cross-sectional estimation noises.
More generally, let
Then equation ((ref)) can be written as
Unlike time series observations, cross-sectional observations have no natural ordering. They can be arranged in different orders. That is why cross-sectional correlation is more difficult to control. The usual cluster standard error makes the following assumption: let $C_1,...,C_G$ be disjoint subsets of $\{1,...,N\}$, so that they are known clusters and that $V_{ij}=0$ when $i$ and $j$ belong to different clusters. So $V$ can be expressed as $$ V= \frac{1}{N}\sum_{g=1}^G\sum_{(i,j)\in C_g} V_{ij}. $$ See liang1986. Suppose the cardinality of each $C_g$ is small (this would be the case if the number of clusters $G$ is large) or grows slowly with $N$, then we only need to estimate $\sum_{g=1}^G\sum_{(i,j)\in C_g}1$ number of $V_{ij}$'s, greatly reducing the number of pair-wise covariances. But as commented in the literature, this requires the knowledge of $C_1,...,C_G$, which in some applications, is not naturally available.
The key assumption we make is that conditionally on $x_{t}, \{u_{it}\}$ is weakly correlated across both $t$ and $i$. Essentially, this means $ V_{ij}$ is zero or nearly so for most pairs of $(i, j)$. There is a partition $\{(i,j) : i,j \leq N\} = S_{s} \bigcup S_{l}$ so that
where the subscript “$s$" indicates “small", and “$l$" indicates “large". We assume that $(i,i) \in S_l$ for all $i\leq N$, and importantly, most pairs $(i,j)$ belong to $S_s$. Yet, we do not need to know which elements belong to $S_{s}$ or $S_{l}$. Then
Furthermore, let $\omega(h, L) = 1-h/(L+1)$ be the Bartlett kernel. Also see andrews1991 for other kernel functions. As suggested by newey1986simple, $V_{ij}$ can be approximated by
Then approximately,
The above approximation plays the fundamental role of our standard error estimator. We estimate $V_{ij}$ using newey1986simple, and estimate $S_l$ using the cross-sectional thresholding.
To apply newey1986simple, we estimate $V_{u,ij}$ by
where $\widehat{u}_{it} = y_{it}-x_{it}'\widehat{\beta}$. For a predetermined threshold value $\lambda_{ij}$, we approximate $S_l$ by
Hence, a “matrix hard-thresholding" estimator of $V$ is
As for the threshold value, we specify
for a constant $M>0$. The converging sequence $\omega_{NT} \rightarrow 0$ is chosen to satisfy:
In practice, the thresholding constant, $M$, can be chosen through multifold cross-validation, which is discussed in the next subsection. In addition, we can obtain $\widehat{V}_{DK}$ from $\widehat{V}_{\text{Hard}}$ by setting $M=0$.
We also recommend a “matrix soft-thresholding" estimator as follow:
where $\widehat{S}_{u,ij}$ is
where the $(k,k')$'s element of $A_{u,ij}$ is ($sgn(x)$ denotes the sign function)
for the threshold value
for some constant $M>0$.\\
Our suggested estimators, $\widehat{V}_{\text{Hard}}$ and $\widehat{V}_{\text{Soft}}$, require the choice of tuning parameters $L$ and $M$, which are the bandwidth and the threshold constant respectively. To choose the bandwidth $L$, we recommend using $L = 4(T/100)^{2/9}$ as newey1994 suggested.
In practice, $M$ can be chosen through multifold cross-validation. After obtaining the estimated residuals $\widehat{u}_{it}$ by OLS, we split the data into two subsets, denoted by $\{\widehat{u}_{it}\}_{t \in J_{1}}$ and $\{\widehat{u}_{it}\}_{t \in J_{2}}$; let $T(J_{1})$ and $T(J_{2})$ be the sizes of $J_{1}$ and $J_{2}$, which are $T(J_{1})+ T(J_{2}) = T$ and $T(J_{1})\asymp T$. As suggested by bickel2008a, we can set $T(J_{1})=T (1-\log(T)^{-1})$ and $T(J_{2})=T/\log(T)$; $J_1$ represents the training data set, and $J_2$ represents the validation data set.
The procedure requires splitting the data multiple times, say $P$ times. At the $p$th split, we denote by $\widehat{V}^{p}$ the sample covariance matrix based on the validation set, defined by \[ \widehat{V}^{p} = \frac{1}{N}\sum_{ij}S_{u,ij}^{p}, \] with $S_{u,ij}^{p}$ defined similarly to $S_{u,ij}$ using data on $J_2$. Let $\widehat{V}_{s}(M)$ be the thresholding estimator with threshold constant $M$ using the entire sample. Then we choose the constant $M^{*}$ by minimizing a cross-validation objective function
and the resulting estimator is $\widehat{V}_{s}(M^*)$. We use $L=4(T/100)^{2/9}$ for both $\widehat{V}_{s}(M)$ and $\widehat{V}^{p}$ and find that setting $M_0=1$ works well. So the minimization is taken over $M\in (0,1)$ through a grid search.
The above procedure modifies that of bickel2008a in two aspects. One is to use the entire sample when computing $\widehat{V}_{s}$ instead of $J_1$. Since $T(J_{1})$ is close to $T$, this modification does not change the result much, but simplifies the computation. The second modification is to use a consecutive block for the validation set because of time series, so that the serial correlation is not perturbed. Hence in view of the time series nature, we first divide the data into $P=\log(T)$ blocks with block length $T/\log(T)$. Each $ J_{2} $ is taken as one of the $P$ blocks when computing $\widehat V^p$, similar to the K-fold cross validation. We have conducted simulations of the cross-validation in the presence of both correlations, and the results show that this procedure performs well. For instance, the cross-validation tends to choose smaller $M$ as the cross-sectional correlation becomes stronger. Due to the page limit, however, those are not reported in this paper.
Below we present assumptions under which $\widehat{V}$ (either $\widehat{V}_{\text{Hard}}$ or $ \widehat{V}_{\text{Soft}}$) consistently estimates $V$. We define
and
where $X=\{x_{it}\}_{i\leq N,t\leq T}$. These coefficients give measures of autocovariances and cross-section covariances.
Assumption (ref) (i) is the standard alpha-mixing condition, adapted to the large-$N$ panel. Condition (ii) is new here. It requires weak cross-sectional correlations. It is similar to the “approximate sparse assumption" in bickel2008a. Note that we actually allow the presence of many “small" but nonzero $\|Ex_{it}u_{it}u_{j,t+h}x_{j,t+h}'\|$. Clusters that have “large" $\|Ex_{it}u_{it}u_{j,t+h}x_{j,t+h}'\|$ are unknown to us. Hence the appealing feature of our method is that we allow for unknown clusters.
Essentially the assumption $\omega_{N,T}^{1-q}\max_{i\leq N}\sum_{j \leq N}(\sum_{h=0}^{L}\rho_{ij,h})^{q} = o(1)$ controls the order of elements in $S_l$. The following example presents a case of cross sectional weak correlations that satisfies condition (ii).
Condition (i) is well satisfied by various kernels for the HAC-type estimator. Condition (ii) ensures the Bernstein-type inequality for weakly dependent data. Note that it requires the underlying distributions to be thin-tailed. Allowing for heavy-tailed distributions is also an important issue. However, it would require a very different estimation method, and is out of the scope of this paper. Nevertheless, we have conducted simulation studies under heavy-tailed distributions (e.g., $t$-distribution with degree of freedom 5). Indeed, the proposed estimator works well in this case, even though the theory requires thin-tailed distributions.\footnote{The simulation results for the heavy-tailed distributions are available upon request from the authors.}
We have the following main theorem and all proofs are contained in Appendix A1.\\
Theorem (ref) allows us to construct a $(1-\tau)\%$ confidence interval for $c'\beta$ for any given $c \in \mathbb{R}^{k}$. The standard error of $c'\hat \beta_{OLS}$ is \[ \Big( \frac 1 {NT} c'(V_{X}^{-1}\widehat{V}V_{X}^{-1})c \Big)^{1/2} \] and the confidence interval for $c'\beta$ is $[c'\widehat{\beta} \pm Z_{\tau}\hat{\sigma}/\sqrt{NT}]$ where $Z_{\tau}$ is the $(1-\tau)\%$ quantile of standard normal distribution and $\hat{\sigma} = (c'(V_{X}^{-1}\widehat{V}V_{X}^{-1})c)^{1/2}$.
In this section we examine the finite sample performance of the robust standard errors using simulation study. The data generating process (DGP) used for the simulation is produced by the fixed effect linear regression model
where the true $\beta_{0} =1$. The DGP allows for serial and cross-sectional correlations in $x_{it}$ as follow:
where the constants $\{a_{i}, b_{i}\}_{i=1}^{N}$ are i.i.d. Uniform$(0,\gamma_{X})$, which introduce cross-sectional correlation. In addition, $\nu_{it}$ is modeled as $AR(1)$ process with the autoregressive parameter $\rho_{X}$. Throughout this simulation study, we set $\rho_X = 0.3$ and $\gamma_{X}=1$.
We generate the error terms, $u_{it}$, in three different cases as follow:
The regressor is uncorrelated with the error term $u_{it}$ each other. In Case 1, we generate the error term similar to $x_{it}$. The constants $\{c_{i}, d_{i}\}_{i=1}^{N}$ are i.i.d. Uniform$(0,\gamma)$, which introduce cross-sectional correlation, and heteroskedasticity when $\gamma>0$. $m_{it}$ is modeled as $AR(1)$ process with the autoregressive parameter $\rho$. Varying $\gamma>0$ allows us to control for the strength of the cross-sectional correlation. Data are generated with four different structures of regressors and error terms: (a) no correlations ($\rho = 0, \gamma = 0$); (b) only serial correlation ($\rho = 0.5, \gamma = 0$); (c) only cross-sectional correlation ($\rho = 0, \gamma = 1$); and (d) both serial and cross-sectional correlations ($\rho = \{0.3,0.9\}, \gamma = 1$). In Case 2, the error terms are modeled as a spatial autoregressive (SAR(1)) process. The matrix $W=(w_{ij})_{N\times N}$ is a rook type weight matrix whose diagonal elements are zero. Note that the rows of $W$ are standardized, hence they sum to one. $\psi$ is the scalar spatial autoregressive coefficient with $|\psi|<1$. In this paper, we report the case of $\psi=0.5$. Importantly, SAR(1) model does not produce the serial correlation on the error term. In Case 3, we consider an error factor structure. Both factors and factor loadings follow AR(1) processes, which introduce both serial and cross-sectional correlations. We set $r=2$, and consider the cases of $\rho_{\lambda} = 0.3$ and $\rho_{F} = 0.9$.
In this simulation study, we examined $t$-statistics for testing the null hypothesis $H_{0}: \beta_{0}=1$ against the alternative $H_{1}:\beta_{0} \neq 1$. In each simulation we compare the proposed estimator with that of other common five types of standard errors for $\widehat{\beta}$: the standard White estimator given by $\widehat{V}_{White} = \frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\widetilde{x}_{it}\widetilde{x}_{it}'\widehat{u}_{it}^2$, where $\tilde{x}_{it}$ is demeaned version of regressor. Two types of clustered standard errors, $\widehat{V}_{CX}$ and $\widehat{V}_{CT}$, as defined in Section 2. In addition, we use two types of Newey and West HAC estimators for the panel version as follows: \\ $$\widehat{V}_{DK} = \frac{1}{NT}\sum_{t=1}^{T}\widetilde{x}_{t}'\widehat{u}_{t}\widehat{u}_{t}'\widetilde{x}_{t} + \frac{1}{NT}\sum_{h=1}^{L}\omega(h,L)\sum_{t=h+1}^{T}[\widetilde{x}_{t}'\widehat{u}_{t}\widehat{u}_{t-h}'\widetilde{x}_{t-h}+\widetilde{x}_{t-h}'\widehat{u}_{t-h}\widehat{u}_{t}'\widetilde{x}_{t}]$$ and $$\widehat{V}_{HAC} = \frac{1}{NT}\sum_{i=1}^{N}\sum_{t=1}^{T}\widetilde{x}_{it}\widetilde{x}_{it}'\widehat{u}_{it}^2 + \frac{1}{NT}\sum_{i=1}^{N}\sum_{h=1}^{L}\omega(h,L)\sum_{t=h+1}^{T}[\widetilde{x}_{it}\widehat{u}_{it}\widehat{u}_{i,t-h}\widetilde{x}_{i,t-h}'+\widetilde{x}_{i,t-h}\widehat{u}_{i,t-h}\widehat{u}_{it}\widetilde{x}_{it}']. $$ Note that $\widehat{V}_{HAC}$ assumes cross-sectional independence, while $\widehat{V}_{DK}$ allows arbitrary cross-sectional dependence. In addition, $\widehat{V}_{DK}$ and $\widehat{V}_{HAC}$ can be obtained from our proposed estimator with $M=0$ and a large constant $M$, respectively.
Results are given for sample sizes $N = 50, 200$ and $T = 100, 200$. For each $\{N, T\}$ combination, we set $L=3,7,11$ as the bandwidth for $\widehat{V}_{HAC}$, $\widehat{V}_{DK}$, and the proposed estimator, $\widehat{V}_{\text{Hard}}$. We also use Bartlett kernel for these three estimators. For the thresholding constant parameters of $\widehat{V}_{\text{Hard}}$, we set $M = 0.10, 0.15, 0.20, 0.25$ in all cases. The simulation is replicated for one thousand times for each case and the nominal significance level is 0.05. Simulation results are reported in Tables (ref) - (ref).
Tables (ref) - (ref) present the simulation results, where each table corresponds to different cases. Each table presents results of null rejection probabilities for 5% level tests based on six different standard errors. As expected, a common feature in all tables is that when both $N$ and $T$ are small, all six estimators have rejection probabilities greater than 0.05. This might happen even when the errors are drawn from i.i.d. standard normal, and this problem becomes more noticeable in the presence of serial, cross-sectional, or both correlations. A number of interesting findings based on tables are summarized below.
Tables (ref)-(ref) shows the results of Case 1. In Table (ref), Panel A indicates that all the estimators perform well due to no correlation. Especially, White standard error estimators give rejection probabilities close to 0.05. In Panel B, when the serial correlation is introduced, the performances of $\widehat{V}_{CX}$ and $\widehat{V}_{HAC}$ are markedly better than others except for small sample size. In addition, our proposed estimator, $\widehat{V}_{\text{Hard}}$, also performs well when we use both larger threshold constant $M$ and bandwidth $L$. Since there is only a serial correlation in the error term, these estimators take this correlation into account and perform well. As the size of bandwidth increases, the standard error estimated by $\widehat{V}_{HAC}$ increases to a level similar to the results of $\widehat{V}_{CX}$ and the tendency to over-reject diminishes. Since the Newey-West technique gives the weight, which is less than one, the estimated standard error may be underestimated. Hence, the traditional cluster standard error, $\widehat{V}_{CX}$, dominates the standard error of Newey-West panel version, $\widehat{V}_{HAC}$. Note that the unreported rejection probabilities of $\widehat{V}_{DK}$ exponentially increases as the bandwidth $L$ increases.
Table (ref) considers the case of cross-sectionally correlated errors and regressors. In Panel A, except the case of small sample size, $\widehat{V}_{CT}$ and $\widehat{V}_{DK}$ with small bandwidth $L$ have rejection probabilities close to 0.05 in the first panel. Also, $\widehat{V}_{\text{Hard}}$ with small $L$ and $M$ performs well. Importantly, notice that the rejection rate of $\widehat{V}_{DK}$ and $\widehat{V}_{\text{Hard}}$ tend to over-reject substantially as the lag length $L$ increases. In addition, as the cross-section size $N$ increases, the over-rejection problem becomes worse, as we mentioned in Section (ref). This tendency is easy to explain. Since $\widehat{V}_{DK}$ is an estimator based on a single time series and it is zero when full weight is given to the sample autocovariance, the bias in $\widehat{V}_{DK}$ initially falls but then increases as the lag length increases, while the variance of $\widehat{V}_{DK}$ is initially increasing but eventually becomes decreasing. Hence, $\widehat{V}_{DK}$ is biased downward substantially, and its t-statistics tends to over-reject when a large bandwidth is used. On the other hand, in the case of the small size of $L$ and $M$, $\widehat{V}_{\text{Hard}}$ gives less bias on the estimated standard error.
Panel B of Table (ref) allows the serial correlation as well as the cross-sectional correlation. Not surprisingly, all estimators except $\widehat{V}_{\text{Hard}}$ and $\widehat{V}_{DK}$ tend to over-reject substantially. In the small sample, these two estimators get worse than the case of the first panel. In the large sample, however, rejection probabilities of $\widehat{V}_{\text{Hard}}$ and $\widehat{V}_{DK}$ are close to 0.05. Importantly, $\widehat{V}_{\text{Hard}}$ outperforms $\widehat{V}_{DK}$ by choosing $M$ properly. Unreported results of $\widehat{V}_{DK}$ with larger bandwidth, $L$, show much larger rejection probabilities than that of $\widehat{V}_{\text{Hard}}$. This indicates that we can obtain unbiased standard error estimator and appropriate rejection rates using our proposed estimators, $\widehat{V}_{\text{Hard}}$. Table (ref) is the result of strong serial correlation with the cross-sectional dependence. When the serial correlation gets stronger, such as $\rho=0.9$, all estimators tend to over-reject exponentially in small samples. However, $\widehat{V}_{\text{Hard}}$ and $\widehat{V}_{DK}$ outperform other estimators as the dimensionality increases.
Table (ref) considers the error with SAR(1) structure, which does not require the serial correlation on the error term. Similar to the results reported in the first panel of Table (ref), $\widehat{V}_{CT}$ gives rejection probabilities close to 0.05. $\widehat{V}_{DK}$ and $\widehat{V}_{\text{Hard}}$ with small bandwidth $L$ also perform well. Moreover, $\widehat{V}_{\text{Hard}}$ with proper thresholding constant $M$ gives less bias than $\widehat{V}_{DK}$ on the estimated standard error.
Finally, Table (ref) presents the results of the error factor structure. Similar to the results of Table (ref), all estimators except $\widehat{V}_{\text{Hard}}$ and $\widehat{V}_{DK}$ tend to over-reject. Rejection probabilities of $\widehat{V}_{\text{Hard}}$ and $\widehat{V}_{DK}$ are relatively close to 0.05 when the sample size is large.
In this section, we re-examine the empirical work of the association between divorce law reforms and divorce rates using our proposed OLS standard error. There are many empirical studies on the effects of divorce law reforms on divorce rates. friedberg1998 found that state law reforms significantly increased divorce rates with controls for state and year fixed effects. wolfers2006did investigated the question of whether law reform continues to have an impact on the divorce rate by including dummy variables for the first two years after the reforms, 3-4 years, 5-6 years, and so on. Specifically, he studied the following fixed effect panel data model
where $y_{it}$ is the divorce rate for state $i$ and year $t$; $\alpha_{i}$ and $\mu_{t}$ are the state and year fixed effects; $X_{it,k}$ is a binary regressor that representing the treatment effect $2k$ years after the reform; $\delta_{i}t$ a linear time trend. wolfers2006did suggested that there might be two sides of the same treatment yield this phenomenon: a number of divorces gradually shifted after the earlier dissolution of bad matches, after the reform.
Both friedberg1998 and wolfers2006did estimated OLS regressions using state population weight for each year. In addition, they estimated standard errors under the assumption that errors are homoskedastic, serially and cross-sectionally uncorrelated. However, ignoring these correlations might lead to bias in the standard error estimators. We re-estimated the model of wolfers2006did using proposed OLS standard error estimators.
The same data as in wolfers2006did are used, but we exclude Indiana, New Mexico and Louisiana due to missing observations around divorce law reforms. As a result, we obtain a balanced panel data contain the divorce rates, state-level reform years and binary regressors from 1956 to 1988 over 48 states. We fit models both with and without linear time trend, and also calculate our standard errors, as well as OLS, White, cluster and HAC standard errors. We set lag choices $L=3$ for HAC and our standard errors as suggested by newey1994 ($L=4(T/100)^{2/9}$). The threshold values $M$ chosen by the cross-validation method is $M=0.2$ for the model without state-specific linear trends, and $M=0.1$ with state-specific linear trends. These $M$ values are relatively small, implying the existence of cross-sectional correlations. The estimated $\beta_{1}, \cdots, \beta_{8}$ with and without linear time trend and their different types of standard errors are presented respectively in Table (ref) below. Note that robust standard errors are not necessarily larger than the usual OLS standard errors, as shown in columns corresponding to $se_{CT}$, $se_{DK}$ and $se_{Hard}$.
In Table (ref), OLS estimates with and without linear time trend are similar to each other. These estimates are also closely comparable to the results obtained in wolfers2006did. The OLS estimates indicate that divorce rates rose soon after the law reform. However, within a decade, divorce rates had fallen over time. Most of the coefficient estimates are statistically significant at the $5\%$ level using usual OLS standard errors. According to the cluster standard errors, however, the only significant estimates are 11-15+ after the reform in the model without linear time trend. We use our method of correcting standard error estimates for heteroskedasticity, serial correlation and also cross-sectional correlation. In the model without linear trend, the estimates for 3-4 and 7-15+ are significant. On the other hand, the estimates for 1-4 are significant when linear trend is added. Our estimated standard errors are close to those of $se_{CT}$ and $se_{DK}$, which allow arbitrary cross-section correlations. The result indicates non-negligible cross-sectional correlations. The result is also consistent with kim2014, who used the interactive fixed effects approach. The latter approach is suitable for models with strong cross-sectional correlations.
This paper studies the standard error problem for the OLS estimator in linear panel models, and proposes a new standard-error estimator that is robust to heteroskedasticity, serial and cross-sectional correlations when clusters are unknown. Simulated experiments demonstrate the robustness of the new standard-error estimator to various correlation structures.