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.
93,270 characters · 11 sections · 37 citation commands
Cluster-Robust Standard Errors for Linear Regression Models with Many Controls
}
\let\markeverypar\everypar \newtoks\everypar \everypar\markeverypar \markeverypar{\the\everypar\looseness=-2\relax} \thispagestyle{empty} \parskip=3.pt \baselineskip=20pt
\setcounter{page}{1}
It is common practice in empirical work to use standard errors and associated confidence intervals that are robust to heteroskedasticity and/or various forms of dependence.
In particular, since \citet*{moulton} highlighted the importance of accounting for dependence arising in data with a group structure, researchers often assume that the data are clustered at some economically relevant level, e.g. by individual unit or geographical location.
The justification for this type of inference procedures is asymptotic, in the sense that their validity relies on the assumption that the sample size is large relative to the number of parameters in the model. In small samples two issues arise: (i) confidence intervals based on the usual Gaussian approximation become invalid, (ii) robust standard errors are biased. A variety of methods that address these issues in the context of the linear regression model have been proposed in the literature. However, these usually alleviate but do not entirely solve the problem and are not always appealing given their ad hoc nature (see imbenskolesar, for a discussion).
Furthermore, while modern datasets usually include a large number of observations, the assumption that the number of parameters in the estimated model is negligible relative to the sample size can still be unattractive, even when the researcher's goal is to conduct inference on a small set of parameters. For example, in many important applications of the linear regression model, the object of interest is $\bm{\beta}$ in a model of the form
where $y_{i,n}$ is a scalar outcome variable, $\mathbf{x}_{i,n}$ is a $d\times1$ vector of regressors of fixed dimension, $\mathbf{w}_{i,n}$ is a vector of covariates of possibly “large" dimension $K_n$, and $u_{i,n}$ is an unobserved scalar error term. In many applications of this model, the assumption that $K_n/n\to0$ is unpalatable or even violated, as researchers often include a large set of covariates in $\mathbf{w}_{i,n}$ in order to control for observed and unobserved confounders (see discussion below).
Motivated by the above observations, this paper develops inference theory for linear regression models with many controls and clustering. In particular, we first show that the usual cluster-robust standard errors by \citet*{liangzeger} are inconsistent in general when $K_n/n\nrightarrow 0$. We then propose a new clustered standard error formula that allows to carry out valid inference on $\bm{\beta}$ under asymptotics in which $K_n$ is allowed (but not required) to grow as fast as the sample size.
The findings of this paper contribute to the long-established literature initiated by \citet*{whitebook} dealing with cluster-robust inference in a variety of models, a recent review of which is given by \citet*{cameronmiller}; see, e.g., \citet*{arellano1987}, \citet*{bellmccaffrey}, \citet*{hansen2007}, cameron2008, \citet*{ibragimovmuller}, \citet*{pustejovskytilton} and canaycluster2008. In particular, our analysis is related to a literature, reviewed in \citet*{imbenskolesar}, in which bias-reduction modifications of standard errors and particular distributional approximations are proposed with the aim of improving the performance of cluster-robust inference procedures in small samples. We contribute to this literature by studying a new general class of cluster-robust variance estimators that allows to fully correct such “small-sample" bias, while also exploiting the particular structure of the model in ((ref)) to circumvent the need for non-Gaussian distributional approximations.
This paper also adds to a sizeable body of literature that deals with inference procedures in models that involve the estimation of many incidental parameters; see, e.g., \citet*{angristhahn2004}, \citet*{hahnewey}, \citet*{stockwatson2008}, belloni2014rev, \citeauthor*{cjn} (cjn, cattaneo2018partially), \citet*{verdier}, and references therein. In particular, our findings can be seen as generalising those of \citet*{cjn}, CJN hereafter, who establish asymptotic normality of the OLS estimator of $\bm{\beta}$ in ((ref)) when $K_n/n\nrightarrow0$, and provide inference methods under such asymptotics when the errors are independent and heteroskedastic.
The results in this paper were derived independently of \citet*{lithesis}, who tackles the problem of cluster-robust variance estimation for the full set of coefficients of a generic high-dimensional linear model and obtains a similar estimator to ours. However, his results are not directly applicable to inference and are silent about the extent to which sufficient conditions for consistent variance estimation restrict the underlying data generating process of the regressors.
The rest of this paper is organized as follows. Section 2 introduces the framework of the paper and illustrates its relevance using three leading examples. Section 3 discusses our assumptions. Section 4 presents our main theoretical results. Section 5 reports the findings of a Monte Carlo study. Section 6 presents an empirical illustration. Section 7 briefly concludes. Proofs and extensions of the results are given in the Appendix.
The main object of interest in our analysis is $\bm{\beta}$ in ((ref)), on which we would like to carry out inference while treating the high-dimensional $\mathbf{w}_{i,n}$ as nuisance covariates. A natural choice of estimator for $\bm{\beta}$ is the OLS estimator, which can be written as
where $M_{ij,n}=\mathbbm{1}\{i=j\}-\mathbf{w}_{i,n}'(\sum_{k=1}^{n}\mathbf{w}_{k,n}\mathbf{w}_{k,n}')^{-1}\mathbf{w}_{j,n}$ is the $(i,j)$ entry of the symmetric and idempotent annihilator matrix $\mathbf{M}_n$, with $\mathbbm{1}\{\cdot\}$ denoting the indicator function. Defining $\bm{\hat{\Gamma}}_n=\sum_{i=1}^n\mathbf{\hat{v}}_{i,n}\mathbf{\hat{v}}_{i,n}'/n$ and $\bm{\Sigma}_n$ the (conditional) variance of $\sum_{i=1}^n\mathbf{\hat{v}}_{i,n}u_{i,n}/\sqrt{n}$, it is well-known that, when $n\to\infty$ and $K_n$ is fixed, the asymptotic distribution of $\bm{\hat{\beta}}_n$ is
When the errors are assumed to be correlated only within $G_n$ clusters of bounded size, $\bm{\Sigma}_n$ can be estimated consistently with the popular cluster-robust variance estimator by Liang and Zeger (1986, LZ hereafter):
where $\mathcal{T}_{g,n}$ denotes the subset of observations contained in cluster $g$ and $\{\mathcal{T}_{g,n}:1\leq g\leq G_n\}$ is a partition of the data. As a result, asymptotically valid inference can be carried out using the usual testing procedures based on the distributional approximation $\bm{\hat{\beta}}_n\stackrel{a}{\sim}\mathcal{N}(\bm{\beta},\bm{\hat{\Gamma}}_n^{-1}\bm{\hat{\Sigma}}_n^{\texttt{LZ}}\bm{\hat{\Gamma}}_n^{-1}/n)$.
The objective of this paper is to establish cluster-robust inference procedures for $\bm{\beta}$ under asymptotics in which $K_n/n\nrightarrow 0$. Allowing the dimension of the nuisance covariates $K_n$ to grow at the same rate as the sample size $n$ enables us to cover many relevant applications of the general model in ((ref)).
To simplify exposition, we present our inference theory for linear regression models with many controls for the case of strictly exogenous regressors. While all the results of this paper can be well-understood for this special case, their generalisation to (potential) misspecification bias in the model is straightforward and is provided in the Appendix.
In this section we present a set of assumptions for the special case of strict exogeneity of the regressors. A more general set of assumptions that allows for misspecification bias is given in the Appendix.
Suppose that $\{(y_{i,n},\mathbf{x}_{i,n}',\mathbf{w}_{i,n}') : 1 \leq i \leq n\}$ is generated by ((ref)) and set $\mathcal{X}_{n}=(\mathbf{x}_{1,n},\dots,\mathbf{x}_{n,n})$ and $\mathcal{W}_{n}=(\mathbf{w}_{1,n},\dots,\mathbf{w}_{n,n})$. We define the following quantities:
where $\mathbf{v}_{i,n}=\mathbf{x}_{i,n}-(\sum_{j=1}^{n}\mathbb{E}[\mathbf{x}_{j,n}\mathbf{w}_{j,n}'])(\sum_{j=1}^{n}\mathbb{E}[\mathbf{w}_{j,n}\mathbf{w}_{j,n}'])^{-1}\mathbf{w}_{i,n}$ is the population counterpart of $\mathbf{\hat{v}}_{i,n}$. Also, letting $\lambda_{\text{min}}(\cdot)$ denote the minimum eigenvalue of its argument, define
where $\mathbf{V}_{i,n}=\mathbf{x}_{i,n}-\mathbb{E}[\mathbf{x}_{i,n}|\mathcal{W}_n]$, $\mathbf{\tilde{\Gamma}}_n=\sum_{i=1}^{n}\mathbf{\tilde{V}}_{i,n}\mathbf{\tilde{V}}_{i,n}'/n$ and $\mathbf{\tilde{V}}_{i,n}=\sum_{j=1}^{n}M_{ij,n}\mathbf{V}_{i,n}$.\\ We impose the following three assumptions:
Assumption 1 defines the sampling structure, in which we allow for arbitrary dependence within clusters of finite but possibly heterogenous size for both the regressors and the errors. In terms of clustering structure, the resulting asymptotics are the same as the usual ones of \citet*{whitebook} and liangzeger in which $n, G_n\to\infty$ and $G_n\propto n$. We expect that the results of this paper would generalize to asymptotics where cluster sizes are allowed to diverge with $n$ and $G_n$, as considered in \citet*{hansen2007} and \citet*{hansenlee}. It is likely that such extension would require imposing more restrictive conditions on the regression design and the distributional properties of the errors, e.g. stationarity and/or mixing, and we leave it to future work.
Assumption 2 allows for asymptotics where $K_n/n\nrightarrow 0$, while imposing standard restrictions on the regression design and some bounds on the (conditional) higher-order moments of the structural residuals $u_{i,n}$ and $\mathbf{V}_{i,n}$.
The condition on $\chi_n$ in Assumption 3 is a requirement on the quality of the linear approximation for the conditional expectation $\mathbb{E}[\mathbf{x}_{i,n}|\mathcal{W}_n]$. The high-level condition $\max_{1\leq i \leq n}\|\mathbf{\hat{v}}_{i,n}\|/\sqrt{n}=o_p(1)$ also places restrictions on the relationship between $\mathbf{x}_{i,n}$ and $\mathbf{w}_{i,n}$ and has a central importance in establishing asymptotic normality of the OLS estimator for $\bm{\beta}$ and consistency of our proposed variance estimator. cjn show that this restriction holds under mild moment conditions when either (i) $K_n/n\to0$, or (ii) $\chi_n=o(1)$ or (iii) $\max_{1\leq i \leq n}\sum_{i=1}^n\mathbbm{1}\{M_{ij,n}\neq0\}=o_p(n^{1/3})$. While condition (i) is not the case of primary interest of this paper, (ii) and (iii) accomodate $K_n/n\nrightarrow 0$ and can be used to verify the high-level condition $\max_{1\leq i \leq n}\|\mathbf{\hat{v}}_{i,n}\|/\sqrt{n}=o_p(1)$ when $\mathbf{w}_{i,n}$ can be interpreted as approximating functions, dummy/discrete variables or fixed effects. See cjn for details.
This section presents our main theoretical results for inference in linear regression models with many controls and clustering under the set of simplified assumptions presented in Section 3. Proofs of the theorems and other auxiliary results are given in the Appendix for the general case that allows for misspecification bias in the model.
Our first result extends the asymptotic normality result for $\hat{\bm{\beta}}$ previously derived by cjn to the case of clustering.
Theorem 1 implies that the asymptotic distribution of $\hat{\bm{\beta}}$ under $K_n/n\nrightarrow0$ resembles the standard one obtainable under fixed-$K_n$.\footnote{From Assumptions 1-3 it also follows that $\bm{\hat{\Omega}}_n=O_p(1)$, implying that $\bm{\hat{\beta}}_n$ is $\sqrt{n}$-consistent.} As a result, confidence intervals can be constructed using the usual Gaussian approximation and the problem of conducting valid inference reduces to finding a consistent estimator for $\mathbf{\Sigma}_n$ under our asymptotics of interest.
For our discussion of variance estimation, we introduce a new class of estimators. Let $\bm{\Omega}_{u,n}=\mathbb{E}[\mathbf{u}_n\mathbf{u}_n'|\mathcal{X}_n,\mathcal{W}_n]$ be the (conditional) variance-covariance matrix of the errors $\mathbf{u}_{n}=(u_{1,n}, \dots, u_{n,n})'$ and $L_n=\sum_{g=1}^{G_n}(\#\mathcal{T}_{g,n})^2$ the number of non-zero elements contained in it. We define a general class of cluster-robust estimators for $\bm{\Sigma}_n$ of the form
where $\kappa_{g_1,g_2,i_1,j_1,i_2.j_2,n}$ is an entry of the $L_n\times L_n$ symmetric matrix $\bm{\kappa}_{n}=\bm{\kappa}_{n}(\mathbf{w}_{1,n}, \dots, \mathbf{w}_{1,n})$.\footnote{In particular, $\kappa_{g_1,g_2,i_1,j_1,i_2.j_2,n}$ corresponds to the $(h(g_1,i_1,j_1), h(g_2,i_2,j_2))$ entry of $\bm{\kappa}_{n}$, where $h(g,i,j)=[\sum_{k=0}^{(g-1)}(\#\mathcal{T}_{k,n})^2+(\#\mathcal{T}_{g,n})(i-1) +j]$ and we adopt the convention that $\#\mathcal{T}_{0,n}=0$.} Notice that by setting $\bm{\kappa}_n=\mathbf{I}_{L_n}$ one obtains the usual cluster-robust estimator by liangzeger:
The next theorem provides an asymptotic representation for this class of estimators.
Heuristically, in Theorem 2 consistency of $\bm{\hat{\beta}}_n$ implies that the estimated residuals $\hat{u}_{i,n}$ asymptotically converge to $\tilde{u}_{i,n}=\sum_{j=1}^{n}M_{ij,n}u_{j,n}$, which are only affected by the estimation noise due to projecting out the high-dimensional covariates $\mathbf{w}_{i,n}$.
The result of Theorem 2 has a central importance in our analysis. First, it immediately provides an explicit characterization for the asymptotic limit of LZ's estimator, as shown in the following corollary.
Corollary 1 implies that inference based on LZ's clustered standard errors is invalid in general under asymptotics where $K_n/n\nrightarrow0$. In fact, $\bm{\hat{\Sigma}}_n^{\textup{\texttt{LZ}}}$ does not converge to the target $\bm{\Sigma}_n$ due to elements of $\mathbf{M}_n$ arising in its asymptotic limit. While the sign of the asymptotic bias of LZ's estimator cannot be determined in general, $\bm{\hat{\Sigma}}_n^{\textup{\texttt{LZ}}}$ will typically underestimate $\bm{\Sigma}_n$.\footnote{A particular case in which $\operatorname*{plim}{\bm{\hat{\Sigma}}_n^{\textup{\texttt{LZ}}}}\leq \bm{\Sigma}_n$ holds in general is when the true residuals are in fact homoskedastic, which can be shown using arguments from Theorem 1 in \citet*{bellmccaffrey}.} Intuitively, the “asymptotic" regression residuals $\tilde{u}_{i,n}$ tend to be smaller than the true residuals as a result of the overfitting due to the high-dimensional controls. In addition, estimated residuals will tend to have lower intra-cluster correlation than the true errors \citep*{bellmccaffrey}. Inference based on $\bm{\hat{\Sigma}}_n^{\textup{\texttt{LZ}}}$ is therefore expected to be asymptotically liberal in most applications.
Furthermore, Theorem 2 suggests that a particular choice of $\bm{\kappa}_n$ might set the leading term in the expansion ((ref)) equal to the target $\bm{\Sigma}_n$. Based on this insight, we define the estimator
where $\bm{\kappa}^{\texttt{CR}}_n$ solves the system of $L_n(L_n-1)/2$ equations
It also turns out that $\bm{\kappa}^{\texttt{CR}}_n$ can be characterized in closed form as
where $\otimes$ denotes the Kronecker product and $\mathbf{S}_n$ is the $n^2 \times L_n$ selection matrix with full column rank such that $\mathbf{S}_n'\text{vec}(\bm{\Omega}_{u,n})$ is the $L_n\times 1$ vector containing the non-zero elements of $\bm{\Omega}_{u,n}$.
In the next theorem we establish consistency of our proposed estimator.
Since $\mathbf{S}_n'(\mathbf{M}_n\otimes\mathbf{M}_n) \mathbf{S}_n$ is observable, the first high-level condition in Theorem 3 is expected to be verified whenever $\mathbf{S}_n'(\mathbf{M}_n\otimes\mathbf{M}_n) \mathbf{S}_n$ is invertible. The second high-level condition could be verified using Theorem 1 of \citet*{varah}, which provides a bound for $\left\lVert\bm{\kappa}_{n}^{\textup{\texttt{CR}}}\right\rVert_\infty$ under the condition that $\mathbf{S}_n'(\mathbf{M}_n\otimes\mathbf{M}_n) \mathbf{S}_n$ is diagonally dominant. In simulations we find that diagonal dominance typically does not hold but our high-level condition is verified in a wide range of models and designs, as shown in Section 5.\footnote{cjn instead develop their theory under the requirement that $\mathbf{M}_n\odot\mathbf{M}_n$ is diagonally dominant. It would be interesting to investigate whether this requirement could be relaxed in practice.}
The structure of our proposed estimator is related to the cluster-robust variance estimator proposed by \citet*{bellmccaffrey}, which corresponds to a particular choice of block-diagonal $\bm{\kappa}_n$ that sets the the bias of the variance estimator to 0 only in the special case in which $\bm{\Omega}_{u,n}=\sigma^2\mathbf{I}_n$, i.e. the true residuals are in fact homoskedatic. Differently from \citet*{bellmccaffrey}, our choice of correction matrix $\bm{\kappa}^{\texttt{CR}}_n$ induces an averaging over cross-products of estimated residuals not just within but also across clusters, thus allowing to set the leading term in expansion ((ref)) equal to $\bm{\Sigma}_n$ in general.
The results of this paper can be easily extended to a more general version of the variance estimators, described in Section (ref) of the Appendix, that allows to impose within-cluster zero restrictions on the variace-covariance matrix of the errors. In such form, our proposed estimator reduces to the one of \citet*{stockwatson2008} in the case of one-way fixed effects panel data models with zero restrictions on the conditional autocovariances of $U_{it}$ within entities. While our results cover a much wider class of models, they also partly improve on \citet*{stockwatson2008} as we do not require $(\mathbf{X}_{i1}',\dots, \mathbf{X}_{iT}', U_{i1}, \dots, U_{iT})$ to be i.i.d. nor we require $(\mathbf{X}_{it}, U_{it})$ to be stationary.
Although consistency of $\bm{\hat{\Sigma}}_n^{\textup{\texttt{CR}}}$ is derived under asymptotic sequences that allow but do not require $K_n/n\nrightarrow0$, it is still desirable to establish consistency of LZ's estimator under some sufficiently slow rate of growth for $K_n$. For this purpose, define $\mathbf{w}^{*}_{i,n}=\mathbf{w}_{i,n}\bm{\hat{\Sigma}}_{\mathbf{w},n}^{-1/2}$, where $\bm{\hat{\Sigma}}_{\mathbf{w},n}^{1/2}$ is the unique symmetric positive definite $K_n\times K_n$ matrix such that $\bm{\hat{\Sigma}}_{\mathbf{w},n}^{1/2}\bm{\hat{\Sigma}}_{\mathbf{w},n}^{1/2}=\frac{1}{n}\sum_{i=1}^{n}\mathbf{w}_{i,n}\mathbf{w}_{i,n}'$. The following theorem provides sufficient conditions for consistency of LZ's cluster-robust estimator.
Although we can only prove consistency of LZ's estimator under $K_n^2/n\rightarrow0$, we speculate that $K_n/n\to0$ might suffice in general. We leave the refinement of this result for future work.\footnote{Theorem 4 also states that $K_n/n\to0$ is sufficient for consistency of LZ's estimator in the special case of homoskedastic errors. }
This section reports the findings of a simulation study that investigates the finite sample behaviour of the cluster-robust variance estimators studied in this paper. We consider three distinct designs motivated by the empirical examples covered by the theoretical framework of this paper: the linear regression models with increasing dimension, the semiparametric partially linear model and the fixed effects panel data regression model.
The chosen designs for our Monte Carlo experiments closely resemble those of cjn, also borrowing from specifications in \citet*{stockwatson2008} and \citet*{mackinnon}. The data generating process (DGP) for the linear regression model with many covariates is:
where $\mathbf{w}_{gi}\stackrel{i.i.d.}{\sim}\mathcal{U}(-1,1)$, $\bm{\iota}=(1,1, \dots,1)'$, $\beta=1$, $\bm{\gamma=0}$, $\rho=0.3 $, the constants $\varkappa_{x}$ and $\varkappa_{u1}$ are chosen so that $\mathbb{V}[x_{gi}]=\mathbb{V}[U_{g1}]=1$ and $t(a)=a\mathbbm{1}(-2\leq a \leq 2)+2\text{sgn}(a)(1-\mathbbm{1}(-2\leq a \leq 2))$. Table (ref) reports the results of our experiment for five dimensions of $\mathbf{w}_{gi}$: $K \in \{1, 71, 141, 211, 281\}$, where the first covariate is an intercept, as well as three different numbers of equal-sized clusters: $G\in \{175, 70, 35 \}$. We consider three different estimators for the variance of the OLS estimator $\hat{\beta}$: the unfeasible estimator based on $\bm{\hat{\Sigma}}^{\texttt{Unf}}_n=\frac{1}{n}\sum_{g=1}^{G_n}\sum_{i,j\in\mathcal{T}_{g,n}}\mathbf{\hat{v}}_{i,n}\mathbf{\hat{v}}'_{j,n}U_{i,n}U_{j,n}$ that makes use of the true error realizations, the classical estimator by LZ and our proposed cluster-robust formula, as previously defined. For each of these estimators, we report the bias (expressed in percentage), the standard deviation (denoted by Std.) and the empirical coverage probability (denoted by $\hat{p};\alpha$) of the Gaussian confidence interval of the form:
where $\Phi^{-1}$ denotes the inverse of the standard normal cumulative distribution function $\Phi$, $\hat{\Sigma}_{\ell}$ with $\ell\in\{\texttt{Unf, LZ, CR} \}$ corresponds to the variance estimators already discussed and we set $\alpha=0.05$.
The findings from this experiment are in line with our theoretical predictions. Firstly, we find that inference based on LZ's clustered standard errors formula is highly inaccurate. In fact, its bias quickly increases with the dimensionality of the model, resulting in substantial undercoverage even for $K/n=0.101$. On the other hand, our proposed estimator performs well, with negligible bias and close-to-correct empirical coverage even for $K/n=0.401$. Such improvement in inference accuracy compared to LZ's estimator is achieved in spite of a decrease in relative precision. As expected, the performance of all estimators is adversely affected by a reduction in the number of clusters. In Tables (ref)-(ref) we also report on the behaviour of $\left\lVert\bm{\kappa}_{n}^{\texttt{CR}}\right\rVert_\infty$ in this design; we find that $\left\lVert\bm{\kappa}^{\texttt{CR}}_{n}\right\rVert_\infty$ does not only seem to be bounded but even decreasing as $n$ grows.\footnote{Notice that diagonal dominance of $\mathbf{S}_n'(\mathbf{M}_n\otimes\mathbf{M}_n) \mathbf{S}_n$ does not hold in any of the simulations carried out in this section.}
Analogous results are found for a different version of this experiment that considers independent and discrete controls constructed as $\mathbbm{1}\{\mathcal{N}(0,1)\geq 1\}$, as reported in Tables (ref) and (ref)-(ref).
The experimental design chosen for the semiparametric partially linear model takes the form:
where $\text{dim}(\mathbf{z}_{gi})=6$, $\mathbf{z}_{gi}=(z_{1,gi},\dots,z_{6,gi})'$ with $z_{\ell,gi}\stackrel{i.i.d.}{\sim}\mathcal{U}(-1,1)$, $\ell=1,\dots,6$. The unknown regressions functions are set to $g(\mathbf{z}_{gi})=\text{exp}\big(-\left\lVert\mathbf{z}_{gi}\right\rVert^{1/2}\big)$ and $h(\mathbf{z}_{gi})=\text{exp}\big(\left\lVert\mathbf{z}_{gi}\right\rVert^{1/2}\big)$, and the constants $\varkappa_{v}$ and $\varkappa_{u1}$ are again chosen so that $\mathbb{V}[x_{gi}]=\mathbb{V}[u_{g1}]=1$. Similarly to the previous simulation, we set $\beta=1$ and $\rho=0.3$. To construct the covariates $\mathbf{w}_{gi}$ entering the estimated linear regression model $y_{gi}=\bm{\beta}'\mathbf{x}_{gi} + \bm{\gamma}_n'\mathbf{w}_{gi} + u_{gi}$, we consider power series expansions. The table below gives a summary of the expansions considered, where $\mathbf{w}_{gi}=\mathbf{p}(\mathbf{z}_{gi};K)$ for $K \in \{1, 7, 13, 28, 34, 84, 90, 210, 216\}$ is defined as follows:
The results for this experiment are given in Table (ref), in which we only report $K\in \{1, 13, 34, 90, 216\}$ for reasons of parsimony. The numerical findings are largely consistent with those reported for the other two simulation models. Although $\left\lVert\bm{\kappa}_{n}^{\texttt{CR}}\right\rVert_\infty$ has bigger magnitude in this setting compared to the other simulation models, it still appears to be bounded (see Tables (ref)-(ref)).
The main difference between this setting and the linear model with increasing dimension considered previously is that the unfeasible estimator that uses realizations of the true structural disturbances is free not just from estimation error but also specification error, which in turn affects LZ's and our proposed estimator when $K$ is small; in addition, the degree of heteroskedasticity and dependence in the errors is invariant with respect to the dimensionality of the model, since it only depends on $x_{gi}$ and $\mathbf{z}_{gi}$ but not $\mathbf{w}_{gi}$.
For fixed effects panel data regression model we consider the following specification:
where $\alpha_i$ is a time-invariant individual effect and $e_{d_{it}}$ are unobserved factors common to all observations sharing the same value of the indexing variable $d_{it}\in\{1,\dots, N_d\}$. This model coincides with the one studied in Verdier (2018), whose theory and simulation results concern the case of two-way clustering. We instead consider the case of one-way clustering at the individual level as we postulate the following DGP:
where $\text{dim}(\mathbf{z}_{it})=6$, $\mathbf{z}_{it}=(z_{1,it},\dots,z_{6,it})'$ with $z_{\ell,gi}\stackrel{i.i.d.}{\sim}\text{Uniform}(-1,1)$, $\ell=1,\dots,6$, the constants $\varkappa_{x}$ and $\varkappa_{u1}$ are chosen so that $\mathbb{V}[x_{it}]=\mathbb{V}[U_{i1}]=1$, the function $t(\cdot)$ is as previously defined and we set $\beta=1$ and $\alpha_i=e_{d_{it}}=0$. For the purpose of estimation, we transform ((ref)) by partialling out the individual fixed effects $\alpha_i$, so that the estimated model $\tilde{y}_{it}=\beta \tilde{x}_{it} + \tilde{e}_{d_{it}}+ \tilde{u}_{it}$ has $\text{dim}(\mathbf{w}_{i})=N_d$.\footnote{The motivation for this transformation is that $(\mathbf{S}_n'(\mathbf{M}_n\otimes\mathbf{M}_n) \mathbf{S}_n)$ is not invertible when the controls $\mathbf{w}_{i,n}$ include indicators for the clusters (see, e.g., \citealp*{stockwatson2008}). Notice that partialling out the fixed effects does not affect the correlation structure of the errors.} We consider $G=N=[700/T]$ for $T\in\{4, 10, 20\}$, as well as $N_d=700/r$ for $r\in\{700, 10, 5, 4, 3\}$, so that the total sample size is always roughly $n=700$. Tables (ref) and (ref)-(ref) report the numerical findings of this experiment, which are consistent with our theoretical predictions and in line with the results obtained for the other simulation models.
In this section we illustrate the use of the inference methods discussed in this paper by revisiting \citet*{dl2001} study of the impact of abortion on crime rates.
\citeauthor*{dl2001} (dl2001, henceforth DL) put forward the hypothesis that the legalization of abortion in the United States in the 1970s played a major role in explaining the sharp decline in crime observed two decades later. In particular, they describe two causal channels through which abortion might affect crime. The first is that abortion reduces the absolute size of a cohort, resulting in lower crime 15-25 years later, when its members are at the highest risk of engaging in criminal activities. The second channel is ascribed to the increased control over fertility that abortion provides to women. In fact, women may use abortion to optimize the timing of childbearing, thus ensuring that the child grows in a more favourable environment, e.g. when a father is present in the family, the mother is better educated and household income is stable. As a result, increased access to abortion is expected to cause a reduction in crime levels even if fertility rates were to remain constant.
In order to estimate the impact of abortion on crime, dl2001 consider state-level yearly data for the period 1985-1997 and propose a model for crime rates whose basic specification is
where $i$ indexes states, $t$ indexes the time period, $c\in\{\text{violent, property, murder}\}$ indexes the type of crime, $y_{cit}$ is the crime-rate for crime type $c$; $a_{cit}$ is measure of abortion rate relevent for crime type $c$; $z_{it}$ is a set of time-varying state-specific controls consisting of the log of lagged prisoners per capita, the log of lagged police per capita, the unemployment rate, per-capita income, the poverty rate, AFDC generosity at time $t-15$, a dummy for concealed weapons law and beer consumptions; $\theta_{ci}$ are state fixed effects; and $\lambda_{ct}$ are time fixed effects. Further details on data definitions and the institutional background can be found in the original paper.
The results from estimating the baseline model in ((ref)) are reported in Table (ref) and resemble those in dl2001, although not identical as we have excluded Washington DC from the sample.\footnote{We exclude Washington DC for simplicity, as it produces similar results to dl2001 and circumvents the need to introduce the estimation weights used in their paper.} Following dl2001, we report standard errors clustered at the state level. These estimates indicate a strong (and statistically significant) negative association between abortion and crime, as they imply that an increase in the abortion rate of 100 per 1,000 live births is associated with a reduction in crime rates between 9 and 13 per cent, depending on the type of crime. However, the extent to which this association can be interpreted as causal crucially depends on the assumption that abortion rates can be taken as random after controlling for a national trend, time-invariant state-specific confounders and $z_{it}$. Even if one believes that abortion rates can be taken as exogenous conditional on the the controls included by dl2001, one can still expect the assumption that they enter the structural equation for crime rates linearly as in ((ref)) to be too restrictive. For example, \citet*{foote} have argued that the results in dl2001 might not be robust to the inclusion of state-specific trends.\footnote{In a response to \citet*{foote}, \citet*{dl2008} reexamine their orginal study and use a longer panel to argue that their original results are robust to the inclusion of state-specific linear time trends.}
For these reasons, we consider a model for crime rates and abortion in which the controls $\mathbf{z}_{it}$ are allowed to enter in a much more flexible way compared to dl2001. In particular, we consider a version of the high-dimensional regression model studied in this paper where in addition to the controls included by dl2001, we include first-order interactions, quadratics, cumulative values and interactions of those variables and their initial values with a quadratic trend; in addition, we also include the interaction between the initial level of abortion and a quadratic trend. Once we stack all these regressors and the time-effects $\lambda_{ct}$ in the vector $\mathbf{w}_{it}$ and absorb the state-effects, we obtain a regression model of the same form as ((ref)):
where the dimension of the high-dimensional controls is $K_n=105$, resulting in $K_n/n\approx0.161$. Estimates for the causal effect of abortion on crime based on ((ref)) are given in Table (ref), where we also report standard errors based on the variance estimators considered in this paper. These estimates are qualitatively similar to those obtained for the baseline model considered in dl2001, and interestingly imply an even more sizeable negative effect of abortion on crime rates for all types of crime. The statistical significance of these effects however crucially depends on the choice of standard errors. In fact, clustered standard errors based on the variance estimator proposed in this paper are between 42 and 74 per cent bigger than the traditional clustered stardard errors by liangzeger, depending on the type of crime. In the case of violent crime, for example, the estimated coefficient for abortion rates has associated p-value below 1 per cent when traditional clustered standard errors are used, while the use of our proposed standard errors leads to failure to reject the hypothesis of no effect of abortion on crime at the 5 per cent level.
This empirical illustration showcases the relevance of the inference methods proposed in this paper. In this particular application, the inclusion of many controls arises naturally as a way to flexibly control for observable state-level characteristics and trends that are allowed to depend on those characteristics. Our approach in this particular application resembles the one adopted by belloni2014rev, who also re-examine the empirical setting in dl2001 to illustrate the use of their proposed inference method for treatment effects with many controls based on LASSO double-selection. They consider a similar specification of the high-dimensional model in ((ref)) but allow for an even more flexible specification that includes higher-order interactions of the variables we consider (and a few additional ones, such as initial differences of $\mathbf{z}_{it}$) with cubic trends, which gives $K_n/n\approx0.500$ in their application.\footnote{Interestingly, their estimates imply statistically non-significant impact of abortion on all types of crime.} While belloni2014rev's (belloni2014rev) method is naturally suited to handle such large number of controls, its validity relies on the assumption that the effect of confounding factors can be controlled for by a small number of variables (“approximate sparsity").
Our proposed inference procedure therefore offers a valuable alternative to selection-based methods in settings where the inclusion of a relatively large number of controls is expected to yield a reasonable approximation of the structural relationship of interest, while circumventing the need to impose requirements of sparsity on the model.
This paper provides inference results for the OLS estimator of a subset of coefficients in linear regression models with many controls and clustering. We show that the usual cluster-robust variance estimator by \citet*{liangzeger} does not deliver consistent standard errors when the number of controls is a non-vanishing fraction of the sample size, typically resulting in confidence intervals with coverage below the nominal size. We then propose a new clustered standard error formula that is robust to the inclusion of many controls. Monte Carlo evidence supports our theoretical results and shows that our proposed variance estimator performs well in finite samples.
While our results are presented for the case of one-way clustering, we expect that they can be easily adapted to the generalisation of our methods to multi-way clustering proposed by \citet*{verdier}. It would also be of interest to investigate whether the analysis of this paper could be extended to cases where variance estimation does not rely on zero restrictions on the covariance matrix of the errors, e.g. when time series or spacial dependence in the errors is assumed.
\setstretch{1.2} \setcounter{section}{0} \setcounter{equation}{0} \setcounter{theorem}{0} \setcounter{coro}{0} \setcounter{assum}{0}