Gradient Wild Bootstrap for Instrumental Variable Quantile Regressions with Weak and Few Clusters

Gradient Wild Bootstrap for Instrumental Variable Quantile Regressions with Weak and Few Clusters

{18pt} {10pt} \belowdisplayskip\abovedisplayskip {5pt} \abovedisplayshortskip \belowdisplayshortskip {8pt} \belowdisplayskip\abovedisplayskip {4pt} \linespread{1.3}

Introduction

The instrumental variable (IV) regression is one of the five most widely used methods for causal inference, as highlighted by Angrist-Pischke(2008), and it is often employed in analyses involving clustered data. For instance, young2022consistency examines 1,359 IV regressions across 31 papers published by the American Economic Association (AEA), with 24 of these papers accounting for the clustering of observations. In the context of quantile regression models, where endogeneity may be present, Chernozhukov-Hansen(2004), Chernozhukov-Hansen(2005), Chernozhukov-Hansen(2006), Chernozhukov-Hansen(2008a) (hereafter referred to as CH) developed an instrumental variables quantile regression (IVQR) method. This approach offers a general IV procedure to address the endogeneity of regressors in quantile regressions, and it has been widely adopted by empirical researchers to capture the distributional effects of endogenous variables.

However, three difficulties arise when running IVQR with clustered data. First, the number of clusters is small in many empirical applications with IVs. For instance, Acemoglu2011 cluster the standard errors at the country/polity level, resulting in 12-19 clusters, Glitz2020 cluster at the sectoral level with 16 sectors, and Rogall2021 clusters at the province (district) level with 11 provinces (30 districts), respectively. ADH2013 study the effects of Chinese imports on local labor markets in the US by clustering at the state level. However, if we focus on the estimation and inference for the effects on specific regions such as the South region designated by the US Census Bureau, there are only 16 states. Furthermore, Bester-Conley-Hansen(2011) and L22 partition spatial and network data into clusters, respectively. Both papers consider the asymptotic setting in which the number of clusters is small and, thus, treated as fixed. When the number of clusters is small, conventional cluster-robust inference procedures may be unreliable for IVQR.

Second, in many applications, the strength of IVs may be relatively heterogeneous across clusters, with a few clusters providing the main identification power. For instance, Figure (ref) reports the estimated first-stage coefficients for each cluster (state) from the South region in ADH2013's (ADH2013) dataset, which suggests that there exists substantial variation in the IV strength among states. Specifically, the first-stage coefficients of some states are relatively large compared with the rest in the region. In contrast, some other states have coefficients that are rather close to zero and, thus, potentially subject to weak identification. Some states even have opposite signs for their first-stage coefficients. However, there is no existing proven-valid inference method for IVQR with few clusters, where identification may be weak in some clusters.

Third, it is also possible that IVs are weak in all clusters, in which case researchers need to use weak-identification-robust inference methods for IVQR Chernozhukov-Hansen(2008a), chernozhukov2009finite,Andrews-Stock-Sun(2019).

Motivated by these challenges, this paper investigates inference for IVQR with a small, fixed number of clusters and weak within-cluster error dependence, accommodating significant cluster-level heterogeneity in IV strength. We define clusters where the structural parameter (i.e., for the endogenous variable) is strongly identified as “strong IV clusters.” We propose a gradient wild bootstrap procedure for IVQR with clustered data. Our findings show that a bootstrap Wald test, whether studentized by the cluster-robust variance estimator (CRVE) or not, asymptotically controls size, given at least one strong IV cluster. The gradient wild bootstrap tests demonstrate power against local alternatives at the 10% and 5% significance levels with at least five and six strong IV clusters, respectively. Additionally, the CRVE-based bootstrap Wald test proves more powerful for distant local alternatives. We also develop a gradient wild bootstrap Anderson-Rubin(1949) test for IVQR that controls size regardless of instrument strength. Compared to analytical methods using HAC estimators, all our bootstrap inference methods are agnostic to within-cluster dependence structure, avoiding complex covariance matrix estimation and making it applicable to datasets with various weak dependence structures, such as cross-sectional, serial, network, and spatial dependence.

The contributions in the present paper relate to several strands of literature. First, it is related to the literature on the cluster-robust inference.\footnote{See Cameron(2008), Conley2011inference, Imbens-Kolesar(2016), abadie2022should, Hagemann(2017), Hagemann2019placebo, Hagemann2020inference, Hagemann(2019), Mackinnon-Webb(2017), Djogbenou-Mackinnon-Nielsen(2019), Mackinnon-Nielsen-Webb(2019), Ferman2019inference, Hansen-Lee(2019), Menzel2021bootstrap, Mackinnon2021, among others, and mackinnon2022cluster for a recent survey.} Hagemann(2017), Djogbenou-Mackinnon-Nielsen(2019), Mackinnon-Nielsen-Webb(2019), and Menzel2021bootstrap show bootstrap validity under the asymptotic framework in which the number of clusters diverges to infinity.\footnote{We refer interested readers to mackinnon2022cluster for detailed discussions on this asymptotic framework and the alternative asymptotic framework that treats the number of clusters as fixed.} Ibragimov-Muller(2010), Ibragimov-Muller(2016), Bester-Conley-Hansen(2011), Canay-Romano-Shaikh(2017), Hagemann2019placebo, Hagemann2020inference, Hagemann(2019), hagemann2024, and Hwang(2020) consider an alternative asymptotic framework in which the number of clusters is treated as fixed, while the number of observations in each cluster is relatively large and the within-cluster dependence is sufficiently weak. However, the inference methods proposed by BCH and Hwang(2020) require an (asymptotically) equal cluster-level sample size,\footnote{See Bester-Conley-Hansen(2011) and Hwang(2020) for details.} while those proposed by IM, CRS, and hagemann2024 would require strong identification for all clusters in the IVQR context. In contrast, our gradient bootstrap Wald tests are more flexible as they do not require an equal cluster size. In addition, they only need one strong IV cluster for size control and five to six for local power, thus allowing for substantial cluster heterogeneity in identification strength for the IVQR model. To our knowledge, no alternative method proposed in the literature remains valid in such a context. Furthermore, we provide gradient bootstrap AR tests, which are fully robust to weak identification.

Second, Canay-Santos-Shaikh(2020) and Wang-Zhang(2024) study wild bootstrap procedures with a few large clusters. In particular, Canay-Santos-Shaikh(2020) first investigates the validity of wild bootstrap by innovatively connecting it with a randomization test with sign changes. Our results for IVQR generalize and complement those in Canay-Santos-Shaikh(2020) and Wang-Zhang(2024) in the following aspects. First, Canay-Santos-Shaikh(2020) focus on the linear regression with exogenous regressors and then extend their analysis to a score bootstrap for the GMM estimator. Wang-Zhang(2024) focus on the linear IV regression and show the validity of a modified version of the cluster wild restricted efficient (WRE) bootstrap procedure (e.g., Davidson-Mackinnon(2010), Finlay-Magnusson(2019), and Roodman-Nielsen-MacKinnon-Webb(2019), among others) in the case with few clusters. Instead, we propose gradient wild bootstrap procedures for IVQR inspired by Hagemann(2017) and Jiang-Liu-Phillips-Zhang(2020), which avoid the estimation of the Hessian matrix that involves a nonparametric density component. In addition, we obtain the bootstrap estimator from the profiled optimization procedure for IVQR developed by Chernozhukov-Hansen(2004). These set us apart from the score bootstrap in the GMM setting. Second, we study the local power for our bootstrap Wald tests both with and without studentized by CRVE. In particular, the power analysis of the Wald test with CRVE is unconventional because, under a fixed number of clusters, the CRVE itself has a random limit. Specifically, we carefully design a gradient bootstrap counterpart for CRVE, which mimics well the original CRVE's randomness under the null and further diverges with the local alternative. The first property leads to the size control while the second leads to an interesting fact that the bootstrap Wald test with CRVE studentization is more powerful than that without in detecting sufficiently distant local alternatives. Such a power advantage is confirmed by our simulation experiments and empirical application. In addition, also different from its unstudentized counterpart, the local power of the Wald test with CRVE is established without the assumption that the IVQR first-stage coefficients have the same sign across all clusters, which may not hold in some empirical studies (e.g., see Figure (ref)).

Third, our paper is related to the literature on QR and IVQR. See, for example, Chernozhukov-Hansen(2004), Chernozhukov-Hansen(2005), Chernozhukov-Hansen(2006), Chernozhukov-Hansen(2008a), Hagemann(2017), and KW18. Furthermore, chernozhukov2020 provides a comprehensive overview of IVQR. We differ from them by considering an alternative asymptotic setting with a fixed number of clusters. hagemann2024 recently proposed a randomization test procedure in the spirit of CRS for inference on entire quantile and (exogenous) regression quantile processes under a small number of large clusters. However, as discussed above, a similar randomization test under the current IVQR setting would require strong identification for all clusters to establish validity. In contrast, our gradient wild bootstrap procedures for the Wald inference only need strong identification for at least one of the clusters.

Fourth, our paper is related to the literature on weak-identification-robust inference, in which various normal approximation-based inference approaches are available, among them Stock-Wright(2000), Kleibergen(2005), Andrews-Cheng(2012), Andrews(2016), Andrews-Mikusheva(2016), Moreira-Moreira(2019), and Andrews-Guggenberger(2019). However, these robust inference methods cannot be directly applied in the current context with few clusters. On the other hand, it is found in the literature that when implemented appropriately, bootstrap approaches may substantially improve the inference for linear IV models, including the cases where IVs may be rather weak,\footnote{See, for example, Davidson-Mackinnon(2008), Davidson-Mackinnon(2010), Moreira-Porter-Suarez(2009), Wang-Kaffo(2016), Kaffo-Wang(2017), Wang-Doko(2018), Finlay-Magnusson(2019), young2022consistency, and Wang-Zhang(2024), among others. In addition, tuvaandorj2021robust develops permutation versions of weak-IV-robust tests with (non-clustered) heteroskedastic errors.} although the related literature for IVQR inference remains sparse. The bootstrap AR test developed in this paper is a bootstrap counterpart of the analytical AR test proposed by Chernozhukov-Hansen(2008a). We show that it controls asymptotic size for IVQR under both weak/non-identification and a small number of large clusters.

The remainder of this paper is organized as follows. Section (ref) presents the setup, the IVQR estimation, and our gradient wild bootstrap procedures. Section (ref) presents assumptions and asymptotic results: Section (ref) gives the main assumptions, Section (ref) provides several specific examples related to our assumptions, Section (ref) provides three different methods to construct the IVs, Section (ref) presents the asymptotic results for the Wald inference, while Section (ref) presents those for the weak-identification-robust inference. Simulations in Section (ref) suggest that our procedures have outstanding finite sample size control and, in line with our theoretical analysis, the bootstrap Wald test studentized by CRVE has power advantages compared with the other bootstrap tests. The empirical application with ADH2013's dataset is presented in Section (ref).

Notation. Throughout the paper, we write $0_{d_1 \times d_2}$, $\iota_d$, and $\mathbb{I}_{d}$ as a $d_1 \times d_2$ matrix of zeros, $d$-dimensional vector of ones, and a $d \times d$ identity matrix, respectively. For any positive integer $d$, we denote $[d] = (1,\cdots,d)$. We further denote $||\cdot||_2$ and $||\cdot||_F$ as the $\ell_2$ norm for a vector and the Frobenius norm for a matrix, respectively.

Setup, Estimation, and Inference Procedure

Setup

Throughout the paper, we observe clustered data where the clusters are indexed by $j \in [J]$ and units in the $j$-th cluster are indexed by $i \in I_{n,j} = \{1, ..., n_j \}$. For the $i$-th unit in the $j$-th cluster, we observe $y_{i,j} \in \textbf{R}$, $X_{i,j} \in \textbf{R}$, and $W_{i,j} \in \textbf{R}^{d_w}$ as an outcome of interest, a scalar endogenous regressor, and exogenous regressors, respectively. Furthermore, we let $Z_{i,j} \in \textbf{R}^{d_z}$ be the exogenous variables that are excluded from the outcome equation defined through conditional CDF:

where $\Upsilon$ is a compact subset of $(0,1)$.

Throughout the paper, we focus on the setting with a single endogenous variable, as it is the most common case in empirical applications involving IVs. For instance, 101 out of 230 specifications in Andrews-Stock-Sun(2019)'s (Andrews-Stock-Sun(2019)) sample and 1,087 out of 1,359 in young2022consistency's (young2022consistency) sample feature one endogenous regressor and one IV. Similarly, lee2021 find that 61 out of 123 IV papers published in AER between 2013 and 2019 use single-IV regressions. While our setting also accommodates multiple IVs, most IVQR applications involve only one endogenous variable and one IV, as seen in studies by Chernozhukov-Hansen(2004), Chernozhukov-Hansen(2006), Chernozhukov-Hansen(2008a), and chernozhukov2013quantile. In our empirical application, we revisit the influential study by ADH2013, which also employs a single endogenous variable and one IV.

We allow the parameter of interest $\beta_n(\tau)$ (and the coefficient $\gamma_n(\tau)$ for the exogenous controls) to shift with respect to (w.r.t.) the sample size, which incorporates the analyses of size and local power in a concise manner: $\beta_n(\tau) = \beta_0(\tau) + \mu_\beta(\tau)/r_n$, and $\gamma_n(\tau) =\gamma_0(\tau) + \mu_\gamma(\tau)/r_n$, where $\mu_{\beta}(\tau) \in \textbf{R}$ and $\mu_{\gamma}(\tau) \in \textbf{R}^{d_w}$ are the local parameters and $r_n$ is the convergence rate of the score defined later. Throughout the paper, for a generic function $g$ of data $D_{i,j} = (y_{i,j},X_{i,j},W_{i,j},Z_{i,j})$, we let $\mathbb{P}_ng(D_{i,j})=\frac{1}{n}\sum_{j \in [J]} \sum_{i \in I_{n,j}} g(D_{i,j})$, $\overline{\mathbb{P}}_ng(D_{i,j})=\frac{1}{n}\sum_{j \in [J]} \sum_{i \in I_{n,j}} \mathbb{E}g(D_{i,j})$, $\mathbb{P}_{n,j}g(D_{i,j})=\frac{1}{n_j}\sum_{i \in I_{n,j}} g(D_{i,j})$, and $\overline{\mathbb{P}}_{n,j}g(D_{i,j})=\frac{1}{n_j}\sum_{i \in I_{n,j}} \mathbb{E}g(D_{i,j})$.

Estimation

Following Chernozhukov-Hansen(2006), we construct instrumental variables $\Phi_{i,j}(\tau) \in \textbf{R}^{d_\phi}$ from $(W_{i,j},Z_{i,j})$, that is, $\Phi_{i,j}(\tau) = \Phi(W_{i,j},Z_{i,j},\tau)$. For the validity of our bootstrap inference with a fixed number of clusters, we further require the instruments to be orthogonal to the control variables in the quantile regression context. Therefore, the function $\Phi(\cdot)$ may be unknown but can be estimated as $\hat{\Phi}(\cdot)$. The corresponding feasible IVs are defined as $\hat{\Phi}_{i,j}(\tau) = \hat{\Phi}(W_{i,j},Z_{i,j},\tau)$. We will provide more details about the construction of $\hat{\Phi}(\cdot)$ in Section (ref). Additionally, a scalar nonnegative weight is defined as $V_{i,j}(\tau)$, which also may be unknown, and its estimator is defined as $\hat{V}_{i,j}(\tau)$. Then, the estimation of $\beta_n(\tau)$ can be implemented via a profiled method described below. For a given value $b$ of $\beta_n(\tau)$, we first compute

where $\rho_\tau(u) = u(\tau - 1\{u \leq 0\})$. Under appropriate conditions for strong identification, which will be made clear later, we can estimate $\beta_n(\tau)$ by $\hat{\beta}(\tau)$ defined as

where $ \mathcal{B}$ is a compact subset of $\textbf{R}$, $\hat{A}_1(\tau)$ is some $d_{\phi} \times d_{\phi}$ weighting matrix, and the notation $||u||_{A}$ for a compatible vector $u$ and matrix $A$ means $(u^\top A u)^{1/2}$. Last, we define $\hat{\gamma}(\tau) = \hat{\gamma}(\hat{\beta}(\tau),\tau)$ and $\hat{\theta}(\tau) = \hat{\theta}(\hat{\beta}(\tau),\tau)$. Following the lead of Chernozhukov-Hansen(2006), in practice, we suggest setting $\hat V_{i,j}(\tau) = 1$ and $\hat A_1(\tau) = \mathbb I_{d_\phi}$ or $\hat A_1(\tau) = \mathbb P_n \hat \Phi_{i,j}(\tau)\hat \Phi_{i,j}^\top(\tau)$. When the instrument is a scalar so that the model is just identified, the choice of $\hat A_1(\tau)$ becomes irrelevant. In addition, we note that the case of just identification is always achievable because even if the original IV $Z_{i,j}$ is multi-dimensional, we can construct $\hat \Phi_{i,j}$ as the prediction of $X_{i,j}$ using $Z_{i,j}$ and $W_{i,j}$ from a first-stage linear regression, which is again a scalar (see Section (ref) for further details on the construction of $\hat \Phi_{i,j}(\tau)$).

Gradient Wild Bootstrap Inference

Inference Procedure for Wald Statistics

In this section, we consider the null and local alternative hypotheses defined as

which is equivalent to

where $\Upsilon$ is a compact subset of $(0,1)$.

Consider the test statistic with a normalization factor $\hat{A}_2(\tau)$, and let

be the test statistic. In the following, we describe the gradient wild bootstrap procedure.

Below, we discuss two cases for the normalization factor $\hat{A}_2(\tau)$: (1) it has a deterministic limit, and (2) it involves the cluster-robust variance estimator (CRVE) for score from IVQR. For case (1), given that $\hat{A}_2(\tau)$ has a deterministic limit, we do not need to bootstrap it and just let $\hat{A}_{2,g}^*(\tau) = \hat{A}_2(\tau)$ in Step 3. By an abuse of notation, the corresponding test statistic, bootstrap statistics, and critical value are still denoted as $T_n$, $T_n^{*}(g)$, and $\hat{c}_n(1-\alpha)$, respectively.

For case (2), we need some extra notation to define the normalization factor formally. Let $\hat G(\tau) \in \Re^{d_\phi}$ have a deterministic limit and

where $\omega = (0_{d_\phi \times d_w}, \mathbb{I}_{d_\phi})$. Then, we define the normalization factor $\hat{A}_{CR}(\tau)$ in case (2) as

and the corresponding CRVE-weighted Wald test statistic is defined as $$T_{CR,n} = \sup_{\tau \in \Upsilon}|| \hat{\beta}(\tau) - \beta_0(\tau)||_{\hat{A}_{CR}(\tau)}.$$

The normalization factor $\hat A_{CR}(\tau)$ and the test $T_{CR,n}$ take a cluster-robust form because the form of $\hat \Omega(\tau,\tau)$ preserves all within-cluster dependence. Naturally, for the choice of $\hat G(\tau)$, we would like to use a consistent estimator of the Jacobian $\Gamma(\tau)$ defined in (ref) below. However, this would involve a nonparametric conditional density estimation and parameter tuning. Instead, we suggest using

where $\hat {\mathcal E}_{X,\Phi}(\tau) = \mathbb P_n X_{i,j}\hat \Phi_{i,j}^\top (\tau) \hat V_{i,j}(\tau)$ and $\hat {\mathcal E}_{\Phi,\Phi}(\tau) = \mathbb P_n \hat \Phi_{i,j}(\tau) \hat \Phi_{i,j}^\top (\tau) \hat V_{i,j}(\tau)$. When $\hat A_1(\tau)$ is set as $\hat {\mathcal E}_{\Phi,\Phi}(\tau)$, we can further simplify $\hat G(\tau)$ as

In addition, even if we use a consistent estimator of the Jacobian $\Gamma(\tau)$ as $\hat G(\tau)$, it will not guarantee the consistency of the CRVE because the number of clusters is fixed in our setting. In fact, $\hat \Omega(\tau,\tau)$, and thus, $\hat A_{CR}(\tau)$ have random limits after a proper normalization. Therefore, in case (2), to construct a valid critical value for $T_{CR,n}$, we also need to bootstrap $\hat A_{CR}(\tau)$ properly to mimic well this randomness in our bootstrap samples.

To define an appropriate gradient bootstrap analogue of $\hat A_{CR}(\tau)$, we let

Then, we let $\hat{A}_{2,g}^*(\tau)$ in Step 3 of the previous bootstrap algorithm equal $\hat{A}_{CR,g}^*(\tau)$, which is defined as

where $\hat{\beta}_g^*(\tau) $ and $\hat{\gamma}_g^*(\tau) $ are defined in (ref). The bootstrap counterpart of $T_{CR,n}$ and the critical value are defined as $$T_{CR,n}^{*}(g) = \sup_{\tau \in \Upsilon}||\hat{\beta}_g^*(\tau) -\hat{\beta}(\tau)||_{\hat{A}_{CR,g}^*(\tau)} \quad \text{and} \quad \hat{c}_{CR,n}(1-\alpha), \quad \text{respectively}.$$

Inference Procedure for Weak-instrument-robust Statistics

This section considers the weak-instrument-robust inference for $\beta_n(\tau)$ when it may be weakly or partially identified. Recall $\beta_n(\tau) = \beta_0(\tau) + \mu_\beta(\tau)/r_n$. Under the null, we have $\mu_\beta(\tau) = 0$, or equivalently, $\beta_n(\tau) = \beta_0(\tau)$. Our test statistic follows the construction by Chernozhukov-Hansen(2008a). Specifically, we let

where $\hat{\theta}(b,\tau)$ is defined in (ref) and $\hat{A}_3(\tau)$ is a $d_\phi \times d_\phi$ weighting matrix, which will be specified later. We differentiate the weighting matrix used here (denoted as $\hat A_3(\tau)$) with that used for the estimation of $\beta_n(\tau)$ in (ref) (denoted as $\hat A_1(\tau)$) because of their different usages: the former is for the construction of the weak-instrument-robust test statistic while the later is for the estimation under strong identification. Theoretically, we require $\hat A_3(\tau)$ to have a deterministic limit, same as $\hat A_1(\tau)$.

Next, the bootstrap procedure for the weak-instrument-robust inference is defined as follows.

It is also possible to studentize $\hat{\theta}(\beta_0(\tau),\tau)$ by a (null-imposed) CRVE, i.e., replace $\hat{A}_3(\tau)$ with $\tilde {A}_{CR}(\tau)$, where

$\hat H(\tau) \in \Re^{d_\phi \times d_\phi}$ is some symmetric matrix, and

Note that different from $\hat{\Omega}(\tau, \tau')$ of the Wald test statistic in the previous section, we impose the null hypothesis in the form of $\tilde \Omega(\tau,\tau')$. This is essential for the validity of the test under weak/non-identification for all clusters. Additionally, note that $ \tilde {A}_{CR}(\tau)$ admits a random limit in our setting, which distinguishes it from $\hat A_3(\tau)$ used for $AR_n$.

We define the corresponding test statistic, its bootstrap counterparts, and critical value as

and $\hat{c}_{AR,CR,n}(1-\alpha)$, respectively. Naturally, for $\hat H(\tau)$, we want to use a consistent estimator for the Jacobian of $\hat \theta(\beta_0(\tau),\tau)$, which would again involve a kernel density estimation and parameter tuning. Therefore, similar to the bootstrap Wald inference, we instead suggest setting

Then, we reject the null hypothesis when $AR_{CR,n} > \hat{c}_{AR,CR,n}(1-\alpha)$. Note that different from the Wald statistic, we do not need to bootstrap $\tilde{A}_{CR}(\tau)$ here because $\tilde \Omega(\tau,\tau)$ is invariant to sign changes. Also, when the IV $\hat \Phi_{i,j}(\tau)$ is a scalar and $\Upsilon$ is a singleton, the choice of $\hat A_3(\tau)$ becomes irrelevant as it gets canceled in both the test statistic and the bootstrap critical value. Therefore, the bootstrap tests based on $AR_n$ and $AR_{CR,n}$ are numerically equivalent in this case.

Computation

Given a value of $b$, we can compute $(\hat \gamma(b,\tau),\hat \theta(b,\tau))$ in (ref) by the standard quantile regression algorithm. Then, we follow the lead of Chernozhukov-Hansen(2006) and implement a one-dimensional grid search to compute $\hat \beta(\tau)$ in (ref).

Furthermore, we note that given a value of $b$, the gradient bootstrap estimator $(\hat{\gamma}_g^*(b,\tau),\hat{\theta}_g^*(b,\tau))$ in (ref) can be formulated as linear programming and solved by well-developed linear optimization solvers. Specifically, we can stack up $y_{i,j} - X_{i,j}b$ and $\hat V_{i,j}$ first within each cluster and then across clusters $j=1,...,J$. Denote them as $\mathcal Y \in \Re^n$ and $\mathcal V \in \Re^n$, respectively. Similarly, we stack up $(W_{i,j}^\top, \hat \Phi_{i,j}^\top(\tau))$ together and construct a $n \times (d_w+d_\phi)$ matrix denoted as $\mathcal X$. Last, denote $\eta^\top = (r^\top,t^\top)$, and

By letting $u_i = \max(0,\mathcal Y_i - \mathcal X_i^\top \eta)$ and $v_i = \max(0,-\mathcal Y_i + \mathcal X_i^\top \eta)$, we can rewrite (ref) as

subject to

where $u = (u_1,\cdots,u_n)^\top$ and $v = (v_1,\cdots,v_n)^\top$. Therefore, the gradient-based wild bootstrap estimator $(\hat{\gamma}_g^*(b,\tau),\hat{\theta}_g^*(b,\tau)) = \hat \eta$, where

Assumptions and Asymptotic Results

Main Assumptions

We make the following assumptions to establish the statistical properties of our bootstrap procedures formally.

Examples

This section considers several examples and discusses why our assumptions are satisfied or violated in different scenarios.

Construction of Instruments

In Assumption (ref)(i) above, we require the IVs to satisfy the following condition: for $j \in [J]$,

where $\varepsilon_{i,j}(\tau) = y_{i,j} - X_{i,j}\beta_n(\tau) - W_{i,j}^\top \gamma_n(\tau)$ and $f_{\varepsilon_{i,j}(\tau)}(0|W_{i,j},Z_{i,j})$ is the conditional PDF of $\varepsilon_{i,j}(\tau)$ given $(W_{i,j},Z_{i,j})$ and evaluated at $0$. This section proposes three ways to construct IVs that satisfy this requirement.

Inference for Wald Statistics

Denote $Q_{\Psi,\Psi}(\tau) = \sum_{j \in [J]}\xi_jQ_{\Psi,\Psi,j}(\tau)$, $Q_{\Psi,X}(\tau) = \sum_{j \in [J]} \xi_j Q_{\Psi,X,j}(\tau)$, $Q_{\Phi,X,j}(\tau) = \omega Q_{\Psi,X,j}(\tau)$, and $Q_{\Phi,X}(\tau) = \sum_{ j \in [J]} \xi_j Q_{\Phi,X,j}(\tau)$, where $\omega = (0_{d_\phi \times d_w}, \mathbb{I}_{d_\phi})$.

We need the following assumption for the Wald inference of IVQR with CRVE.

Inference for Weak-instrument-robust Statistics

Monte Carlo Simulation

To examine the performance of the gradient wild bootstrap inference for IVQR, we consider two data-generating processes (DGPs). The first DGP is inspired by those designed by KS17 and KW18 and considers the setting with clustered data and within-cluster dependence. The second DGP is inspired by that designed by BDF09 with network dependence, and then, the clusters are obtained by spectral clustering as proposed by L22. For both designs, we set the number of observations, simulation repetitions, and bootstrap repetitions as 500, 500, and 300, respectively.

Simulation Designs

DGP 1. Let $a_{k,j} = (a_{k,1,j},\cdots,a_{k,n_j,j})^\top$ for $k=1,\cdots,d_z$, where $d_z$ is the dimension of $Z$, and $u_{j} = (u_{l,1,j},\cdots,u_{l,n_j,j})^\top$ for $l=1,2$. Then, $(a_{1,j},\cdots,a_{d_z,j},u_{1,j},u_{2,j})$ are independent, and each of them is independent across $j \in [J]$ and follows an $n_j \times 1$ multivariate normal distribution with mean zero and covariance $\Sigma(\rho_j)$, where $\Sigma(\rho_j)$ is an $n_j \times n_j$ toeplitz matrix with coefficient $\rho_j$ and $\rho_j = 0.2+0.5j/J$ for $j \in [J]$. We define the original instruments as

where $F_N(\cdot)$ is the standard normal CDF. Then, the endogenous variable $X_{i,j}$ is generated as

where $\Pi_{k,j}$ determines the identification strength. To allow for first-stage heterogeneity in the identification strength across clusters, for all $k \in [d_z]$, we let

where $\pi \in (1,1/2,1/4)$.

The outcome $Y_{i,j}$ is generated as

where $\underline{W}_{i,j}$ is distributed following $\frac{1}{2}\chi^2_{1}$ and independent of $(a_{1,j},\cdots,a_{d_z,j},u_{1j}, u_{2j})$, and $W_{i,j} = (1,\underline{W}_{i,j})^\top$. We set $\beta(u) = 1 + F_N(u)$, $\gamma = 0$, and the total number of observations $n=500$. The correlation between $X_{i,j}$ and $u_{i,j}$ is about 0.77, indicating the endogeneity level.

To allow for unbalanced clusters, we follow Djogbenou-Mackinnon-Nielsen(2019) and Mackinnon2021 and set the cluster sizes as

and $n_J = n - \sum_{j \in [j]} n_j$. We let $r=4$ when $J=9$ and $J=18$ to generate substantial heterogeneity in cluster sizes. The parameter $r$ corresponds to the heterogeneity parameter considered in the simulations of Djogbenou-Mackinnon-Nielsen(2019) and mackinnon2022cluster, and $r=4$ is the most heterogeneous setting considered in those papers. When $J=9$, the cluster sizes are $(5,8,12,19,30,48,75,117,186)$. When $J=18$, the resulting cluster sizes are $(2,2,3,4,5,7,8,10,13,17,21,26,33,41,52,65,81,110)$. We test the following hypothesis for $d_z \in \{1,3\}$ and $\tau \in \{0.1,0.25,0.5,0.75,0.9\}$:

DGP 2. We consider the linear-in-mean social interaction model detailed in Example (ref). Specifically, following L22, we generate the network $\mathcal A$ with $n=500$ nodes as $\mathcal A_{\ell,\ell'} = 1\{||\eta_{\ell} - \eta_{\ell'}||_2 \geq (7/(\pi n))^{1/2}\}$, where $\eta_\ell \stackrel{i.i.d.}{\sim} \text{Uniform}[0,1]^2$.

Then, we generate $B_\ell$ in (ref), the $\ell$-th node's background characteristic, as follows:

where $\varepsilon_\ell' \stackrel{i.i.d.}{\sim}\mathcal{N}(0,1)$. Then, the control variables $W$ and IV $Z$ can be constructed as in Example (ref).

Further denote $\tilde {\mathcal A}$ a normalized adjacency matrix with a typical entry $\tilde {\mathcal A}_{\ell,\ell'} = \mathcal A_{\ell,\ell'}/d_\ell$, where $d_\ell = \sum_{\ell' \neq \ell} {\mathcal A}_{\ell,\ell'}$ is the degree of node $\ell$. Then, we generate the outcome variable as (ref) in which $U = \{U_\ell\}_{\ell \in [n]}$ is a sequence of i.i.d. uniform (0,1) random variables and

where $\text{diag}(v)$ for a $n \times 1$ vector $v$ denotes an $n \times n$ diagonal matrix with $v$ as the diagonal, $y = (y_1,\cdots,y_n)^\top$, $B = (B_1,\cdots,B_n)^\top$, $\beta(U) = (\beta(U_1),\cdots,\beta(U_n))^\top$, and $(\delta_0(U), \delta_1(U), \delta_2(U))$ are similarly defined. In addition, we set $\delta_0(u) = 0.7683+0.25(u-0.5)$, $\beta(u) = 0.4666+0.2(u-0.5)$, $\delta_1(u) = 0.0834+0.1(u-0.5)$, and $\delta_2(u) = 0.1507+0.2(u-0.5)$. For comparison, the values $(\delta_0(0.5),\beta(0.5),\delta_1(0.5),\delta_2(0.5)) = (0.7683,0.4666,0.0834,0.1507)$ are set according to the calibration study by BDF09. Given the outcome variable, we can set the endogenous variable as $X = \tilde{\mathcal A}y$, where $X = (X_1,\cdots,X_n)^\top$. Finally, we note that (ref) holds because coefficients $(\delta_0(\cdot),\beta(\cdot))$ are monotone increasing and their corresponding regressors are positive and $\{U_\ell\}_{\ell \in [n]} \perp\!\!\!\perp \{B_\ell\}_{\ell \in [n]}$. The correlation between $X_\ell$ and $U_\ell$ is about 0.16, which indicates the endogeneity level. As discussed in Example (ref), we let $W = (\iota_n, B, \tilde A B)$ and $Z = \tilde A^2 B$.

We follow the classification procedure proposed by L22 to obtain the clusters.

We let $L$ be $(10,20)$. The number of clusters ($J$) depends on $L'$, which varies across simulation replications. Note that different from DGP 1, there is no cluster heterogeneity in the identification strength under DGP 2. We test the following hypothesis for $\tau \in \{0.1,0.25,0.5,0.75,0.9\}$:

Inference Procedure

We compare the performance of our three bootstrap inference methods $T_{CR,n}$, $T_n$, and $AR_n$ (denoted as T_CR, T, and AR in Tables (ref)-(ref)) with the conventional asymptotic Wald inference based on CRVE and two alternative methods based on cluster-level IVQR estimators available in the literature.

We construct the instrument using the nonparametric approach outlined in Remark (ref). The tuning parameters are set based on the rule-of-thumb provided in Remark (ref). Because this approach eventually produces a scalar $\hat \Phi_{i,j}(\tau)$ regardless of the dimension of the original IVs $Z_{i,j}$, our $AR_n$ and $AR_{CR,n}$ tests are numerically equivalent. Therefore, we only report the performance of $AR_n$. In addition, because we focus on inference point-wise in $\tau$ and $\hat \Phi_{i,j}(\tau)$ is a scalar, the choices of $\hat A_1(\tau)$ in the estimation of $\hat \beta(\tau)$, $\hat A_2(\tau)$ in the Wald inference based on $T_n(\tau)$, $\hat G(\tau)$ in the Wald inference based on $T_{CR,n}(\tau)$, and $\hat A_3(\tau)$ in the AR inference based on $AR_n(\tau)$ will not affect the performance of corresponding inference methods, and are set to 1 without loss of generality.

Simulation Results

Tables (ref)--(ref) and (ref) collect the simulation results for DGPs 1 and 2, respectively, when the nominal null rejection rate is set at $10\%$.\footnote{The rejection probabilities using $AR_n$ and $AR_{CR,n}$ are numerically the same with one IV.} Several key observations emerge from the results. First, T_CR, T, and AR effectively control size across various settings, while T_STD, IM, and CRS have size distortions, particularly in DGP 1, which includes both weak and strong IV clusters (i.e., there exists substantial cluster heterogeneity in the IV strength under DGP 1). The size distortions of T_STD, IM, and CRS also typically increase when $\pi$ becomes small in Tables (ref)-(ref). Second, for both DGPs, T_CR demonstrates greater power than T across various settings, aligning with our theoretical expectations. Third, we observe that T_STD typically has larger size distortion when the number of clusters is small, while the size distortions of IM and CRS tend to increase with the number of clusters. Fourth, the power of CRS and IM is comparable, and in DGP 2, their power is similar to that of T_CR, although they display slightly higher size distortion. Overall, T_CR has the best size and power performance among the six methods, and is thus our recommended inference procedure.

Empirical Applications

In an influential study, ADH2013 analyzes the effect of rising Chinese import competition on wages and employment in US local labor markets between 1990 and 2007, when the share of total US spending on Chinese goods increased substantially from 0.6% to 4.6%. In this section, we further analyze the region-wise distributional effects of such import exposure by applying IVQR and the proposed gradient wild bootstrap procedures to the Census Bureau-designated South region with 16 states and total number of observations equal to 578.\footnote{We focus on the South region because it has the highest IV strength.}

For the IVQR model, we let the outcome variable $(y_{i,j})$ denote the decadal change in the average individual log weekly wage in a given CZ. The endogenous variable $(X_{i,j})$ is the change in Chinese import exposure per worker in a CZ, which is instrumented by $(Z_{i,j})$ Chinese import growth in other high-income countries.\footnote{See Sections I.B and III.A in ADH2013 for a detailed definition of these variables.} We follow the nonparametric approach in Remark (ref) to construct $\hat \Phi_{i,j}$ used in IVQR. In addition, the exogenous variables $(W_{i,j})$ include the characteristic variables of commuting zones (CZs) and decade specified in ADH2013 as well as state fixed effects. Our IV quantile regressions are based on the CZ samples in the South region, and the samples are clustered at the state level, following ADH2013. Besides the results for the full sample, we also report those for male and female samples separately.

The main results are given in Table (ref). Specifically, for $\tau \in \{0.1, 0.25, 0.5, 0.75, 0.9\}$, we report the point IVQR estimate and the 90% bootstrap confidence sets (CSs) constructed by inverting the corresponding $AR_n$, $T_n$, and $T_{CR,n}$-based tests. The computation of the bootstrap CSs was conducted over the parameter space $[-3,1]$ with a step size of 0.01, and the number of bootstrap draws is set at 300 for each step. We can draw three key observations from the results in Table (ref). First, the impact of Chinese imports on wages shows distributional heterogeneity. Specifically, the three types of bootstrap CSs reveal that the effects are relatively significant at the high quantiles, followed by the median quantile, but not at the lower quantiles, across the full, male, and female samples. Second, the $T_{CR,n}$-based CSs are generally shorter than those of $AR_n$ and $T_n$. This aligns with our theory and simulation results, which suggest that the $T_{CR,n}$-based bootstrap test is more effective at detecting distant local alternatives. Third, the distributional effects of Chinese imports on wages are fairly consistent between males and females. In addition, we observe that the effects of Chinese imports are relatively more substantial for the male samples.