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.
52,508 characters · 10 sections · 32 citation commands
Two-way Clustering Robust Variance Estimator in Quantile Regression Models
\thispagestyle{empty}
Quantile regression (QR), introduced by koenker1978regression, is a widely used tool for studying heterogeneous effects and tail risks in economics and finance. In many empirical environments, however, observations are indexed by multiple clustering dimensions and exhibit dependence along each of them. A canonical example is a two-way array $\{(y_{gh},X_{gh}):g=1,\ldots,G,\;h=1,\ldots,H\}$ in which observations can be correlated within the $g$-dimension and within the $h$-dimension because of latent shocks shared by units in the same row or column (e.g., worker $\times$ firm, exporter $\times$ destination).
This paper develops a unified large-sample theory and feasible inference procedures for linear QR under two-way clustering. We study the conditional quantile regression model and allow for rich two-way dependence using an Aldous--Hoover--Kallenberg (AHK) representation for separately exchangeable arrays aldous1981representations,hoover1979relations,kallenberg1989representation. This framework has become a standard device for modeling multi-way clustered dependence and for deriving projection-based asymptotics for array data davezies2021empirical,menzel2021bootstrap,chiang2023standard,graham2024sparse. Building on this structure, we establish a self-normalized central limit theorem that accommodates regime-dependent rates and delivers asymptotic normality.
We then propose a feasible two-way cluster-robust variance estimator (CRVE) for QR of the familiar sandwich form \[ \widehat{\Sigma}(\tau)=\widehat{D}(\tau)^{-1}\,\widehat{\Omega}(\tau)\,\widehat{D}(\tau)^{-1}. \] The “bread” $\widehat{D}(\tau)$ is a kernel-based estimator of the conditional density at the target quantile, adapted here to two-way clustering, while the “meat” $\widehat{\Omega}(\tau)$ aggregates row- and column-cluster covariance contributions along with a residual component in a manner that mirrors the underlying projection decomposition.
Four features fundamentally complicate establishing the consistency of $\widehat{\Sigma}(\tau)$ relative to standard two-way clustered mean regression (e.g., cameron2011robust,mackinnon2021wild) and to one-way clustered quantile regression (e.g., parente2016quantile,hagemann2017cluster). First, unlike mean regression, the quantile score is non-smooth, which makes uniform control of score fluctuations in neighborhoods of $\beta_{0}(\tau)$ more delicate. Second, the Jacobian depends on the conditional density at zero and is estimated nonparametrically, so the proof must control the bias and stochastic error of a kernel-based “bread” under clustering. Third, unlike one-way clustered quantile regression, the effective convergence rate of $\widehat{\beta}(\tau)$, denoted $r_{GH}$, can vary across dependence regimes: depending on the relative magnitudes of the row, column, and interaction components of the score, different projection terms may dominate the leading stochastic fluctuation. Fourth, two-way dependence precludes reducing the sample score to a sum of independent (or weakly dependent) terms along either dimension without explicitly isolating the row, column, and interaction components. More importantly, these four difficulties are not additive. In our setting, non-smooth scores and kernel Jacobian estimation must be handled simultaneously with regime-dependent rates and genuinely two-way dependence, requiring a uniform analysis of both the score and the Jacobian that remains valid across dependence regimes.
Monte Carlo results confirm that conventional QR standard errors can severely understate uncertainty when two-way clustering is present, whereas the proposed CRVE delivers reliable coverage across a wide range of dependence configurations, including one-way clustering and cluster-independent settings.
Furthermore, we characterize the boundary of uniform inference. When the interaction component remains asymptotically non-negligible while the row and column components are weak, the limiting distribution of $\widehat{\beta}(\tau)$ may be non-Gaussian. In this regime, the distribution of the normalized estimator depends sensitively on the underlying DGP. We show that over a natural class of two-way clustered triangular arrays, no procedure can uniformly consistently approximate the asymptotic distribution of $\widehat{\beta}(\tau)$. Uniform inference over the full model class is therefore unattainable without additional structure.
In an empirical application, we study how teacher-licensing stringency relates to the supply of high-quality teachers. Consistent with prior evidence, we find little indication that stricter licensing affects high-quality candidates on average. However, this average pattern masks substantial heterogeneity: a negative effect is concentrated in the lower part of the distribution of high-quality outcomes, suggesting that some relatively strong candidates are close to the margin between teaching and other careers and are therefore sensitive to increases in licensing costs. In contrast, for the higher quantiles we find little evidence that increased stringency discourages right-tail teacher.
Recently and independently, chiang2024extremal study extremal quantiles under two-way clustered dependence, focusing on rare-event estimation in two-way clustered data. While menzel2021bootstrap show that sample means may exhibit non-Gaussian limits under two-way clustering, chiang2024extremal demonstrate that extremal quantiles can remain robust even in degenerate dependence regimes. Our paper complements this line of work by focusing on interior quantiles: while chiang2024extremal analyze $\widehat{\beta}(\tau)$ as $\tau\to0$, we consider fixed $\tau\in(0,1)$, where the non-smooth quantile score and the interaction of row and column components have fundamentally different implications for the limiting distribution and the validity of inference.
More broadly, our results contribute to the literature on robust quantile regression inference under dependence, including kernel-based theory under weak dependence kato2012asymptotic, CRVE and pigeonhole bootstrap results for GMM (e.g. quantile IV) under multiway clustering davezies2018asymptotic, and recent advances in weak-dependence-robust covariance estimation for quantile regression galvao2024hac.
Relative to these papers, our contributions are threefold. First, we establish asymptotic normality of Powell's kernel estimator under two-way clustering, and derive a feasible optimal bandwidth rule. Second, to the best of our knowledge, we provide the first feasible two-way cluster-robust variance estimator for quantile regression at a fixed $\tau\in(0,1)$, and prove its uniform validity whenever the Gaussian limit arises. Although davezies2018asymptotic propose multiway variance estimation for GMM, their approach is not directly applicable here because it relies on a plug-in Jacobian that requires knowledge of the true conditional density, which is typically unavailable in practice. Moreover, unlike the setting emphasized in davezies2018asymptotic, we do not impose nondegeneracy of the asymptotic variance: the rate of convergence is allowed to vary with the strength of clustering, and our variance estimator is designed to adapt across these regimes. Third, when the limiting distribution in two-way clustered quantile regression is non-Gaussian, we show that uniform consistency of inference is impossible.
The remainder of the paper is organized as follows. Section (ref) introduces the two-way QR framework and the projection decomposition and develops the regime-adaptive limit theory for $\widehat{\beta}(\tau)$. Section (ref) proposes $\widehat{D}(\tau)$ and $\widehat{\Omega}(\tau)$ and establishes the consistency of $\widehat{\Sigma}(\tau)$ and the validity of inference based on the $t$-statistic. Section (ref) presents Monte Carlo evidence. Section (ref) studies an application to teacher licensing. Section (ref) concludes. Technical proofs and additional results are deferred to the appendix.
Let $\{(y_{ghi},X_{ghi}^{\top}):g=1,\dots,G,h=1,\dots,H,i=1,\dots,N_{gh}\}$ be an array of observations, where $y_{ghi}\in\mathbb{R}$ is the scalar response and $X_{ghi}\in\mathbb{R}^{d}$ is a vector of regressors. The first index $g$ identifies the cluster in the first dimension (the $g$-cluster), and the second index $h$ identifies the cluster in the second dimension (the $h$-cluster). The index pair $(g,h)$ therefore labels a cell formed by the intersection of a $g$-cluster and an $h$-cluster (e.g., unit $\times$ time). Let $N_{gh}\in\mathbb{N}$ denote the number of observations within cell $(g,h)$.
Fix a quantile index $\tau\in(0,1)$. We consider the quantile regression model
where $Q_{y_{ghi}}(\tau\vert X_{ghi})$ denotes the conditional $\tau$-quantile of $y_{ghi}$ given $X_{ghi}$. The quantile error $e_{ghi}(\tau)$ is defined as $e_{ghi}(\tau):=y_{ghi}-X_{ghi}\beta_{0}(\tau)$. Let $\rho_{\tau}(u):=u\bigl(\tau-\mathbf{1}\{u\le0\}\bigr)$ denote the check loss. For two-way clustered data, the QR estimator minimizes \[ \hat{\beta}(\tau):=\arg\min_{\beta\in\Theta}\frac{1}{GH}\sum_{g=1}^{G}\sum_{h=1}^{H}\sum_{i=1}^{N_{gh}}\rho_{\tau}\!\bigl(y_{ghi}-X_{ghi}^{\top}\beta\bigr), \] with respect to $\beta\in\Theta\subset\mathbb{R}^{d}$, where $\Theta$ is compact.
For later use, define the quantile score
The function $\Psi_{ghi}(\tau)$ is nonlinear due to the indicator function, which plays a central role in the asymptotic analysis. For each cell $(g,h)$, let $X_{gh}$ be the $N_{gh}\times d$ matrix with $i^{\text{th}}$ row $X_{ghi}$, and let $y_{gh}$ and $e_{gh}(\tau)$ be the corresponding $N_{gh}\times 1$ vectors with $i^{\text{th}}$ elements $y_{ghi}$ and $e_{ghi}(\tau)$. We impose the conditional quantile restriction $Q_{e_{gh}}(\tau\vert X_{gh})=0$, i.e., the conditional $\tau$-quantile of $e_{gh}(\tau)$ given $X_{gh}$ equals zero. For simplicity, we focus on the case where each cell contains exactly one observation, that is, $N_{gh}=1$ for all $g,h$, and suppress the replicate index $i$. Extensions to heterogeneous $N_{gh}$ are provided in the Internet Appendix.
In the two-way clustering literature, the Aldous--Hoover--Kallenberg (AHK) representation is widely used; see, for example, davezies2021empirical, mackinnon2021wild, and chiang2023standard.
Under Assumption (ref), the array $(y_{gh},X_{gh})$ is separately exchangeable across $(g,h)$, and hence marginally identically distributed, though generally dependent. The quantile index $\tau$ affects the model only through the conditional quantile restriction and does not enter the regressor process. There exists a measurable function $\Psi(U,V,W;\tau)$ such that $\Psi_{gh}(\tau)=\Psi(U_{g},V_{h},W_{gh};\tau).$ The score then admits the Hoeffding type decomposition
where
This decomposition follows from $L^{2}$ projection theory for separately exchangeable arrays and is unique in $L^{2}$. Closely related decompositions for nonlinear statistics under AHK dependence have been developed recently for U-statistics on bipartite and row--column exchangeable arrays; see le2025hoeffding. When convenient, we write $\Psi_{\bullet}^{(j)}$ for $\Psi^{(j)}(\cdot;\tau)$, $j=\text{I},\ldots,\text{IV}$. We suppress the dependence on $\tau$ to conserve space.
By construction, $E[\Psi^{(j)}]=0$ for each $j$ and $E[\Psi_{\bullet}^{(j)}\Psi_{\bullet}^{(j')\top}]=0$ for $j\neq j'$. Although $(U_{g},V_{h},W_{gh})$ are independent, the components $\Psi_{g}^{(\text{I})}$, $\Psi_{h}^{(\text{II})},$ $\Psi_{gh}^{(\text{III})}$, and $\Psi_{gh}^{(\text{IV})}$ need not be. These components are, however, pairwise orthogonal in $L^{2}$, which suffices to characterize asymptotic variances and limit distributions.
Let $f_{e\vert X}(e\vert x)$ denote the conditional density of $e_{gh}$ given $X_{gh}=x$. $f_{e\vert X}^{(1)}(e\vert x)$ and $f_{e\vert X}^{(2)}(e\vert x)$ denote the corresponding first and second derivatives, respectively. We impose a natural two-way array analogue of the standard moment, smoothness, and nonsingularity conditions used in i.i.d. quantile regression.
For $j\in\{\text{I},\text{II},\text{III},\text{IV}\}$, define the component variances \[ \sigma_{j,\Gamma}^{2}:=E\!\left[\Psi^{(j)}\Psi^{(j)\top}\right]. \] The subscript $\Gamma$ emphasizes that these quantities depend on the underlying DGP and may vary with $(G,H)$. To simplify notation, we suppress the explicit $(G,H)$ dependence. For each $j$, we further assume, for expositional convenience, that the diagonal elements of $\sigma_{j,\Gamma}^{2}$ are of the same order; we use the first diagonal entry, $\sigma_{j,1\Gamma}^{2}$, to represent this order.
A standard argument yields the Bahadur representation \[ \hat{\beta}-\beta_{0}(\tau)\;+\;o_{P}\!\left(\bigl\|\hat{\beta}-\beta_{0}(\tau)\bigr\|\right)\;=\;D(\tau)^{-1}\,\bar{\Psi}_{GH}, \] where \[ D(\tau):=E\!\left[f_{e\vert X}(0\vert X_{gh})\,X_{gh}X_{gh}^{\top}\right],\qquad\bar{\Psi}_{GH}:=\frac{1}{GH}\sum_{g=1}^{G}\sum_{h=1}^{H}\Psi_{gh}. \] Using (ref), we decompose $\bar{\Psi}_{GH}$ as
Conditional on the latent variables, the arrays $\{\Psi_g^{(\mathrm{I})}\}_{g=1}^G$ and $\{\Psi_h^{(\mathrm{II})}\}_{h=1}^H$ are i.i.d. across clusters, and, conditional on $\{U_g,V_h\}$, $\{\Psi_{gh}^{(\mathrm{IV})}\}_{g,h}$ are i.i.d. across cells. Consequently, after appropriate normalization, the sums associated with $\bar{\Psi}^{(\text{I})}$, $\bar{\Psi}^{(\text{II})}$, and $\bar{\Psi}^{(\text{IV})}$ are asymptotically Gaussian, whereas $\bar{\Psi}^{(\text{III})}$ may admit a non-Gaussian limit.
Let the asymptotic variance of $\hat{\beta}$ be \[ \Sigma_{GH}:=D(\tau)^{-1}\,\Omega_{GH}(\tau)\,D(\tau)^{-1},\qquad\Omega_{GH}(\tau):=Var\!\left(\bar{\Psi}_{GH}\right). \] By orthogonality of the ANOVA components,
Assumption (ref)(i) guarantees that the asymptotic variance of $\hat{\beta}$ is not identically zero, although some components of the variance decomposition may be absent. Assumption (ref)(ii) rules out non-Gaussian limits driven by the interaction component. Specifically, either (a) clustering along at least one dimension is sufficiently strong so that the Gaussian components $\bar{\Psi}^{(\mathrm{I})}+\bar{\Psi}^{(\mathrm{II})}$ dominate, or (b) the interaction variance $\sigma_{\mathrm{III},1\Gamma}^{2}$ is asymptotically negligible, which suppresses the potentially non-Gaussian contribution of $\bar{\Psi}^{(\mathrm{III})}$. In either case, the normalized score admits a Gaussian limit. We impose these conditions along any convergent subsequence, since the original sequence need not converge. This allows us to establish uniform validity along subsequence. Note that we do not restrict the relative growth rate between $G$ and $H$.
Theorem (ref) establishes asymptotic normality under self-normalization. This normalization accommodates the possibility that the convergence rate of $\hat{\beta}$ varies with the clustering structure. In particular, Theorem (ref) and equation (ref) together imply that the (infeasible) convergence rate of $\widehat{\beta}(\tau)$ is $r_{GH}^{1/2}$, where \[ r_{GH}\;:=\;\min\left\{ \frac{G}{\sigma_{\mathrm{I},1\Gamma}^{2}},\frac{H}{\sigma_{\mathrm{II},1\Gamma}^{2}},GH\right\} . \] Here, the rate is determined by the projection component that dominates the variance decomposition in (ref). In particular, under standard one-way clustering (e.g., along the first dimension), where $\sigma_{\mathrm{I},1\Gamma}^{2}$ is fixed and positive definite and $\sigma_{\mathrm{II},1\Gamma}^{2}=0$, the rate reduces to $G$. Under i.i.d.\ sampling, where $\sigma_{\mathrm{I},1\Gamma}^{2} =\sigma_{\mathrm{II},1\Gamma}^{2}=0$, it reduces to $GH$.
The two-way cluster-robust variance estimator for quantile regression takes the usual sandwich form \[ \widehat{\Sigma}=\widehat{D}^{-1}\widehat{\Omega}\,\widehat{D}^{-1}, \] where $\widehat{D}$ is a consistent estimator of $D(\tau)$, and $\widehat{\Omega}$ is consistent for the deterministic target ${\Omega}_{GH}$.
The matrix $D(\tau)=E\!\left[f_{e\vert X}(0\vert X_{gh})X_{gh}X_{gh}^{\top}\right]$ captures the impact of conditional heteroskedasticity through the conditional density at the target quantile. We estimate $D(\tau)$ using Powell's (nonparametric) kernel estimator, \[ \widehat{D}=\frac{1}{GH\,\ell}\sum_{g=1}^{G}\sum_{h=1}^{H}K\!\left(\frac{y_{gh}-X_{gh}^{\top}\widehat{\beta}}{\ell}\right)X_{gh}X_{gh}^{\top}, \] where $\ell>0$ is a bandwidth and $K(u)=\tfrac{1}{2}\,\mathbf{1}\{|u|\le1\}$ is the uniform kernel. Notably, the form of $\widehat{D}$ is identical to that used under i.i.d. sampling; the difference lies entirely in the dependence structure that governs its asymptotic behavior.
kato2012asymptotic establishes consistency of Powell's estimator under weak dependence. Extending this result to two-way clustered arrays is non-trivial for three reasons. First, the convergence rate of $\widehat{\beta}(\tau)$, denoted $r_{GH}$, can vary across dependence regimes, and this rate enters $\widehat{D}$ in an essential way. Second, $\widehat{D}$ itself may converge at a different regime-dependent rate, say $r_{GH,D}$, and its leading asymptotic component may change with the regime. The rates $r_{GH}$ and $r_{GH,D}$ need not coincide. If $r_{GH,D}$ is relatively small, the nominal leading term in the expansion of $\widehat{D}$ may be dominated by remainder terms driven by $r_{GH}$. Third, dependence arises along both cluster dimensions, so the analysis must disentangle the row- and column-cluster components.
Let $Q_{gh}:=vech\!\left(X_{gh}X_{gh}^{\top}\right)\in\mathbb{R}^{d(d+1)/{2}},$ and denote the condition density of $e_{gh}=e$ given subvectors of $(X_{gh}^\top,U_g,V_h)$ by $f_{e\vert X,U}(e\vert X_{gh},U_g)$, $f_{e\vert X,V}(e\vert X_{gh},V_h)$, and $f_{e\vert X,U,V}(e\vert X_{gh},U_g,V_h)$. Define
Similarly, for $j\in\{\mathrm{I},\mathrm{II},\mathrm{IV}\}$, we assume for convenience that the diagonal elements of $\sigma_{j,Q}^{2}$ share the same order and may depend on $G$ and $H$, and we use $\sigma_{j,1Q}^{2}$ to denote this order. Let $R:=\min\{G,H\}$ and define the (infeasible) rate for $\widehat{D}$ \[ r_{GH,D}:=\min\left\{ \frac{G}{\sigma_{\mathrm{I},1Q}^{2}},\frac{H}{\sigma_{\mathrm{II},1Q}^{2}},GH\ell\right\} . \] The rate $r_{GH,D}$ takes a form reminiscent of $r_{GH}$, since $\widehat{D}$ also admits a three-way decomposition into row, column, and interaction components. However, the two rates can behave quite differently, because there is no direct link between $\sigma_{\mathrm{I},1\Gamma}^{2}$ and $\sigma_{\mathrm{I},1Q}^{2}$, nor between $\sigma_{\mathrm{II},1\Gamma}^{2}$ and $\sigma_{\mathrm{II},1Q}^{2}$. Moreover, the variance components enter the two rates in different ways. For instance, if the data is i.i.d. over each $(g,h)$ cell, one can have $r_{GH}=GH$ while $r_{GH,D}=GH\ell$, with a ratio of $\ell$. In contrast, the first two components of $r_{GH,D}$ do not involve $\ell$.
We now impose the density, stronger moment, and bandwidth conditions that ensure consistency of $\widehat{D}$.
Assumptions (ref)(i)-(ii) require uniform boundedness of the conditional density around $e=0$ and the conditional fourth moments of the regressors. Assumption (ref)(iii) imposes bounded second derivative to ensure the dominated convergence. Assumption (ref)(iv) is a standard bandwidth restriction; it is the two-way clustered analogue of Assumption 3 in kato2012asymptotic. Finally, Assumption (ref)(v) imposes a nonsingularity condition to ensure that $E\!\big[Q_{gh}Q_{gh}^{\top}f_{e\vert X}(0\vert X_{gh})\big]$ is positive definite, and hence the limiting variance is not identically zero in the worst case.
Theorem (ref)(1) establishes that $\widehat{D}$ is a consistent estimator of $D(\tau)$ and provides its uniform convergence rate. Theorem (ref)(2) further gives an asymptotic linear expansion and a central limit theorem for $vech(\widehat{D})$ at the rate $r_{GH,D}^{1/2}$. The additional conditions (2)(i)-(ii) ensure that the leading stochastic term is not dominated by the remainder terms (whose size may be governed by $r_{GH}$ through $\widehat{\beta}$).
From Theorem (ref), we can see that the approximated MSE is \[ \text{AMSE}\left(\ell\right)=\frac{\ell^{4}}{36}\left\Vert Bias\right\Vert ^{2}+tr\left\{ Var\left(vech(\widehat{D})\right)\right\} , \] where $Bias:=\frac{\ell^{2}}{6}\left(E\!\big[f_{e\vert X}^{(2)}(0\vert X_{gh})\,Q_{gh}\big]\right)$. The optimal $\ell$ that minimizes AMSE is given by \[ \ell_{\text{opt}}=\left(GH\right)^{-1/5}\left(\frac{4.5\cdot tr\left(E\!\big[Q_{gh}Q_{gh}^{\top}f_{e\vert X}(0\vert X_{gh})\big]\right)}{E\!\big[f_{e\vert X}^{(2)}(0\vert X_{gh})\,Q_{gh}\big]^{\top}E\!\big[f_{e\vert X}^{(2)}(0\vert X_{gh})\,Q_{gh}\big]}\right)^{1/5} \] and we apply a rule-of-thumb bandwidth for the Gassuian location model \[ \widehat{\ell}_{\text{opt}}=\widehat{\sigma}\left(GH\right)^{-1/5}\left(\frac{4.5\cdot\frac{1}{GH}\sum_{g,h}\left\Vert Q_{gh}\right\Vert ^{2}}{\alpha\left(\tau\right)\left\Vert \frac{1}{GH}\sum_{g,h}Q_{gh}\right\Vert ^{2}}\right)^{1/5}, \] with $\widehat{\sigma}=\text{MAD}\left(\left\{ \widehat{e}_{gh}\right\} \right)/0.6745$ and $\alpha\left(\tau\right)=\left(1-\Phi^{-1}\left(\tau\right)\right)^{2}\phi\left(\Phi^{-1}\left(\tau\right)\right)$. Here, $\text{MAD}$ is the median absolute deviation, and $\Phi$ and $\phi$ are the distribution function and the density function of the standard normal distribution. In our simulations, we find that this rule-of-thumb bandwidth adapts well.
In contrast to $\widehat{D}$, the construction of $\widehat{\Omega}$ must account explicitly for two-way clustering and therefore differs from the i.i.d.\ case. Recall the (estimated) quantile score \[ \widehat{\Psi}_{gh}=X_{gh}\Bigl(\tau-\mathbf{1}\{y_{gh}\le X_{gh}^{\top}\widehat{\beta}\}\Bigr). \] We estimate $\Omega_{GH}(\tau)$ by aggregating row-, column-, and idiosyncratic components: \[ \widehat{\Omega}:=\widehat{\Omega}_{\mathrm{I}}+\widehat{\Omega}_{\mathrm{II}}+\widehat{\Omega}_{\mathrm{III,IV}}, \] where
This estimator is the quantile-regression analogue of the two-way CRVE for simple OLS estimator, cameron2011robust. The operator $\text{EVC}(\cdot)$ denotes the eigenvalue correction (e.g., projection onto the cone of positive semidefinite matrices) applied to ensure a positive semidefinite estimate.
Let $f(e_{gh},e_{gh'}\vert X_{gh},X_{gh'},U_{g},V_{h})$ and $f(e_{gh},e_{g'h}\vert X_{gh},X_{g'h},U_{g},V_{h})$ denote the conditional joint densities of $(e_{gh},e_{gh'})$ and $(e_{gh},e_{g'h})$, respectively. For integers $l,m\ge0$, define the mixed partial derivatives
We impose the following conditions for validity of $\widehat{\Omega}$.
Assumption (ref)(i) imposes standard boundedness conditions on the maximum and norm of the regressors. Assumptions (ref)(ii)--(iii) require smoothness and uniform boundedness of the relevant conditional densities and their derivatives. These conditions facilitate uniform expansions and concentration arguments under two-way dependence, and are not needed in the i.i.d.\ case. See galvao2024hac.
Theorem (ref) establishes the uniform validity of the proposed two-way CRVE. Consequently, standard large-sample inference procedures can be implemented using the quantile regression estimator $\widehat{\beta}$ together with the variance estimator $\widehat{\Sigma}$.
Note that if Assumption (ref)(ii) fails, the limiting distribution may be non-Gaussian. This case is substantially more delicate and has only recently begun to be analyzed in a systematic way; see, for example, menzel2021bootstrap, hounyo2025projection, and davezies2025analytic. In particular, menzel2021bootstrap (c.f. Proposition 4.1) provides a sharp and highly influential characterization of the asymptotic distribution for sample means. Building on this insight, we show that a closely related impossibility phenomenon extends beyond sample means to uniform inference in two-way clustered quantile regression.
Proposition (ref) establishes an impossibility result for a non-Gaussian regime. In this case, no procedure can achieve uniform consistency. Consequently, without Assumption (ref)(ii), there exists no procedure such that uniform consistency can hold over the full class of DGPs under consideration.
In this simulation section, we assess the robustness of the proposed two-way clustered quantile regression inference procedure across a range of clustering configurations. We evaluate the finite-sample performance of the proposed two-way CRVE and compare it with alternatives that only account for dependence along the $g$-dimension, the $h$-dimension, or the $(g,h)$ intersection, respectively.
For each replication, we generate a two-way array $\{(y_{gh},X_{gh})\}_{g\le G,\,h\le H}$ from
The latent components are mutually independent and i.i.d.\ standard normal. Hence both the regressor and the regression error exhibit additive two-way dependence through $(U_{g},V_{h})$ plus an idiosyncratic component. We set $\beta_j=1$ for each $j=1,\ldots,d$ and conduct inference on the null hypothesis $\mathcal{H}_0:\beta_d(\tau)=1$. We also consider specifications in which $\beta_d(\tau)$ varies with $\tau\in(0,1)$ and and test a range of $\tau$-specific null hypotheses. The results are qualitatively similar, and we therefore relegate them to the Internet Appendix.
We compute the quantile regression estimator $\widehat{\beta}_{d}(\tau)$ and the associated two-way clustered variance estimator. All results are based on $10,000$ Monte Carlo replications. By default, we set $d=10$, $G=H=50$, and $\omega_{\bullet}^{X}=\omega_{\bullet}^{e}=1$. Nominal level is $5\%$.
We compare the proposed two-way procedure (denoted CTW) with four alternatives, described in detail in the Internet Appendix:
The one-way clustered quantile bootstrap of hagemann2017cluster exhibits qualitatively similar behavior to CG and CH in our simulations. For clarity, we therefore relegate the corresponding results to the Internet Appendix.
In this DGP, both $X_{gh,j}$ and $e_{gh}$ contain additive $g$- and $h$-level components. Consequently, the score contributions relevant for inference inherit dependence in both dimensions. The proposed estimator targets this structure by combining the $g$-level, $h$-level, and $(g,h)$ components. In contrast, CG, CH, and CI omit at least one of these components. Under the present scaling, the omitted component does not vanish as $G,H$ increase and may become relatively more important as the array grows, which leads to progressively more distorted standard errors and hence worsening size (typically over-rejection) as $G,H$ increases. Rejection is based on the usual two-sided $t$-test with standard normal critical values.
\afterpage{
}
Figure (ref) reports rejection frequencies under varying clustering structures. In Panel (a), the data exhibit two-way clustering. The two-way CRVEs, CTW and CTW$_{\mathrm{II}}$, deliver stable and accurate size control as $G$ and $H$ increase, whereas the one-way CRVEs, CG and CH, substantially overreject, with rejection frequencies around $0.15$. Ignoring clustering altogether leads to the worst performance: CI overrejects increasingly as $G$ and $H$ grow. Between the two two-way procedures, CTW$_{\mathrm{II}}$ yields slightly lower rejection frequencies because it does not correct for the double-counting term, which inflates the estimated variance and therefore makes rejection harder.
Panel (b) considers one-way clustering along the second (the $H$) dimension only. In this case, CH, CTW, and CTW$_{\mathrm{II}}$ perform well, as each accounts for dependence in the $H$ dimension.
Panel (c) considers the cluster-independent design. For readability, we rescale the vertical axis because all methods yield rejection frequencies below $0.10$. Here, all procedures except CTW$_{\mathrm{II}}$ provide satisfactory size control. This indicates that, while CTW$_{\mathrm{II}}$ works well under dependence, the resulting variance inflation renders it invalid (overly conservative) when clustering is absent.
Panel (d) varies the strength of clustering dependence in the second dimension. When dependence in the $H$ dimension is weak (small $\omega^X_V,\omega_V^e$), accounting for dependence in the first dimension is more important, and CG performs well. As dependence in the $H$ dimension strengthens (large $\omega^X_V,\omega_V^e$), CH becomes more appropriate. In both settings, CI fails, whereas both CTW and CTW$_{\mathrm{II}}$ remain reliable across the full range of dependence strengths.
Figure (ref), Panel (a), further reports results for an unbalanced design in which we fix $G=50$ and vary $H$ from $20$ to $100$. We find that CH performs slightly better than CG when $H$ is small, whereas CG performs better when $H$ is large. The intuition is that when $H$ is small, each $h$-cluster contains a larger number of observations (i.e., a larger cluster size along the second dimension), so a substantial portion of the dependence is concentrated within the $H$ dimension and must be controlled; consequently, CH is more appropriate. As $H$ increases, clusters along the second dimension become smaller and less dominant, making it relatively more important to account for dependence along the first dimension, so CG improves. Panel (b) varies the number of regressors, $d$. The qualitative patterns remain essentially unchanged, indicating that the results are not sensitive to the dimension of the covariate vector.
These patterns highlight that accounting for both clustering dimensions is essential in two-way array settings. Procedures that ignore any one dimension systematically under-estimate sampling variability and over-reject. CI performs worst because it effectively treats observations as independent across $(g,h)$ and therefore misses the dominant row/column correlation. The one-way cluster methods (CG and CH) partially correct the problem by capturing dependence in a single direction, which explains why they perform better than CI, but they remain misspecified because the neglected dimension contributes non-negligibly to the score covariance. By construction, CTW targets the full two-way covariance structure, which yields stable size and a clear improvement toward the nominal level as $G$ and $H$ increase. CTW$_\mathrm{II}$ is robust to two-way clustering dependence as well, but is overly conservative when clustering is absent.
This section uses a QR framework to study how teacher-licensing restrictions affect teacher quality. Policy views on licensing are mixed. Some states have increased licensing stringency, motivated by the idea that tighter requirements can screen out lower-ability candidates and raise the left tail of the quality distribution (e.g., kraft2020teacher). Other states decreased licensing stringency, a policy choice that speaks directly to our focus on the right tail. One argument for reducing stringency is that it may attract more competitive candidates who would otherwise choose other professions (e.g., hanushek1995chooses,ballou1998case). By contrast, other work suggests that licensing requirements may have little effect on high-quality candidates (e.g., angrist2004teacher,larsen2020effect).
Let $s$ index states and $t$ index years. For each state--year cell, let $y_{st}$ denote the 90th percentile of college SAT scores among teachers in that cell, which we interpret as a measure of the right-tail (high-quality) teacher workforce. We consider the QR model
where $X_{st}$ is a measure of licensing stringency and $W_{st}$ collects controls, including school characteristics, teacher-market conditions, non-teacher labor-market conditions, education-policy controls, and political conditions. The parameter of interest is $\beta(\tau)$: a negative value, $\beta(\tau)<0$, indicates that greater stringency is associated with a lower right-tail outcome at quantile $\tau$. We use the publicly available data from larsen2020effect, who report that, on average (based on OLS), licensing stringency does not affect high-quality candidates.
Table (ref) reports $\widehat{\beta}(\tau)$ for a grid of quantiles together with $p$-values computed under several CRVE choices. The main evidence of an right-tail effect arises at low $\tau$. At $\tau=0.10$, $\widehat{\beta}(0.10)=-0.0998$, and the CTW $p$-value is $0.0004$, indicating a statistically significant negative association at the 1% level. At $\tau=0.20$, $\widehat{\beta}(0.20)=-0.0870$ with a CTW $p$-value of $0.0039$, again significant at conventional levels. At $\tau=0.30$, the point estimate remains negative ($\widehat{\beta}(0.30)=-0.0668$), but inference becomes sensitive to the variance estimator: CI and CG reject at 5%, whereas CH and CTW are borderline (around the 10% level) and CTW$_{\mathrm{II}}$ is more conservative. For quantiles $\tau\in\{0.40,\ldots,0.90\}$, the estimates are close to zero and none of the CRVEs yield statistically significant effects.
Overall, emphasizing the two-way robust CTW inference, the results suggest that licensing stringency may not affect the right tail on average, consistent with larsen2020effect, but the effect can be heterogeneous across quantiles. In particular, the negative association is concentrated in the lower part of the conditional distribution of $y_{st}$ (roughly $\tau\le 0.20$), consistent with the presence of a margin of high-quality candidates whose occupational choice is sensitive to licensing costs. For higher quantiles, we find little evidence that stringency discourages right-tail teacher quality at 5% significance level.
This paper develops a unified large-sample theory and practical inference procedures for linear quantile regression under two-way clustering. The key challenge is that both the non-smooth quantile score and the two-way dependence invalidate standard arguments, and, moreover, the effective convergence rate of the quantile regression estimator can vary across dependence regimes. To address these issues, we work within a separately exchangeable array framework and employ a projection-based decomposition that isolates row, column, interaction, and idiosyncratic components. This structure yields a transparent variance identity and an asymptotic distribution theory that adapts to regime-dependent normalizations.
Building on the limit theory, we propose a feasible two-way cluster-robust sandwich covariance estimator. We show that both the “bread” component (a kernel estimator of the conditional density at the target quantile) and the “meat” component (an estimator of the covariance of the sample score that aggregates row and column contributions) are consistent under appropriate smoothness and moment conditions. The resulting procedure is asymptotically valid in the Gaussian regimes, with a proof that explicitly tracks how regime-dependent rates and two-way dependence alter the relative magnitude of leading terms and remainder terms.
Moreover, we clarify the intrinsic limits of uniform inference under two-way clustering. When the interaction component remains asymptotically non-negligible while clustering variation along both dimensions is bounded, the limiting distribution can be non-Gaussian, and uniform consistency over the full model class may be unattainable without additional restrictions.
The simulation results further demonstrate the necessity of using a two-way cluster-robust variance estimator when two-way clustering is present. They also highlight the robustness of the two-way procedure across a range of dependence structures: it remains valid under varying levels of clustering dependence in two dimensions, and even in the absence of within-cluster dependence. In an empirical application, we find that the effect of teacher-licensing stringency on teacher quality is heterogeneous across the distribution. Specifically, tighter licensing requirements primarily affect teachers in the bottom $20\%$, who are plausibly closer to the margin of selecting into alternative careers. In contrast, we find little evidence that licensing stringency discourages high-quality teachers at higher quantiles.
Overall, the paper closes a theoretical gap for quantile regression with two-way clustered data and offers easy-to-implement inference tools that are directly applicable in empirical settings where multi-dimensional clustering is unavoidable.