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.
21,654 characters · 6 sections · 53 citation commands
Wild Bootstrap Inference for Linear Regressions with Many Covariates
An important econometric method for estimating the causal effect of a certain variable is based on the linear regression model with the assumption that the variable of interest is exogenous after controlling for a sufficiently large set of covariates (e.g., panel data models with fixed effects or partial linear models). However, inference with many covariates can be challenging under (conditional) heteroskedasticity. cattaneo2018alternative provide a unified framework for nonstandard distributional approximation, which accommodates a variety of asymptotics studied in the literature, including the small bandwidth, many instruments, and many covariates asymptotics. Furthermore, as shown by cattaneo2018inference, under the asymptotic framework where the number of covariates, $q_n$, grows with the sample size, $n$, the consistency of the traditional heteroskedasticity-robust variance estimators would require $q_n/n \rightarrow 0$. In their seminal studies, cattaneo2018inference and jochmans2022heteroscedasticity propose new variance estimators that remain consistent even when $q_n/n \nrightarrow 0$, and provide asymptotically valid inferences that are based on normal critical values.
However, it is well known that asymptotic normal approximations may have finite sample distortions, especially when the number of observations is small. Bootstrap is found in the econometric literature to provide an improvement over asymptotic approximations in many situations; e.g., see Hall-Horowitz(1996), Horowitz(2001), Davidson-Mackinnon(2010), Djogbenou-Mackinnon-Nielsen(2019), and the references therein. However, it may fail in high-dimensional settings; e.g., see el2018can and the references therein. For linear regressions with heteroskedastic errors, mammen1993bootstrap established that the validity of nonparametric bootstrap and wild bootstrap requires $q^{1+\delta}_n/n \rightarrow 0, \delta>0$. Similarly, for linear instrumental variable (IV) models with many IVs and homoskedastic errors, Wang-Kaffo(2016) show bootstrap inconsistency for the standard residual-based bootstrap when the number of IVs is of the same order of magnitude as the sample size $n$, and propose a valid alternative bootstrap procedure. However, their procedure cannot be directly extended to the case with both many IVs and heteroskedasticity, even with a fixed number of controls.
To our knowledge, there is no existing proven-valid bootstrap method for linear regression models under both $q_n/n \nrightarrow 0$ and heteroskedasticity. Simulations in cattaneo2018inference also suggest that the nonparametric bootstrap seems to fail in their settings (e.g., see Remark 2 in their paper). In this paper, following the construction of the (approximate) cross-fit variance estimator in jochmans2022heteroscedasticity, we propose a simple modification to the standard wild bootstrap and show its asymptotic validity when $q_n/n \nrightarrow 0$. Monte Carlo experiments illustrate good finite sample performance of our procedure.
Consider the linear regression model
where $y_{i}$ is a scalar outcome variable, $x_{i}$ is a scalar explanatory variable such as a certain treatment, $w_{i}$ is a $q_n \times 1$ vector of covariates, and $u_{i}$ is an error term. We allow $q_n$ to be a non-negligible fraction of the sample size $n$. The ordinary least squares (OLS) estimator of $\beta$ is then defined as
where $\hat{v}_{i} = \sum_{j=1}^n (M_n)_{ij} x_{j}$, $(M_n)_{ij} = \{ i=j \} - w_i' (\sum_{k=1}^n w_kw_k')^{-1} w_j$, and $\{ \cdot \}$ denotes the indicator function.
Recently, cattaneo2018inference propose a new heteroskedasticity-robust variance estimator, which is consistent when $\limsup_{n} q_n/n < 1/2$. Furthermore, jochmans2022heteroscedasticity proposes an alternative variance estimator, which remains consistent as long as $\limsup_{n} q_n/n < 1$, and its corresponding $t$-ratio can be defined as $t_n = (\hat{\beta}_n - \beta_0)/\sqrt{\acute{\Omega}_n}$, where
with $\acute{u}_i = \hat{u}_i/(M_n)_{ii}$ and $\hat{u}_i = \sum_{j=1}^n (M_n)_{ij} (y_j - x_j \hat{\beta}_n)$. Then, $t$-tests based on normal critical values have correct size in large samples. However, simulations suggest that asymptotic normal approximation may not have satisfactory finite sample performance (Section (ref)), especially when the sample size is small and/or $q_n$ equals a large fraction of $n$. In Section (ref), we propose a modified wild bootstrap procedure, which follows the (approximate) cross-fit variance estimator of jochmans2022heteroscedasticity.
Our procedure is defined as follows:
Step 1: Given the null hypothesis $H_0: \beta=\beta_0$, generate the residuals $\{ \tilde{u}_i(\beta_0) \}_{i=1}^n$, where $\tilde{u}_i(\beta_0) = y_i - x_i \beta_0 - w_i' \tilde{\gamma}_n,$ and $\tilde{\gamma}_n$ is the null-restricted OLS estimator of $\gamma$.
Step 2: Generate $\{ u_i^* \}_{i=1}^n$, where $u_i^* = a_n(\beta_0) \omega^*_i \tilde{u}_i(\beta_0)$, $\{ \omega^*_i \}_{i=1}^n $ is an i.i.d. sample of random weights that are independent of the data with zero mean and unit variance while
is an adjustment factor, which takes into account the high dimensionality of the covariates. Specifically, $\acute{\Sigma}_n(\beta_0) = \sum_{i=1}^n \hat{v}_i^2 y_i \acute{u}_i(\beta_0)$, $\hat{\Sigma}_n(\beta_0) = \sum_{i=1}^n \hat{v}_i^2 \tilde{u}^2_i(\beta_0)$, $\tilde{u}_i(\beta_0)$ is defined in Step 1, and $\acute{u}_i(\beta_0) = \tilde{u}_i(\beta_0)/(M_n)_{ii}$. Notice that $\acute{\Sigma}_n(\beta_0)$ is similar to the “meat" part of $\acute{\Omega}_n$ in ((ref)), but here we impose $H_0$, as we generated the bootstrap samples with the null imposed in Step 1. Additionally, $max\{ \cdot, 1/n \}$ guards against the case that $\acute{\Sigma}_n(\beta_0)$ could take a non-positive value in finite samples.\footnote{This may also occur for the variance estimators of cattaneo2018inference and jochmans2022heteroscedasticity.}
Step 3: Generate $y^*_i = x_i \beta_0 + w_i' \tilde{\gamma}_n + u_i^*, i=1, ..., n$. Then, compute the bootstrap OLS estimator $\hat{\beta}^*_n = \left( \sum_{i=1}^n \hat{v}^2_{i} \right)^{-1} \left( \sum_{i=1}^n \hat{v}_{i} y^*_{i} \right)$, and the bootstrap residual $\hat{u}_i^* = \sum_{j=1}^n (M_n)_{ij}(y^*_j - x_j \hat{\beta}^*)$.
Step 4: Compute $t_n^* = (\hat{\beta}^*_n - \beta_0)/\sqrt{\acute{\Omega}^*_n}$, where $\acute{\Omega}^*_n = (\sum_{i=1}^n \hat{v}_i^2)^{-1} (\sum_{i=1}^n \hat{v}_i^2 (y^*_i \acute{u}^*_i)) (\sum_{i=1}^n \hat{v}_i^2)^{-1}$, with $\acute{u}_i^* = \hat{u}_i^*/(M_n)_{ii}$, i.e., $\acute{\Omega}^*_n$ has the same formula as $\acute{\Omega}_n$ but uses bootstrap samples.
Step 5: Repeat Steps 2-4 B times, and compute the bootstrap $p$-value $p_n^* = B^{-1}\sum_{b=1}^B 1\{ |t_n| > |t_n^{*(b)}| \}$. We reject $H_0$ if $p_n^*$ is smaller than the nominal level $\alpha$.
Several comments are in order.
Comment 1. In Step 1, we impose null when computing the residuals, which is advocated by Cameron(2008), Davidson-Flachaire(2008), Roodman-Nielsen-MacKinnon-Webb(2019), among others. As a consequence, we also impose null for the adjustment factor $a_n(\beta_0)$ in Step 2. When $q_n/n \rightarrow 0$, $\acute{\Sigma}_n(\beta_0)$ and $\hat{\Sigma}_n(\beta_0)$ are asymptotically equivalent so that our procedure reduces to the standard (null-imposed) wild bootstrap.
Comment 2. Instead of the above percentile-$t$ type procedure, one may consider a percentile type procedure, whose bootstrap $p$-value equals $B^{-1}\sum_{b=1}^B 1\{ | \sqrt{n}(\hat{\beta}_n - \beta_0)| > | \sqrt{n}(\hat{\beta}^{*(b)}_n - \beta_0)|\}$. However, in line with the recommendation in the literature (e.g., see the papers cited above), we find in simulations that the percentile-$t$ has better behaviours and thus also recommend using it.
Comment 3. Here we focus on the case where $x_{i}$ is a scalar variable. It is possible to extend the analysis to the case where $x_i$ is a (fixed) $d_x$-dimensional vector by considering a modification to the score bootstrap (e.g., Kline-Santos(2012)). Specifically, the procedure for testing $H_\lambda: c'\beta = \lambda$, where $c \in R^{d_x}$ and $\lambda \in R$, can be defined as follows:
(1) Obtain the null-restricted OLS estimator $\tilde{\beta}_n$ and $\tilde{\gamma}_n$, and employ it to the score contributions $\{ S_i(\tilde{\beta}_n) \}_{i=1}^n$, where $S_i(\tilde{\beta}_n) = \hat{v}_i\tilde{u}_i(\tilde{\beta}_n)$ and $\tilde{u}_i(\tilde{\beta}_n) = y_i - x'_i\tilde{\beta}_n - w_i'\tilde{\gamma}_n$;
(2) Compute the perturbed score contributions $\{ S_i^*(\tilde{\beta}_n) \}_{i=1}^n$, where $S_i^*(\tilde{\beta}_n) = \omega^*_i \hat{v}_i \tilde{u}_i(\tilde{\beta}_n)$;
(3) Compute the bootstrap statistic $T_n^* = c' \left( n^{-1} \sum_{i=1}^n \hat{v}_i\hat{v}_i' \right)^{-1} \hat{A}_n(\tilde{\beta}_n) ( n^{-1/2} \sum_{i=1}^n S_i^*(\tilde{\beta}_n))$, where $\hat{A}_n(\tilde{\beta}_n) = \acute{\Sigma}^{1/2}_n(\tilde{\beta}_n)\hat{\Sigma}^{-1/2}_n(\tilde{\beta}_n)$ is an adjustment matrix with $\acute{\Sigma}_n(\tilde{\beta}_n) = \sum_{i=1}^n \hat{v}_i \hat{v}_i' y_i \acute{u}_i(\tilde{\beta}_n)$ and $\hat{\Sigma}_n(\tilde{\beta}_n) = \sum_{i=1}^n \hat{v}_i \hat{v}_i' \tilde{u}^2_i(\tilde{\beta}_n)$;
(4) Use the distribution of $T_n^*$ conditional on the data as an estimate of the null distribution of $T_n = \sqrt{n} ( c'\hat{\beta}_n - \lambda )$.
When $x_i$ is a scalar, this modified score bootstrap is equivalent to the percentile-type procedure in Comment 2. We recommend the current procedure in Steps 1-5 as it has better finite sample performance and might also be more user-friendly to practitioners (e.g., simply need to add an adjustment factor to the residuals for the standard wild bootstrap).
Below we introduce some regularity conditions needed to establish the bootstrap validity. Assumptions (ref)-(ref) follow closely those in cattaneo2018inference and jochmans2022heteroscedasticity. Assumption (ref) contains conditions for the bootstrap procedure. Appendix (ref) gives further discussions on these assumptions.
Let $\mathcal{W}_n$ denote a collection of random variables such that $E[w_{i} | \mathcal{W}_n] = w_{i}$. Let $\epsilon_{i} = u_{i} - e_{i}$, where $e_{i} = E [u_{i} | X_n, \mathcal{W}_n]$, and $V_{i} = x_{i} - E [x_{i} | \mathcal{W}_n]$. Let $\sigma_{i}^2 = E [\epsilon_{i}^2 | X_n, \mathcal{W}_n]$, and $\tilde{V}_{i} = \sum_{j=1}^n (M_n)_{ij} V_{j}$.
Theorem (ref) shows the asymptotic validity of the modified wild bootstrap under many covariates and heteroskedasticity. For the proof of the bootstrap central limit theorem (CLT), we apply the conditional multiplier CLT (e.g., Chapter 2.9 of van1996weak).
The model is similar to that considered in cattaneo2018inference and jochmans2022heteroscedasticity:
where $x_{i} \sim i.i.d. N(0,1)$, $w_{i}$ contains a constant term and a collection of $q_n-1$ zero/one dummy variables, and $\epsilon_{i} \sim i.i.d. N(0,1)$. The dummy variables are drawn independently with success probability $\pi$ and $\gamma=0$. The sample size is fixed to $n=100$, and $q_n/n \in \{0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9\}$. The number of Monte Carlo replications is set to $10,000$ and the number of bootstrap replications is set to $199$ throughout the simulations. Following jochmans2022heteroscedasticity, we consider three designs that vary in $\beta$ and $\pi$: $\beta=1$ and $\pi=0.02$ (Design A), $\beta=1$ and $\pi=0.01$ (Design B), and $\beta=2$ and $\pi=0.02$ (Design C). The results for Design A are presented in Table (ref). The empirical null rejection frequencies of the two-sided $t$-test are given for the Eicker-White variance estimator (“HC0"), the variance estimator of cattaneo2018inference (“HCK"), the variance estimator of jochmans2022heteroscedasticity (“HCA"), and the proposed wild bootstrap with Gaussian and Rademacher random weights (“Wild-G" and “Wild-R"), respectively.
We highlight several findings below. The tests based on “HC0", “HCK", and “HCA" all tend to over-reject as $q_n$ increases, with “HCA" having the smallest distortions among the three (e.g., when $q_n/n=0.9$, their null rejection frequencies are $58.1\%$, $58.1\%$, and $17.2\%$, respectively). By contrast, both bootstraps control size across various values of $q_n$. Also, both become more conservative when $q_n$ increases, and the one with Gaussian weights seems to be slightly less conservative than that with Rademacher weights. The results for Design B and Design C are reported in Table (ref) and Table (ref), respectively. The overall patterns are very similar to those observed in Table (ref). In Section (ref) of the Supplementary Appendix, we present further simulation results for a panel data model.
In this paper, we propose a simple modification to the standard wild bootstrap procedure for linear regression models with many covariates and heteroskedastic errors. We establish its asymptotic validity when the number of covariates is of the same order of magnitude as the sample size. Our construction of the adjustment factor in the bootstrap procedure follows that of the (approximate) cross-fit variance estimator proposed by jochmans2022heteroscedasticity. Monte Carlo simulations show that the modified wild bootstrap has excellent finite sample performance compared with alternative methods that are based on standard normal critical values, especially when the sample size is small and/or the number of covariates is of the same order of magnitude as the sample size. For potential directions for future research, we note that there is a growing literature on robust inference under many (weak) instruments and non-homoskedastic errors, possibly also with many controls.\footnote{E.g., see evdokimov2018inference, crudu2021, MS22, matsushita2024jackknife, LWZ(2023), DKM24, boot-ligtenberg(2023), N23, boot2024, lim2024dimension, and yap2024inference, among others.} Due to the complexity of data structure, it is possible that asymptotic normal approximations may have less than satisfactory performance in this setting as well. On the other hand, it is found that when implemented appropriately, bootstrap approaches may substantially improve the inference accuracy for IV models, including the cases where IVs may be rather weak.\footnote{E.g., see Moreira-Porter-Suarez(2009), Davidson-Mackinnon(2008), Davidson-Mackinnon(2010), Davidson-Mackinnon(2014b), wang2015bootstrap, Wang-Kaffo(2016), Kaffo-Wang(2017), Wang-Doko(2018), Finlay-Magnusson(2019), mackinnon2023fast, and Wang-Zhang2024.} Therefore, it may be interesting to consider new bootstrap-based inference procedures that are robust to many IVs, many controls, and heteroskedastic errors simultaneously.