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,088 characters · 13 sections · 59 citation commands
Feasible Generalized Least Squares for Panel Data with Cross-sectional and Serial Correlations
\thispagestyle{empty}
\onehalfspacing
\setcounter{page}{1} \pagenumbering{arabic}
Heteroskedasticity, cross-sectional and serial correlations are important problems in the error terms of panel regression models. There are two approaches to deal with these problems. The first approach is to use the ordinary least squares (OLS) estimator but with a robust standard error that is robust to heteroskedasticity and correlations, for example, white1980heteroskedasticity; newey1986simple; liang1986; arellano1987; driscoll1998; hansen2007; vogelsang2012, among others. A widely used class of robust standard errors are clustered standard errors, for example, petersen2009, wooldridge2010econometric and cameron2015. bai2019olsse proposed a robust standard error with unknown clusters. In an interesting paper by abadie2017should, they argued for caution in the application of clustered standard errors since they may give rise to conservative confidence intervals. The second approach is to use the generalized least squares estimator (GLS) that directly takes into account heteroskedasticity, and cross-sectional and serial correlations in the estimation. It is well known that GLS is more efficient than OLS.
This paper focuses on the second approach. For panel models, the underlying covariance matrix involves a large number of parameters. It is important to make GLS operational. We thus consider feasible generalized least squares (FGLS). hansen2007fgls studied FGLS estimation that takes into account serial correlation and clustering problems in fixed effects panel and multilevel models. His approach requires the cluster structure to be known. This gives motivation to our paper. We assume the unknown cluster structure, and control heteroskedasticity, both serial and cross-sectional correlations by estimating the large error covariance matrix consistently. In cross-sectional setting, romano2017resurrecting obtained asymptotically valid inference of the FGLS estimator, combined with heteroskedasticity-consistent standard errors without knowledge of the conditional heteroskedasticity functional form. Moreover, miller2018feasible adapted machine learning methods (i.e., support vector regression) to take into account the misspecified form of heteroskedasticity.
In this paper, we consider (i) balanced panel data, (ii) the case of large-$N$ large-$T$, and (iii) both serial and cross-sectional correlations, but unknown stucture of clusters. We introduce a modified FGLS estimator that eliminates the cross-sectional and serial correlation bias by proposing a high-dimensional error covariance matrix estimator. In addition, our proposed method is applicable when the knowledge of clusters is not available. Following an idea suggested in bailaio2017, in this paper, the FGLS involves estimating an $NT\times NT$ dimensional inverse covariance matrix $\Omega^{-1}$, where $$ \Omega=(Eu_{t}u_{s}') $$ where each block $Eu_{t}u_{s}'$ is an $N\times N$ autocovariance matrix. Here parametric structures on the serial or cross-sectional correlations are not imposed. By assuming weak dependences, we apply nonparametric methods to estimate the covariance matrix. To control the autocorrelation in time series, we employ the idea of Newey-West truncation. This method, in the FGLS setting, is equivalent to “banding", previously proposed by bickel2008b for estimating large covariance matrices. We apply it to banding out off-diagonal $N\times N$ blocks that are far from the diagonal block. In addition, to control for the cross-sectional correlation, we assume that each of the $N\times N$ block matrices are sparse, potentially resulting from the presence of cross-sectional correlations within clusters. We then estimate them by applying the thresholding approach of bickel2008a. We apply thresholding separately to the $N\times N$ blocks, which are formed by time lags $Eu_{t}u_{t-h}': h=0,1,2,...$. This allows the cluster-membership to be potentially changing over-time. A contribution of this paper is the theoretical justification for estimating the large error covariance matrix.
For the FGLS, it is crucial for the asymptotic analysis to prove that the effect of estimating $\Omega$ is first-order negligible. In the usual low-dimensional settings that involve estimating optimal weight matrix, such as the optimal GMM estimations, it has been well known that consistency for the inverse covariance matrix estimator is sufficient for the first-order asymptotic theory, e.g., hansen1982, newey1990efficient, newey1994large. However, it turns out that when the covariance matrix is of high-dimensions, not even the optimal convergence rate for estimating $ \Omega^{-1}$ is sufficient. In fact, proving the first-order equivalence between the FGLS and the infeasible GLS (that uses the true $\Omega^{-1}$) is a very challenging problem under the large $N$, large $T$ setting. We provide a new theoretical argument to achieve this goal.
The banding and thresholding methods, which we employ in this paper, are two of the useful regularization methods. In the recent machine learning literature, these methods have been extensively exploited for estimating high-dimensional parameters. Moreover, in the econometric literature, nonparametric machine learning techniques have been verified to be powerful tools: baing2017; chernozhukov2016, chernozhukov2017; wager2018estimation, etc.
The rest of the paper is organized as follows. In Section (ref), we describe the model and the large error covariance matrix estimator. Also we introduce the implementation of FGLS estimatior and its limiting distribution. Section (ref) presents Monte Carlo studies evaluating the finite sample performance of the estimators. In Section (ref), we apply our methods to study the US divorce rate problem. Conclusions are provided in Section (ref). All proofs are given in Appendix (ref).
Throughout this paper, let $\nu_{\min}(A)$ and $\nu_{\max}(A)$ denote the minimum and maximum eigenvalues of matrix $A$ respectively. Also 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, $\ell_1$-norm and the Frobenius norm of a matrix A, respectively. Note that if $A$ is a vector, $\|A\| =\|A\|_{F}$ is equal to the Euclidean norm.
We consider a linear model \footnote{For technical simplicity we focus on a simple model where there are no fixed effects. It is straightforward to allow additive fixed effects $\alpha_i+\mu_t$ by applying the de-meaning first. The theories would be slightly more sophisticated, though such extensions are straightforward.}
The model ((ref)) can be stacked and represented in full matrix notation as
where $Y = (y_1',\cdots, y_T')'$ is the $NT \times 1$ vector of $y_{it}$ with each $y_t$ being an $N\times 1$ vector; $X = (x_1',\cdots, x_T')'$ is the $NT \times d$ matrix of $x_{it}$ with each $x_t$ being an $N\times d$; $U = (u_1',\cdots, u_T')'$ is the $NT \times 1$ vector of $u_{it}$ with each $u_t$ being an $N\times 1$ vector.
Let $\Omega = (Eu_{t}u_{s}')$ be an $NT \times NT$ matrix, consisting of many ‘‘blocks’’ matricies. The $(t,s)$th block is an $N \times N$ covariance matrix $Eu_tu_s'$. We consider the following (infeasible) GLS estimator of $\beta$:
Note that $\Omega$ is a high-dimensional conditional covariance matrix, which is very difficult to estimate. We aim to achieve the following: (i) obtain a “good" estimator of $\Omega^{-1}$, allowing an arbitrary form of weak dependence in $u_{it}$, and (ii) show that the effect of replacing $\Omega^{-1}$ by $\widehat{\Omega}^{-1}$ is asymptotically negligible.
We start with a population approximation for $\Omega$ in order to gain the intuitions. Then, we suggest the estimator for $\Omega$ that takes into account both correlations problem.
We start with a “banding" approximation to control serial correlations. Recall that $\Omega = (Eu_{t}u_{s}')$, where the $(t, s)$ block is $Eu_{t}u_{s}'$. By assuming serial stationarity and strong mixing condition, $Eu_{t}u_{s}'$ depends on $(t, s)$ only through $h=t-s$. Specifically, with slight abuse of notation, we can write $\Omega_{t,s} = \Omega_{h} = Eu_{t}u_{t-h}'$. Note for $i\neq j$, it is possible that $Eu_{it}u_{j,t-h}\neq Eu_{i,t-h}u_{jt}$, so $\Omega_{h}$ is possibly non-symmetric for $h>0$. On the other hand, $\Omega$ is symmetric due to $\Omega_{s,t} = \Omega_{t,s}'$. The diagonal blocks are the same, and all equal $\Omega_{0}=Eu_{t}u_{t}'$, while magnitudes of the elements of the off-diagonal blocks $\Omega_{h}=Eu_{t}u_{t-h}'$ decay to zero as $|h|\to\infty$ under the weak serial dependence assumption.
In the Newey-West spirit, $\Omega$ can be approximated by $\Omega^{NW} = (\Omega_{t,s}^{NW})$, where each block can be written as $\Omega_{t,s}^{NW}= \Omega_{h}^{NW}$ for $h=t-s$. Here $\Omega_{h}^{NW}$ is an $N \times N$ block matrix, defined as:
for some pre-determined $L \rightarrow \infty$. For instance, as suggested by newey1994, we can set $L$ equal to $4(T/100)^{(2/9)}$. Note that $\Omega_{h}^{NW} =\Omega_{-h}^{NW'}$. We regard $\Omega^{NW} = (\Omega_{h}^{NW})$ as the “population banding approximation”.
Next, we focus on the $N\times N$ block matrix $\Omega_{h}=Eu_{t}u_{t-h}' $ to control cross-sectional correlations. Under the intuition that $u_{it}$ is cross-sectional weakly dependent, we assume $\Omega_{h}$ is a sparse matrix, that is, $\Omega_{h,ij}=Eu_{it}u_{j,t-h}$ is “small" for “many" pairs $(i,j)$. Then $\Omega_{h}$ can be approximated by a sparse matrix $\Omega_{h}^{BL} = (\Omega_{h,ij}^{BL})_{N \times N}$ (bickel2008a), where
for some pre-determined threshold $\tau_{ij} \rightarrow 0$. We regard $\Omega_{h}^{BL}$ as the “population sparse approximation”.
In summary, we approximate $\Omega$ by an $NT\times NT$ matrix $(\widetilde \Omega^{NT}_{t,s})$, where each block $\widetilde \Omega^{NT}_{t,s}$ is an $N\times N$ matrix, defined as: for $h=t-s$, $$ \widetilde \Omega^{NT}_{t,s}:=
\quad $$ Therefore, we use “banding” to control the serial correlation, and “sparsity” to control the cross-sectional correlation. Note that an advantage of the method proposed in this paper is that it does not assume known cluster information (i.e., the number of clusters and the membership of clusters). Moreover, this method could also be modified to take into account the clustering information when available.
Given the intuition of the population approximation, we construct the large covariance estimator as follows. First, we denote the OLS estimator of $\beta$ by $\widehat{\beta}_{OLS}$ and the corresponding residuals by $\widehat{u}_{it} = y_{it} - x_{it}'\widehat{\beta}_{OLS}$.
Now we estimate the $N\times N$ block matrix $\Omega_{h}=Eu_{t}u_{t-h}'$. To do so, let $$ \widetilde{R}_{h,ij} =
,\quad and \widetilde{\sigma}_{h,ij} =
$$ where $s_{ij}(\cdot) : \mathbb{R} \rightarrow \mathbb{R}$ is a ``soft-thresholding function" with an entry dependent threshold $\tau_{ij}$ such that
where $(x)_{+} = x$ if $x\geq 0$, and zero otherwise. Here $\text{sgn}(\cdot)$ denotes the sign function, and other thresholding functions, e.g., hard thresholding, are possible. For the threshold value, we specify
for some pre-determined value $M>0$, where $\gamma_{T} = \sqrt{\frac{\log(LN)}{T}}$ is such that $\max_{h \leq L}\max_{i,j \leq N}|\widetilde{R}_{h,ij}-Eu_{it}u_{i,t-h}| = O_{P}(\gamma_{T})$. Note that the constant thresholding parameter could be allowed as bickel2008a. In practice, however, it is more desirable to have entry dependent threshold, $\tau_{ij}$. $M$ can be chosen by multifold cross-validation, which is explained in Section (ref). Then define
Next, we define the $(t, s)$th block $\widehat{\Omega}_{t,s}$ as an $N \times N$ matrix: for $h=t-s$,
Here $\omega(h,L)$ is the kernel function (see andrews1991 and newey1994). We let $\omega(h,L) = 1-h/(L+1)$ be the Bartlett kernel function, where $L$ is the bandwidth. In addition, the choice of $L$ is detailed in Section (ref). Our final estimator of $\Omega$ is an $NT\times NT$ matrix: $$ \widehat{\Omega}=(\widehat{\Omega}_{t,s}). $$ Here $\widehat{\Omega}$ is a nonparametric estimator, which does not require an assumed parametric structure on $\Omega$. Note that, for the large sample size, the proposed estimator may require a huge computational cost due to use of an $NT \times NT$ matrix.
Finally, given $\widehat{\Omega}$, we propose the feasible GLS (FGLS) estimator of $\beta$ as
Our suggested covariance matrix estimator, $\widehat{\Omega}$, requires the choice of tuning parameters $L$ and $M$, which are the bandwidth and the threshold constant respectively. We write $\widehat{\Omega}(M,L)=\widehat{\Omega}$, where the covariance estimator depends on $M$ and $L$. First, to choose the bandwidth $L$, we suggest using $L^{*} = 4(T/100)^{2/9}$, which is proposed by newey1994. For a small size of $T$, we also recommend $L \leq 3$.
The thresholding constant, $M$, can be chosen through multifold cross-validation. We randomly split the data $P$ times. We divide the data into $P=\log(T)$ blocks $J_1,...,J_P$ with block length $T/\log(T)$ and take one of the $P$ blocks as the validation set. At the $p$th split, we denote by $\widetilde{\Omega}_{0}^{p}$ the sample covariance matrix based on the validation set, defined by $\widetilde{\Omega}_{0}^{p} = |J_{p}|^{-1}\sum_{t \in J_{p}}\widehat{u}_{t}\widehat{u}_{t}'$. Let $\widetilde{\Omega}_{0}^{S,p}(M)$ be the thresholding estimator with threshold constant $M$ using the training data set $\{\widehat{u}_{t}\}_{t \notin J_{p}}$. Finally, we choose the constant $M^{*}$ by minimizing the cross-validation objective function
where $\bar C$ is a large constant such that $\widetilde{\Omega}_{0}^{S}(\bar C)$ is a diagonal matrix, and $c$ is a constant that guarantees the positive definiteness of $\widehat{\Omega}(M,L)$ for $M>c$: for each fixed $L$, $$ c= \inf[M>0: \lambda_{\min}\{\widehat{\Omega}(C,L)\}>0, \forall C>M]. $$ Here $\widetilde{\Omega}_{0}^{S}(M)$ is the soft-thresholded estimator as defined in the equation ((ref)). Then the resulting estimator of $\Omega$ is $\widehat \Omega(M^*,L^*)$.
A key step of proving the asymptotic property for $\widehat{\beta}_{FGLS}$ is to show that it is asymptotically equivalent to $\widetilde{\beta}_{GLS}^{inf}$, that is:
In the usual low-dimensional settings that involve estimating optimal weight matrix, such as the optimal GMM estimations, it has been well known that consistency for the inverse covariance matrix estimator is sufficient for the first-order asymptotic theory, e.g., hansen1982, newey1990efficient, newey1994large. It turns out, when the covariance matrix is of high-dimensions, not even the optimal convergence rate of $\|\widehat{\Omega} - \Omega\|$ is sufficient. In fact, proving equation ((ref)) is a very challenging problem. In the general case when both cross-sectional and serial correlations are present, our strategy is to use a careful expansion for $\frac{1}{\sqrt{NT}}X'(\widehat{\Omega}^{-1}-\Omega^{-1})U$. We shall proceed in two steps:\\ Step 1: Show that $\frac{1}{\sqrt{NT}}X'(\widehat{\Omega}^{-1}-\Omega^{-1})U = \frac{1}{\sqrt{NT}}W'(\widehat{\Omega}-\Omega)\varepsilon + o_{P}(1),$ where $W = \Omega^{-1}X$, and $\varepsilon = \Omega^{-1}U$.\\ Step 2: Show that $\frac{1}{\sqrt{NT}}W'(\widehat{\Omega}-\Omega)\varepsilon = o_{P}(1)$.\\
Now we suppose $\omega(h,L) =1, \Omega \approx \Omega^{NW}$ and let $A_{b_h} = \{(i,j) : |Eu_{it}u_{j,t-h}| \neq 0\}, A_{s_h} = \{(i,j) : |Eu_{it}u_{j,t-h}| = 0\}$. As for Step 2, we shall show,
Here $w_{it}$ is defined such that, we can write $W = (w_1',\cdots, w_T')'$ with $w_t$ being an $N\times d$ matrix of $w_{it}$; $\varepsilon_{it}$ is defined similarly. We then further argue that the right hand side of ((ref)) is $o_{P}(1)$ by applying a high-level Assumption (ref), which essentially saying the right hand side of ((ref)) is $o_P(1)$.
To appreciate the need of this high-level condition, let us consider a simple example as follows.
A simple example. To illustrate the key technical issue, consider a simple and ideal case where $u_{it}$ is known, and independent across both $i$ and $t$, but with cross-sectional heteroskedasticity. In this case, the covariance matrix of the $NT \times 1$ vector $U$ is a diagonal matrix, with diagonal elements $\sigma_{i}^2 = Eu_{it}^2$: $$ \Omega =
, where D =
. $$ Then a natural estimator for $\Omega$ is $$ \widehat{\Omega} =
, where \widehat{D} =
, $$ and $\widehat{\sigma}_{i}^2=\frac{1}{T}\sum_{t=1}^{T}u_{it}^2$, because $u_{it}$ is known. Then the GLS becomes:
A key step is to prove that the effect of estimating $D$ is asymptotically negligible:
It can be shown that the problem reduces to proving: $$A \equiv \frac{1}{\sqrt{NT}}\sum_{i=1}^{N}\sum_{t=1}^{T}x_{it}u_{it}\sigma_{i}^{-2}(\frac{1}{T}\sum_{s=1}^{T}(u_{is}^2-Eu_{is}^2)) {\sigma}_{i}^{-2} = o_{P}(1).$$ In fact, straightforward calculations yield
Generally, if $u_{it}|x_{it}$ is non-Gaussian and asymmetric, $E(u_{it}^3|x_{it}) \neq 0$. Hence we require $N/T \rightarrow 0$ to have $EA \rightarrow 0$. Hence, to allow for non-Gaussian and asymmetric conditional distributions, in the GLS setting it turns out $N=o(T)$ is required.
We shall not explicitly impose $N=o(T)$ in this paper as a formal assumption, but instead impose Assumption (ref). On one hand, when the distribution of $u_{it}$ is symmetric, we do not require $N=o(T)$ because as is shown in the above example, $E(u_{it}^3|x_{it})=0$ is sufficient and holds for symmetric distributions. On the other hand, when $u_{it}$ is non-symmetric, Assumption (ref) then implicitly requires $N=o(T)$. Note that $N=o(T)$ is a strong assumption in many microeconomic applications for panel data models. But as illustated in the above simple example, if $u_{it}|x_{it}$ is not symmetric, it is required for feasible GLS even if $\Omega$ is diagonal. One possible approach to weakening this assumption is to remove the higher order bias from $\widehat\Omega$. Higher order debiasing is a complicated procedure in the presence of general weak dependences. This is left for future research.
We impose the following conditions, regulating the sparsity and serial weak dependence.
Condition (ii) requires that $\Omega_{h}$ be well conditioned. Condition (iii) ensures the Bernstein-type inequality for weakly dependent data, which requires the underlying distributions to be thin-tailed. Condition (iv) is the standard $\alpha$-mixing condition, adapted to the large-$N$ panel. In addition, we impose the following regularity conditions.
Assumption (ref) allows us to prove the convergence rate of the covariance matrix estimator. Condition (ii) is an extension of the standard weak serial dependence condition to the high-dimensional case in panel data literature. It allows us to employ banding or Newey-West trunction procedure. Condition (iii) is well satisfied by various kernel functions for the HAC-type estimator. For the Bartlett kernel, for example, $$ \max_{t \leq T}\sum_{|h|\leq L}\|Eu_{t}u_{t-h}'(1-\omega(|h|,L))\| \leq \frac{1}{L}\max_{t \leq T}\sum_{|h|=0}^{\infty}\|Eu_{t}u_{t-h}'\||h| $$ converges to zero as $L\rightarrow \infty$ as long as $\max_{t \leq T}\sum_{|h|=0}^{\infty}\|Eu_{t}u_{t-h}'\||h| < \infty$.
In this paper, we assume $\Omega_{h}$ to be a sparse matrix for each $h$ and impose similar conditions as those in bickel2008a and fan2013large: write $\Omega_{h} = (\Omega_{h,ij})_{N\times N}$, where $\Omega_{h,ij}=Eu_{it}u_{j,t-h}$. For some $q \in [0,1)$, we define $$ m_{N}=\max_{|h| \leq L}\max_{i\leq N}\sum_{j=1}^{N}|\Omega_{h,ij}|^q, $$ as a measurement of the sparsity. We would require that $ m_{N}$ should be either fixed or grow slowly as $N \rightarrow \infty$. In particular, when $q=0$, $m_{N} = \max_{|h|\leq L}\max_{i \leq N}\sum_{j=1}^{N}1(\Omega_{h,ij}\neq 0)$, which corresponds to the exact sparsity case.
Let $$\gamma_{T} = \sqrt{\log(LN)/T}.$$
The following theorem shows the convergence rate of the estimated large covariance matrix. For technical simplicity, we assume that there is no fixed effects so that we do not take the de-meaning procedure. Extending to the more complete estimators with de-meaning is straightforward, but should require more technical arguments to show that the effect from added dependences due to the de-meaning is negligible.
The following conditions are required to prove Step 1 in the previous section.
Conditions (i)-(ii) require the weak cross-sectional correlations. Conditions (iii)-(iv) are the sparsity assumptions. In addition, the sparsity assumptions assume that $m_{N}$ should not be too large.
In addition, we impose the following assumption, which allows us to prove that the second term on the right hand side in the above equation is $o_{P}(1)$.
The right hand side of equation ((ref)) can be written as the equation ((ref)). Then we have the following limiting distribution by using the result of Theorem (ref).
The asymptotic variance of the FGLS estimator is $\text{Avar}(\widehat{\beta}_{FGLS})= \Gamma^{-1}/NT$, and an estimator of it is $(X'\widehat{\Omega}^{-1}X)^{-1}$. Asymptotic standard errors can be obtained in the usual fashion from the asymptotic variance estimates.
In this section we compare the proposed FGLS estimator with OLS estimator. We consider the fixed effect linear regression model, although this paper focuses on the simple linear model for technical simplicity. Hence the de-meaning procedure is applied first. The data generating process (DGP) used for the simulations is given by
where the true $\beta_{0} =1$ and fixed effects $\alpha_{i}, \mu_{t}$ are generated from $\mathcal{N}(0,0.5)$. The DGP allows for serial and cross-sectional correlation in both $x_{it}$ and $u_{it}$, which are generated by $(NT) \times (NT)$ covariance matrices, $\Omega_{X}$ and $\Omega_{U}$, as follows: let $R_{\eta} = (R_{\eta,ij})$ denote an $N \times N$ block diagonal correlation matrix. We fix the number of clusters as $G=25$. Hence, each diagonal block is a $N/G \times N/G$ matrix with the off-diagonal entries $(i,j)$ in the same cluster, $R_{\eta,ij}$ for $i \neq j$, which are generated from i.i.d. Uniform$(0,\gamma)$. In this study, we set the level of cross-sectional correlation in each cluster as $\gamma=0.3,$ or $0.7$. For the cross-sectional heteroskedasticity, let $D=\text{diag}\{d_{i}\}$, where $\{d_{i}\}_{i \leq N}$ are i.i.d. Uniform(1,$m$). Finally, we define the $N \times N$ covariance matrix of $u_{t}$ as $\Sigma_{u}=DR_{\eta}D$. In this case, we report results when $m=\sqrt{5}$. For the covariance matrix of the regressor, we simply set $\Sigma_{x} = R_{\eta}$, which does not have heteroskedasticity.
Now we introduce $i$-dependent serial correlation for the regressor and the error as follows: first let $\sigma_{ii} = \rho_{i}$ if $i=j$ and $\sigma_{ij} = \rho_{i}\rho_{j}$ if $i\neq j$. Then we define the $(NT) \times (NT)$ covariance matrix, $\Omega_{U} = (\Omega_{t,s})$. The $(t,s)$th block is an $N \times N$ covariance matrix, given by $\Omega_{t,s}= (\Omega_{t,s}(i,j))$, where $\Omega_{t,s}(i,j) = \Sigma_{u,ij}\sigma_{ij}^{|t-s|}$. The large covariance matrix of the regressor, $\Omega_{X}$, is generated similarly. The level of $i$-dependent $\rho_{i}$ of the regressor and the error is generated from i.i.d. Uniform$(0,0.6)$, seperately.
Note that the $(t,s)$th block covariance decays exponentially as $|t-s|$ increases. Finally we generate the $NT \times 1$ vectors $(u_{1}',...,u_{T}')' = \Omega_{U}^{1/2}\zeta$, where $\zeta$ is an $NT \times 1$ vector, whose entries are generated from i.i.d. $\mathcal{N}(0,5)$. Similarly, the regressor is generated by $(x_{1}',...,x_{T}')' = \Omega_{X}^{1/2}\xi$, where $\xi$ is an $NT \times 1$ vector, whose entries are generated from i.i.d. $\mathcal{N}(0,1)$. Note that $x_{it}$ is uncorrelated with $u_{it}$.
In this numerical study, we use sample sizes $N=50, 100$ and $T=50, 100, 150$, and the simulation is replicated for one thousand times in all cases.\footnote{The procedure of proposed estimators require use of an $NT \times NT$ matrix as discussed in Section (ref). Indeed, when $NT$ is large, the procedure appears to be computationally demanding. Hence, we focus on the small sample size in this study.} For each $\{N, T\}$ combination, we set the bandwidth $L = 3$ in all cases. The threshold constant, $M$, is obtained by the cross-validation method as suggested in Section (ref). For instance, when $T=100$, the number of folds to split is $\log(100) \approx 5$. In general, the cross-validation chooses $M$ between 1.4 and 1.8. Interestingly, as the level of cross-section correlation increases, the cross-validation tends to choose smaller M, so that the number of non-thresholded elements increases. Hence it takes into account the strength of cross-sectional correlation. We use the Bartlett kernel for our FGLS estimator. Results are summarized in Tables (ref)-(ref).
Tables (ref)-(ref) present the simulation results, where each table corresponds to a different level of cross-sectional correlation, $\gamma = \{0.3, 0.7\}$. In each table, the mean and standard deviation of the estimators are reported. FGLS(Diag) refers to the FGLS estimator using the diagonal covariance matrix, which only takes into account heteroskedasticity. RMSE is the ratio of the mean squared error of FGLS to that of OLS. The mean and standard deviation of the estimated standard errors for OLS and FGLS are also reported. The robust unknown clustered standard error, suggested by bai2019olsse, is used for OLS. For FGLS, we report the results of the standard error as introduced in Theorem (ref). The difference between the standard deviation of the estimators and the mean of standard errors can be explained as the bias of estimated standard errors. In addition, we present null rejection probabilities for the 5% level tests using the traditional $\mathcal{N}(0,1)$ critical value based on each standard errors.
According to Tables (ref)-(ref), we see that both methods are almost unbiased, while our proposed FGLS has indeed smaller standard deviation of $\widehat{\beta}$ than that of OLS and FGLS(Diag). In all cases, the RMSE of our proposed FGLS is significantly smaller than one. Hence the results confirm that the FGLS estimator is more efficient than the OLS and the FGLS(Diag) estimators in presence of heteroskedasticity, serial and cross-sectional correlations. Regarding the $t$-test, in Table (ref), the rejection probabilities of FGLS and OLS are close to 0.05 when $T$ is large, while those of FGLS(Diag) tend to over-reject. Since the FGLS(Diag) estimator does not take into account the serial and the cross-sectional correlations, its standard errors are underestimated. On the other hand, in Table (ref), we find that the standard errors of all estimators are underestimated and the t-test rejection probabilities are much larger than 0.05, especially when $T$ is relatively smaller than $N$ (e.g., $N=100$ and $T=50$). This is due to the strong cross-sectional correlation within clusters. However, the rejection probabilities of FGLS and OLS are much smaller than those of FGLS(Diag). In summary, FGLS does improve efficiency in terms of mean squared error; also we obtain unbiased standard error estimator and appropriate rejection rate as $T$ increases.
In the literature, the cause of the sharp increase in the U.S. divorce rate in the 1960-1970s is an important research question. During 1970s, more than half of states in the U.S. liberalized the divorce system, and the effects of reforms on divorce rates have been investigated by many such as allen1992marriage and peters1986marriage. With controls for state and year fixed effects, friedberg1998 suggested that state law reforms siginificantly increased divorce rates. Also, she assumed that unilateral divorce laws affected divorce rates permanently. However, divorce rates from 1975 have been subsequently decreasing according to empirical evidence. Therefore the question of whether law reforms also affect the divorce rate decrease has arisen. wolfers2006did revisited this question by using a treatment effect panel data model, and identified only temporal effects of reforms on divorce rates. In particular, he used dummy variables for the first two years after the reforms, 3-4 years, 5-6 years, and so on. More specifically, the following fixed effect panel data model was considered:
where $y_{it}$ is the divorce rate for state $i$ and year $t$, $\alpha_{i}$ a state fixed effect, $\mu_{t}$ a time fixed effects, and $\delta_{i}t$ a linear time trend with unknown coefficient $\delta_{i}$. $X_{it}$ is a binary regressor which denotes the treatment effect $2k$ years after the reform. wolfers2006did suggested that “the divorce rate rose sharply following the adoption of unilateral divorce laws, but this rise was reversed within about a decade". He also concluded that “15 years after reform the divorce rate is lower as a result of the adoption of unilateral divorce, although it is hard to draw any strong conclusions about long-run effects’’.
Both friedberg1998 and wolfers2006did used a weighted model by muliplying all variables by the square root of state population. In addition, they used ordinary OLS standard error, which does not take into account heteroskedasticity, serial and cross-sectional correlations. However, standard errors might be biased when one disregards these correlations. Therefore, we re-estimated the model of wolfers2006did using the proposed FGLS method and OLS with the heteroskedastic standard errors of white1980heteroskedasticity, the clustered standard error of arellano1987, and the robust standard error of bai2019olsse.
The same dataset as in wolfers2006did is used, which includes the divorce rate, state-level reform years, binary regressors, and state population. Due to missing observations around divorce law reforms, we exclude Indiana, New Mexico and Louisiana. As a result, we obtain balanced panel data from 1956 to 1988 for 48 states. We fit the models both with and without linear time trend, and use OLS and FGLS in each model to estimate $\beta$. In the FGLS estimation, we set bandwidth $L=3$ as proposed by newey1994 ($L=4(T/100)^{2/9}$). The thresholding values are chosen by the cross-validation method as discussed in Section (ref), more specifically, $M=1.8$ and $M=1.9$ for the model with and without linear time trends, respectively. The Bartlett kernel is used in the OLS robust standard error and FGLS estimation. The estimated $\beta_{1}, \cdots, \beta_{8}$ with and without linear time trend and standard errors are summarized in Table (ref) below.
The OLS and FGLS estimates in both models are similar to each other. The results show that divorce rates rose soon after the law reform. However, within a decade, divorce rates had fallen over time. Interestingly, FGLS confirms the negative effects of the law reforms on the divorce rates, specifically, 11-15+ years after the reform in the model with state-specific linear time trends, and 9-15+ years after the reform in the model without state-specific linear time trends. In addition, the FGLS estimates for 1-6 and 1-4 years are positive and statistically significant in the models with and without linear time trends, respectively. For OLS, the coefficient estimates for 3-4 and 7-15+ are significant in the model without linear time trends based on $se_{BCL}$. In contrast, the OLS estimates are statistically significant only for 1-4 years when a linear time trend is added. According to the clustered standard error, $se_{CX}$, note that only 11-15+ are statistically significant in the model without trends.
According to OLS and FGLS estimation results with and without a linear time trend, we make the following conlcusion: in the first 8 years, the overall trend of divorce rate is increasing, but the law reform reduces the divorce rate after 3-4 years. However, 8 years after the reform, we observe that the law reform has a negative effect on divorce rate. Note that wolfers2006did de-emphasized the negative coefficient at the end of the periods, as these are not robust to inclusion of state-specific quadratic trends, which we did not employ in this paper. Overall, the results of FGLS estimates are consistent with wolfers2006did.
In this paper, we propose a large covariance matrix estimator and a modified version of FGLS that takes into account both serial and cross-sectional correlations in linear panel models that are robust to heteroskedasticity, serial and cross-sectional correlations. The covariance matrix estimator is asymptotically unbiased with an improved convergence rate. It is shown to be more efficient than other existing methods in panel data literature. From simulated experiments, we confirmed that our FGLS estimates are more efficient than OLS estimates.