The exact contents of citations.db main_text.text for this paper — one flattened LaTeX string, title through conclusion, appendix excluded, unmodified except for removing email addresses. This is what our citation measures are computed over.
52,508 characters
Two-way Clustering Robust Variance Estimator in Quantile Regression Models
\maketitle
\begin{abstract}
We study inference for linear quantile regression with two-way clustered data. Using a separately exchangeable array framework and a projection decomposition of the quantile score, we characterize regime-dependent convergence rates and establish a self-normalized Gaussian approximation. We propose a two-way cluster-robust sandwich variance estimator with a kernel-based density ``bread'' and a projection-matched ``meat'', and prove consistency and validity of inference in Gaussian regimes. We also show an impossibility result for uniform inference in a non-Gaussian interaction regime.
\bigskip
\textbf{JEL Classification}: C15, C23, C31, C80
\medskip
\textbf{Keywords}: Clustered
data, cluster-robust variance estimator, two-way clustering, quantile regression.
\end{abstract}
\vspace*{-0.5cm}
\vfill
\thispagestyle{empty} \pagebreak
\section{Introduction}
\label{sec:introduction}
Quantile regression (QR), introduced by \citet{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 \citep{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 \citep[e.g.,][]{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., \citealt{cameron2011robust,mackinnon2021wild}) and to one-way clustered quantile regression
(e.g., \citealt{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 \emph{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, \citet{chiang2024extremal} study extremal quantiles
under two-way clustered dependence, focusing on rare-event estimation in
two-way clustered data. While \citet{menzel2021bootstrap} show that sample means may exhibit
non-Gaussian limits under two-way clustering, \citet{chiang2024extremal}
demonstrate that extremal quantiles can remain robust even in degenerate
dependence regimes.
Our paper complements this line of work by
focusing on \emph{interior} quantiles: while \citet{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 \citep{kato2012asymptotic}, CRVE and pigeonhole bootstrap results for GMM (e.g. quantile IV) under multiway clustering \citep{davezies2018asymptotic}, and recent advances in weak-dependence-robust covariance estimation for quantile regression \citep{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 \emph{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 \citet{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 \citet{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{sec:model}
introduces the two-way QR framework and the projection decomposition
and develops the regime-adaptive limit theory for $\widehat{\beta}(\tau)$.
Section~\ref{sec:variance} 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{sec:mc} presents Monte
Carlo evidence. Section~\ref{sec:empirical} studies an application to teacher licensing. Section~\ref{sec:conclusion} concludes. Technical
proofs and additional results are deferred to the appendix.
\section{Two-Way Clustering in Quantile Regression}
\label{sec:model}
\subsection{Model Setting}
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
\begin{equation}
Q_{y_{ghi}}(\tau\vert X_{ghi})=X_{ghi}^{\top}\beta_{0}(\tau),\qquad g=1,\dots,G,\;h=1,\dots,H,\;i=1,\dots,N_{gh},\label{eq:qr-model}
\end{equation}
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
\begin{equation}
\psi_{ghi}(\beta,\tau):=X_{ghi}\Bigl(\tau-\mathbf{1}\{y_{ghi}\le X_{ghi}^{\top}\beta\}\Bigr),\qquad\Psi_{ghi}(\tau):=\psi_{ghi}\bigl(\beta_{0}(\tau),\tau\bigr).\label{eq:qr-score}
\end{equation}
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, \citet{davezies2021empirical},
\citet{mackinnon2021wild}, and \citet{chiang2023standard}.
\begin{assumption}[Two-way clustered data with the AHK representation]
\label{ass:ahk} There exist measurable functions $\Gamma$ such that
\[
(y_{gh},X_{gh})=\Gamma(U_{g},V_{h},W_{gh}),
\]
where $\{U_{g}\}_{g\ge1}$, $\{V_{h}\}_{h\ge1}$, and $\{W_{gh}\}_{g,h\ge1}$
are mutually independent sequences of i.i.d.\ random variables. Without
loss of generality, each latent variable is uniformly distributed
on $[0,1]$. The function $\Gamma$ may vary with $(G,H)$, allowing for triangular-array sequences of DGPs.
\end{assumption}
Under Assumption~\ref{ass:ahk}, 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
\begin{equation}
\Psi_{gh}(\tau)=E[\Psi_{gh}(\tau)]+\Psi^{(\mathrm{I})}(U_{g},\tau)+\Psi^{(\mathrm{II})}(V_{h},\tau)+\Psi^{(\mathrm{III})}(U_{g},V_{h},\tau)+\Psi^{(\mathrm{IV})}(U_{g},V_{h},W_{gh},\tau),\label{eq:anova_decomposition_main}
\end{equation}
where
\begin{align*}
\Psi^{(\text{I})}(U_{g},\tau) & :=E[\Psi_{gh}(\tau)\vert U_{g}]-E[\Psi_{gh}(\tau)],\\
\Psi^{(\text{II})}(V_{h},\tau) & :=E[\Psi_{gh}(\tau)\vert V_{h}]-E[\Psi_{gh}(\tau)],\\
\Psi^{(\text{III})}(U_{g},V_{h},\tau) & :=E[\Psi_{gh}(\tau)\vert U_{g},V_{h}]-E[\Psi_{gh}(\tau)]-\Psi^{(\text{I})}(U_{g},\tau)-\Psi^{(\text{II})}(V_{h},\tau),\\
\Psi^{(\text{IV})}(U_{g},V_{h},W_{gh},\tau) & :=\Psi_{gh}(\tau)-E[\Psi_{gh}(\tau)\vert U_{g},V_{h}].
\end{align*}
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 \citet{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.
\begin{assumption}[Moments, smoothness, and nonsingularity]\label{as:moment}
(i) $E[\Psi_{gh}]=0$, $E\|X_{gh}\|^{4}<\infty$, and $E(X_{gh}X_{gh}^{\top})$
is nonsingular. (ii) The map $e\mapsto f_{e\vert X}(e\vert x)$ is
twice continuously differentiable for every $x$, and $\sup_{e,x}\big|f_{e\vert X}(e\vert x)\big|<\infty$
as well as $\sup_{e,x}\big|f_{e\vert X}^{(1)}(e\vert x)\big|<\infty$.
(iii) The conditional density at zero is uniformly bounded away from
zero: $\inf_{x}f_{e\vert X}(0\vert x)>0$. (iv) The Jacobian matrix
$E\!\left[f_{e\vert X}(0\vert X_{gh})\,X_{gh}X_{gh}^{\top}\right]$
and the score variance $E\!\left[\Psi_{gh}\Psi_{gh}^{\top}\right]$
are positive definite. (v) $\beta_{0}(\tau)$ lies in the interior
of a compact parameter space $\Theta$. \end{assumption}
\subsection{Asymptotic Distribution}
\label{sec:asymptotics} 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 \eqref{eq:anova_decomposition_main}, we decompose $\bar{\Psi}_{GH}$
as
\begin{align*}
\bar{\Psi}_{GH} & =\frac{1}{G}\sum_{g=1}^{G}\Psi_{g}^{(\text{I})}+\frac{1}{H}\sum_{h=1}^{H}\Psi_{h}^{(\text{II})}+\frac{1}{GH}\sum_{g=1}^{G}\sum_{h=1}^{H}\Bigl(\Psi_{gh}^{(\text{III})}+\Psi_{gh}^{(\text{IV})}\Bigr)\\
& :=\bar{\Psi}^{(\text{I})}+\bar{\Psi}^{(\text{II})}+\bar{\Psi}^{(\text{III})}+\bar{\Psi}^{(\text{IV})}.
\end{align*}
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,
\begin{equation}
\Omega_{GH}(\tau)=\frac{1}{GH}\Bigl(H\sigma_{\text{I},\Gamma}^{2}+G\sigma_{\text{II},\Gamma}^{2}+\sigma_{\text{III},\Gamma}^{2}+\sigma_{\text{IV},\Gamma}^{2}\Bigr).\label{eq: variance of score}
\end{equation}
\begin{assumption}[Orders of variance components]\label{as:order of variance}
(i) The total variance does not vanish:
\[
\liminf_{G,H\to\infty}\Bigl(H\sigma_{\mathrm{I},1\Gamma}^{2}+G\sigma_{\mathrm{II},1\Gamma}^{2}+\sigma_{\mathrm{III},1\Gamma}^{2}+\sigma_{\mathrm{IV},1\Gamma}^{2}\Bigr)\;>\;0.
\]
(ii)Along any subsequence indexed by $(G_n,H_n)$ for which
$\bigl(H_n\sigma_{\mathrm{I},1\Gamma}^{2},G_n\sigma_{\mathrm{II},1\Gamma}^{2},\sigma_{\mathrm{III},1\Gamma}^{2},\sigma_{\mathrm{IV},1\Gamma}^{2}\bigr)$
converges in $[0,\infty]^{4}$, at least one of the following holds:
\[
\text{(a)}\quad H_n\sigma_{\mathrm{I},1\Gamma}^{2}+G_n\sigma_{\mathrm{II},1\Gamma}^{2}\to\infty,
\qquad\text{or}\qquad
\text{(b)}\quad \sigma_{\mathrm{III},1\Gamma}^{2}\to 0 .
\]
\end{assumption}
Assumption \ref{as:order of variance}(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{as:order of variance}(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$.
\begin{theorem}\label{thm:1} Let $\mathcal{B}_{0}$ denote the collection
of DGPs $\Gamma$ that satisfy Assumptions \ref{ass:ahk}--\ref{as:order of variance}.
Then
\[
\Sigma_{GH}^{-1/2}\bigl(\hat{\beta}-\beta_{0}(\tau)\bigr)\overset{d}{\to}\mathcal{N}\!\left(0,\mathbf{I}_{d}\right)
\]
uniformly over $\Gamma\in\mathcal{B}_{0}$. \end{theorem}
Theorem \ref{thm:1} 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{thm:1} and equation \eqref{eq: variance of score} 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 \eqref{eq: variance of score}. 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$.
\section{Cluster-Robust Variance Estimator (CRVE)}
\label{sec:variance} 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}$.
\subsection{Estimating $D(\tau)$.}
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.
\citet{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
\begin{align*}
\sigma_{\mathrm{I},Q}^{2} & :=Var\!\Big(E\!\big[\,Q_{gh}f_{e\vert X,U}(0\vert X_{gh},U_{g})\,\big|\,U_{g}\big]\Big),\\
\sigma_{\mathrm{II},Q}^{2} & :=Var\!\Big(E\!\big[\,Q_{gh}f_{e\vert X,V}(0\vert X_{gh},V_{h})\,\big|\,V_{h}\big]\Big).
\end{align*}
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}$.
\begin{assumption}[Density and bandwidth]\label{as:bandwidth}
(i) There exist $\varepsilon_{0}>0$ and constants $0<c<C<\infty$
such that, uniformly over $(x,u,v)$ and all $|e|\le\varepsilon_{0}$,
$c\le f_{e\vert X,U,V}(e\vert x,u,v)\le C.$ (ii) $E\!\left(\|X_{gh}\|^{4}\vert U_{g},V_{h}\right)<\infty$
uniformly over $(U_{g},V_{h})$. (iii) $\sup_{e,x}\bigl|f_{e\vert X}^{(2)}(e\vert x)\bigr|<\infty$
(iv) As $R\to\infty$, $\ell\to0$ and $R\ell^{2}/\log R\to\infty$.
(v) $E\!\big[Q_{gh}Q_{gh}^{\top}\big]$ is positive definite. \end{assumption}
Assumptions \ref{as:bandwidth}(i)-(ii) require uniform boundedness
of the conditional density around $e=0$ and the conditional fourth
moments of the regressors. Assumption \ref{as:bandwidth}(iii) imposes
bounded second derivative to ensure the dominated convergence. Assumption
\ref{as:bandwidth}(iv) is a standard bandwidth restriction; it is
the two-way clustered analogue of Assumption~3 in \citet{kato2012asymptotic}.
Finally, Assumption \ref{as:bandwidth}(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.
\begin{theorem}\label{thm:Jacobian} Let $\mathcal{B}_{1}$ denote
the collection of DGPs $\Gamma$ that satisfy Assumptions \ref{ass:ahk}--\ref{as:bandwidth}.
Then the following statements hold uniformly over $\Gamma\in\mathcal{B}_{1}$.
\begin{enumerate}[label=(\arabic*)]
\item \textbf{Consistency and rate.}
\[
\widehat{D}-D(\tau)=O_{P}\!\left(r_{GH}^{-1/2}\ell^{-1/2}+\ell^{2}\right)=o_{P}(1).
\]
\item \textbf{Asymptotic normality.} Suppose, in addition, that at least
one of the following conditions holds: (i) $\sigma_{\mathrm{I},1\Gamma}^{2}/(\ell\sigma_{\mathrm{I},1Q}^{2})=O(1)$
and $\sigma_{\mathrm{II},1\Gamma}^{2}/(\ell\sigma_{\mathrm{II},1Q}^{2})=O(1)$;
or (ii) $H\sigma_{\mathrm{I},1\Gamma}^{2}+G\sigma_{\mathrm{II},1\Gamma}^{2}=O\left(1\right)$.
Then
\begin{align*}
V_D^{-1/2}\Bigg(vech(\widehat{D})-vech\!\big(D(\tau)\big)-\frac{\ell^{2}}{6}E\!\big[f_{e\vert X}^{(2)}(0\vert X_{gh})\,Q_{gh}\big]+o(\ell^{2})\Bigg)\overset{d}{\to}\mathcal{N}\!\left(\bm{0}_{\frac{d(d+1)}{2}\times1},\,\mathbf{I}_{\frac{d(d+1)}{2}}\right),
\end{align*}
where
\[
V_{D}=\frac{\sigma_{\mathrm{I},Q}^{2}}{G}+\frac{\sigma_{\mathrm{II},Q}^{2}}{H}+\frac{1}{2GH\ell}\,E\!\big[Q_{gh}Q_{gh}^{\top}f_{e\vert X}(0\vert X_{gh})\big].
\]
\end{enumerate}
\end{theorem}
Theorem \ref{thm:Jacobian}(1) establishes that $\widehat{D}$ is
a consistent estimator of $D(\tau)$ and provides its uniform convergence
rate. Theorem \ref{thm:Jacobian}(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}$).
\begin{remark}
Unlike the score-based limit theory, no extra restriction
on the potentially non-Gaussian interaction component is required
here for the kernel-based estimator, since its contribution is already of smaller order relative
to the dominant Gaussian terms in $\widehat{D}$.
\end{remark}
From Theorem \ref{thm:Jacobian}, 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.
\subsection{Consistency of Quantile Regression CRVE}
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
\begin{align*}
\widehat{\Omega}_{\mathrm{I}} & :=\text{EVC}\!\left(\frac{1}{G^{2}H^{2}}\sum_{g=1}^{G}\sum_{h=1}^{H}\sum_{\substack{h'=1\\
h'\neq h
}
}^{H}\widehat{\Psi}_{gh}\widehat{\Psi}_{gh'}^{\top}\right),\\
\widehat{\Omega}_{\mathrm{II}} & :=\text{EVC}\!\left(\frac{1}{G^{2}H^{2}}\sum_{h=1}^{H}\sum_{g=1}^{G}\sum_{\substack{g'=1\\
g'\neq g
}
}^{G}\widehat{\Psi}_{gh}\widehat{\Psi}_{g'h}^{\top}\right),\\
\widehat{\Omega}_{\mathrm{III,IV}} & :=\frac{1}{G^{2}H^{2}}\sum_{g=1}^{G}\sum_{h=1}^{H}\widehat{\Psi}_{gh}\widehat{\Psi}_{gh}^{\top}.
\end{align*}
This estimator is the quantile-regression analogue of the two-way
CRVE for simple OLS estimator, \citet{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
\begin{align*}
f^{(l,m)}(e_{gh},e_{gh'}\vert X_{gh},X_{gh'}) & :=\frac{\partial^{\,l+m}}{\partial e_{gh}^{\,l}\,\partial e_{gh'}^{\,m}}\,f(e_{gh},e_{gh'}\vert X_{gh},X_{gh'}),\\
f^{(l,m)}(e_{gh},e_{g'h}\vert X_{gh},X_{g'h}) & :=\frac{\partial^{\,l+m}}{\partial e_{gh}^{\,l}\,\partial e_{g'h}^{\,m}}\,f(e_{gh},e_{g'h}\vert X_{gh},X_{g'h}).
\end{align*}
We impose the following conditions for validity of $\widehat{\Omega}$.
\begin{assumption}[Strong moments and smoothness]\label{as:moment-strong}
There exist a constant $C_{1}>0$ and integrable envelope
functions $D_{1}(\cdot)$ and $D_{2}(\cdot)$ such that:
\begin{enumerate}[label=(\roman*)]
\item \label{as:ms-i} $\max_{g\le G}\max_{h\le H}\|X_{gh}\|\le C_{1}R^{1/8}$ a.s. and $\sup_{g,h}E(\|X_{gh}\|^6\vert U_g,V_h)<\infty$ a.s.
\item \label{as:ms-iii} The conditional joint densities are uniformly bounded:
\[
\sup_{e_{1},e_{2},x_{1},x_{2},U_{g},V_{h},U_{g'},V_{h'}}\bigl|f(e_{1},e_{2}\vert x_{1},x_{2},U_{g},V_{h},U_{g'},V_{h'})\bigr|<\infty,
\]
where $(e_{1},e_{2},x_{1},x_{2})$ denotes either $(e_{gh},e_{gh'},X_{gh},X_{gh'})$
or $(e_{gh},e_{g'h},X_{gh},X_{g'h})$.
\item \label{as:ms-iv} For $l,m\in\{1,2\}$,
\[
\sup_{e_{2},x_{1},x_{2}}\bigl|f^{(l,0)}(e_{1},e_{2}\vert x_{1},x_{2})\bigr|\le D_{1}(e_{1}),\qquad\sup_{e_{2},x_{1},x_{2}}\bigl|f^{(0,m)}(e_{1},e_{2}\vert x_{1},x_{2})\bigr|\le D_{2}(e_{1}),
\]
for both pairs $(e_{1},e_{2},x_{1},x_{2})=(e_{gh},e_{gh'},X_{gh},X_{gh'})$
and $(e_{gh},e_{g'h},X_{gh},X_{g'h})$.
\end{enumerate}
\end{assumption} Assumption \ref{as:moment-strong}(i) imposes standard
boundedness conditions on the maximum and norm of the regressors. Assumptions
\ref{as:moment-strong}(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 \citet{galvao2024hac}.
\begin{theorem}\label{thm:2-1} Let $\mathcal{B}_{2}$ denote the
collection of DGPs $\Gamma$ that satisfy Assumptions \ref{ass:ahk}--\ref{as:moment-strong}.
Then, uniformly over $\Gamma\in\mathcal{B}_{2}$,
\[
\Omega_{GH}(\tau)^{-1}\widehat{\Omega}\ \overset{P}{\to}\ \mathbf{I}_{d},\qquad\text{and}\qquad\widehat{\Sigma}^{-1/2}\bigl(\hat{\beta}-\beta_{0}(\tau)\bigr)\ \overset{d}{\to}\ \mathcal{N}\!\left(0,\mathbf{I}_{d}\right).
\]
\end{theorem} Theorem \ref{thm:2-1} 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{as:order of variance}(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, \citet{menzel2021bootstrap}, \citet{hounyo2025projection}, and \citet{davezies2025analytic}. In particular, \citet{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.
\begin{proposition}[Impossibility of uniform consistency]\label{prop:impossibility_uniform}
Let $\mathcal{B}_{3}$ denote the class of DGPs $\Gamma$ satisfying
Assumptions \ref{ass:ahk}--\ref{as:moment-strong}, except that
Assumption \ref{as:order of variance}(ii) is replaced by
\[
H\sigma_{\mathrm{I},1\Gamma}^{2}+G\sigma_{\mathrm{II},1\Gamma}^{2}=O(1)\quad\text{and}\quad\limsup_{G,H\to\infty}\sigma_{\mathrm{III},1\Gamma}^{2}>0.
\]
Let $\mathcal{E}$ be the collection of all measurable maps of the
observed sample $\{y_{gh}^{(\Gamma)},X_{gh}^{(\Gamma)}\}_{g\le G,h\le H}$ generated
by $\Gamma$. Then there exists $\varepsilon>0$ and $\delta>0$ such that
\[
\liminf_{G,H\to\infty}\ \inf_{\widehat{E}\in\mathcal{E}}\ \sup_{\Gamma\in\mathcal{B}_{3}}P_\Gamma\left(\sup_{t\in\mathbb{R}^{d}}\left\vert P_{\Gamma}\!\left(\sqrt{GH}\,(\widehat{\beta}-\beta)\le t\right)-\ \widehat{E}\left(\{y_{gh}^{(\Gamma)},X_{gh}^{(\Gamma)}\}_{g\le G,h\le H};t\right)\right\vert>\varepsilon\right) \ge \delta,
\]
\end{proposition}
Proposition \ref{prop:impossibility_uniform} establishes an impossibility
result for a non-Gaussian regime. In this case, no procedure can achieve
uniform consistency. Consequently, without Assumption \ref{as:order of variance}(ii),
there exists no procedure such that uniform consistency can hold over
the full class of DGPs under consideration.
\section{Monte Carlo simulation}
\label{sec:mc}
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
\begin{align*}
y_{gh} & =\beta_{1}+\sum_{j=2}^{d}\beta_{j}X_{gh,j}+e_{gh},\\
X_{gh,j} & =\omega_{U}^{X}U_{g}^{X,j}+\omega_{V}^{X}V_{h}^{X,j}+\omega_{W}^{X}W_{gh}^{X,j},\\
e_{gh} & =\omega_{U}^{e}U_{g}^{e}+\omega_{V}^{e}V_{h}^{e}+\omega_{W}^{e}W_{gh}^{e}.
\end{align*}
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 \textbf{CTW}) with
four alternatives, described in detail in the Internet Appendix:
\begin{itemize}
\item \textbf{CG (cluster-$g$ only).} A one-way clustered inference method
that treats $g$ as the only clustering dimension and ignores dependence
across $h$.
\item \textbf{CH (cluster-$h$ only).} A one-way clustered inference method
that treats $h$ as the only clustering dimension and ignores dependence
across $g$.
\item \textbf{CI (intersection-only).} An i.i.d.-style inference method
that effectively uses only the $(g,h)$ intersection component and
ignores both two-way additive components.
\item \textbf{CTW$_{\text{II}}$ (two-way cluster without intersection correction).
}A two-way clustered inference procedure that enforces positive semidefiniteness
without using EVC, but does not correct for the ``double-counting''
of the intersection component.
\end{itemize}
The one-way clustered quantile bootstrap of \citet{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 \emph{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{ \begin{landscape}
\begin{figure}[h!]
\centering \begin{subfigure}[t]{0.49\hsize} \subcaption{Two-way Clustering}\includegraphics[width=1\textwidth]{diff_GH_two_cluster.png}
\end{subfigure} \begin{subfigure}[t]{0.49\hsize} \subcaption{One-way Clustering, $\omega_U^X=\omega_U^e=0$} \includegraphics[width=1\textwidth]{diff_GH_one_cluster.png}
\end{subfigure}
\begin{subfigure}[t]{0.49\hsize} \subcaption{ Independence, $\omega_U^X=\omega_U^e=\omega_V^X=\omega_V^e=0$}\includegraphics[width=1\textwidth]{diff_GH_no_cluster.png}
\end{subfigure} \begin{subfigure}[t]{0.49\hsize} \subcaption{Varying Clustering Dependence $\omega_V^X$ and $\omega_V^e$} \includegraphics[width=1\textwidth]{diff_omega_two_cluster.png}
\end{subfigure}
\caption{\textbf{Rejection frequency under varying levels of clustering dependence}
The default setting is $d=10$, $G=H=50$, and $\omega_{\bullet}^{X}=\omega_{\bullet}^{e}=1$. Results are based on 10,000 Monte Carlo replicates.
The predetermined significance level is 5\%.}
\label{fig: rej frequency 1}
\end{figure}
\end{landscape} }
Figure~\ref{fig: rej frequency 1} 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{fig: rej frequency 2}, 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 \emph{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.
\begin{figure}[t!]
\centering \begin{subfigure}[t]{0.49\hsize} \subcaption{$G=50$, Varying $H$}\includegraphics[width=1\textwidth]{diff_GG_two_cluster.png}
\end{subfigure} \begin{subfigure}[t]{0.49\hsize} \subcaption{Varying $d$} \includegraphics[width=1\textwidth]{diff_K_two_cluster.png}
\end{subfigure}
\caption{\textbf{Rejection frequency under different structures}
The default setting is $d=10$, $G=H=50$, and $\omega_{\bullet}^{X}=\omega_{\bullet}^{e}=1$. Results are based on 10,000 Monte Carlo replicates.
The predetermined significance level is 5\%.}
\label{fig: rej frequency 2}
\end{figure}
\section{Empirical Studies}
\label{sec:empirical}
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., \citealt{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., \citealt{hanushek1995chooses,ballou1998case}).
By contrast, other work suggests that licensing requirements may have
little effect on high-quality candidates (e.g.,
\citealt{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
\begin{equation}\label{eq:qr_top_tail}
Q_{y_{st}\mid X_{st},W_{st}}(\tau)
=
\alpha(\tau)
+ X_{st}\beta(\tau)
+ W_{st}^{\top}\gamma(\tau),
\qquad \tau\in(0,1),
\end{equation}
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 \citet{larsen2020effect}, who
report that, on average (based on OLS), licensing stringency does not affect
high-quality candidates.
\begin{table}[t]
\centering
\caption{Effects of licensing stringency and $p$-values under different CRVEs.}
\label{tab:qr_pvalues}
\begin{tabular}{lccccccccc}
\hline\hline
$\tau$ & $0.10$ & $0.20$ & $0.30$ & $0.40$ & $0.50$ & $0.60$ & $0.70$ & $0.80$ & $0.90$ \\
\hline
$\hat\beta(\tau)$ & -0.0998 & -0.0870 & -0.0668 & -0.0295 & -0.0277 & -0.0107 & 0.0033 & 0.0143 & 0.0164 \\
CI & 0.0001 & 0.0014 & 0.0138 & 0.2467 & 0.2725 & 0.7407 & 0.9657 & 0.8561 & 0.9641 \\
CG & 0.0000 & 0.0000 & 0.0171 & 0.2747 & 0.3864 & 0.7723 & 0.9735 & 0.8919 & 0.9724 \\
CH & 0.0035 & 0.0125 & 0.0847 & 0.4493 & 0.4961 & 0.8447 & 0.9790 & 0.9095 & 0.9769 \\
CTW & $\bm{0.0004}$ & $\bm{0.0039}$ & $\bm{0.0898}$ & 0.4610 & 0.5399 & 0.8524 & 0.9812 & 0.9207 & 0.9797 \\
CTW$_{\mathrm{II}}$ & 0.0092 & 0.0320 & 0.1624 & 0.5340 & 0.5925 & 0.8711 & 0.9835 & 0.9305 & 0.9823 \\
\hline\hline
\end{tabular}\label{tab:empirical}
\end{table}
Table~\ref{tab:empirical} 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
\emph{on average}, consistent with \citet{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.
\section{Conclusion}
\label{sec:conclusion}
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.
\clearpage{}