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.
64,216 characters · 11 sections · 26 citation commands
Extremal Quantiles under Two-Way Clustering
Quantile regression has become increasingly important in economics and related fields for analyzing the effects of covariates across different parts of the conditional distribution. Many key applications focus on the study of rare events, or quantiles in tails. Notable examples include value-at-risk analysis in financial economics Chernozhukov2001, Engle2004, the analysis of highest bids in auctions chernozhukov2004likelihood, Hirano2003, the investigation of birth weights among live infants \citep*{Chernozhukov2011, Zhang2018, Sasaki2022, Sasaki2023}, and the growth-at-risk (GaR) and related analyses adrian2019vulnerable, adrian2022term, kiley2022unemployment, kiley2024growth, lopez2024inflation, among others.
The conventional large-sample theory for quantile regression encounters difficulties in the extreme tails of the distribution due to data sparsity. As a result, traditional inference tools, such as the analytical variance estimator (koenker2005quantile, Section 3.4) and the pivotal bootstrap \citep*{chernozhukov2009finite}, are not directly applicable in these regions. However, advances in modern extreme value theory have helped overcome this limitation for unconditional quantiles \citep*{Leadbetter1983, Resnick1987, Embrechts1997}. Chernozhukov2005 extends this theoretical framework to estimation and inference for extremal conditional quantiles. Further developments include asymptotic theory and inference methods for quantile treatment effects in the tails, as explored by \citet*{d2018extremal}, Zhang2018, \citet*{Deuber2021}, and \citet*{Sasaki2023}.
At the same time, applied economic research increasingly involves data exhibiting complex clustering structures. In particular, two-way clustering is frequently encountered in empirical work. Common examples include panel data, matched employer–employee data, matched student–teacher data, scanner data indexed by stores and products, market share data indexed by markets and products, and growth or development data indexed by ethnicity and geographic units. These examples underscore the prevalence of multi-way clustering in real-world datasets. It is essential to account for the complications that arise under two-way clustering in order to conduct reliable and valid econometric analysis.
The seminal work by \citet*{Cameron2011} introduces methods for two-way cluster-robust inference, which has had a substantial influence on the field. Menzel2021 provides a formal analysis of the asymptotic properties and bootstrap validity in the presence of two-way clustering. \citet*{Davezies2021} advance the development of empirical process theory for such settings, and \citet*{MacKinnon2021} propose wild bootstrap procedures that explicitly account for clustering effects. In a related direction, Graham2020 studies sparse network asymptotics, with a focus on logistic regression models.
Our work builds on these foundational contributions and focuses on the estimation of rare events in data sampled under two-way clustering. This focus is not only novel and empirically relevant, but also theoretically crucial, as it addresses the distinctive challenges and provides new insights. Two-way clustering introduces dependence structures that complicate statistical analysis, particularly in distributional tails. By targeting tail estimation, we deepen understanding of the behavior of extreme observations in such settings. Notably, while Menzel2021 shows that the asymptotic distribution of sample means and related statistics may deviate from Gaussianity under two-way clustering, we find that extremal quantiles exhibit robustness to such deviations, even in degenerate cases.
The techniques underlying our asymptotic theory are closely related to prior research developed under alternative asymptotic frameworks. Our approach draws parallels with the analysis of network subgraph moments by \citet*{Bickel2011}, dyadic density estimation by \citet*{Graham2022, chiang2023empirical, cattaneo2023uniform}, and the semiparametric asymptotic regimes studied by \citet*{cattaneo2014small, cattaneo2018alternative}, in which both linear and quadratic components play essential roles. These studies provide valuable insights and foundational frameworks upon which our work builds. In addition, this paper contributes to the growing literature on inference for quantile regression with cluster-sampled data, including hagemann2017cluster, hagemann2024inference, and connects with the asymptotic theory for various regression models in dyadic data contexts, such as graham2020dyadic and sassi2023ordinary.
The remainder of the paper is organized as follows. Section 2 analyzes unconditional quantiles under two-way clustering, presenting both the theoretical framework and supporting simulation results. Section 3 extends these results to quantile regression. Section 4 provides an empirical application to growth-at-risk (GaR). Section 5 concludes. All mathematical proofs are collected in the appendices.
We begin our analysis by examining the limiting distribution of unconditional quantiles in the tails. In the subsequent section, we extend this baseline framework to extremal quantile regression. The core ideas and methodological approach remain consistent across both settings. We first focus on unconditional quantiles because this setting allows us to more transparently illustrate how Gaussianity is achieved.
Consider a random variable $Y$ taking values in $\mathbb{R}$, and let its distribution function be denoted by $F_Y$. Our focus is on the extremal quantile $$ F_Y^{-1}(\tau) = \inf \left\{ y : F_Y(y) > \tau \right\} $$ for small values of $\tau$ (close to $0$). A symmetric argument applies to the case when $\tau$ is close to $1$.
Suppose we observe data $\{Y_{it} : 1 \le i \le N,\ 1 \le t \le T\}$, where $Y_{it} \overset{d}{=} Y$ is generated according to $$ Y_{it} = f(\alpha_i, \gamma_t, \varepsilon_{it}) $$ for some Borel measurable function $f$, where $\alpha_1, \ldots, \alpha_N$, $\gamma_1, \ldots, \gamma_T$, and $\varepsilon_{11}, \ldots, \varepsilon_{NT}$ are mutually independent random vectors. Each of the latent components, $\alpha_i$, $\gamma_t$, and $\varepsilon_{it}$, is allowed to be vector-valued with possibly countably infinite dimensions. This setup accommodates dependence within each row $\left(Y_{i1}, \ldots, Y_{iT}\right)$ and within each column $\left(Y_{1t}, \ldots, Y_{Nt}\right)$ via the shared latent factors.
The unconditional $\tau$-th quantile $\beta(\tau)$ of $Y$ is estimated by the sample quantile
where $\rho_\tau$ defined by $\rho_\tau(u) = u(\tau - \mathbbm{1}\{u < 0\})$ denotes the check function.
To conduct inference for $\beta(\tau)$ based on $\widehat{\beta}(\tau)$, we study the limiting distribution of the statistic $$ \widehat{Z}_{NT} = r_{NT} \left(\widehat{\beta}(\tau) - \beta(\tau)\right), $$ where $\{r_{NT}\}$ is a sequence that governs the rate of convergence. Our focus is on the asymptotic behavior of $\widehat{Z}_{NT}$, and consequently of $\widehat{\beta}(\tau)$, in the left tail of the distribution of $Y$, as $\tau \to 0$.
For ease of writing, define
and denote its variance by $$ \sigma^2_{NT} = \mathrm{Var} \left( W_{NT}(\tau) \right). $$ Let $m > 1$ be fixed. With these notations, write the statistic as $$ \widehat{Z}_{NT} \equiv \frac{a_{NT}}{\sigma_{NT}} \left( \widehat{\beta}(\tau) - \beta(\tau) \right), \quad \text{where} \quad a_{NT} \equiv \frac{\sqrt{\tau NT}}{\beta(m \tau) - \beta(\tau)}. $$
By equation (ref), the transformation $\widehat{Z}_{NT}$ of $\widehat{\beta}(\tau)$ minimizes the objective $Q_{NT}(\cdot,\tau)$, where
To rewrite (ref) in a form more convenient for asymptotic analysis, we use Knight's identity: $$ \rho_\tau(u - v) - \rho_\tau(u) = -v(\tau - \mathbbm{1}\{u < 0\}) + \int_0^v \left( \mathbbm{1}\{u \le s\} - \mathbbm{1}\{u \le 0\} \right) ds. $$ Applying this identity, we can express the objective (ref) as
where $W_{NT}(\tau)$ is defined in (ref) and $G_{NT}(z, \tau)$ is given by $$ G_{NT}(z, \tau) \equiv \frac{a_{NT}}{\sigma^2_{NT} \sqrt{\tau NT}} \sum_{i=1}^N \sum_{t=1}^T \int_0^{\sigma_{NT} \cdot z / a_{NT}} \left[ \mathbbm{1}\left(Y_{it} - \beta(\tau) \le s\right) - \mathbbm{1}\left(Y_{it} - \beta(\tau) \le 0\right) \right] ds. $$
To analyze the asymptotic behavior of $\widehat{Z}_{NT}$, therefore, we first study the asymptotic properties of the components, $W_{NT}(\tau)$ and $G_{NT}(z, \tau)$, in the decomposition of the objective function $Q_{NT}(z, \tau)$ in (ref). To this end, we introduce a set of assumptions.
Recall that $F_Y$ denotes the distribution function of $Y$. Assume, without loss of generality, by location normalization that the lower endpoint of its support is $s_Y = 0$ or $s_Y = -\infty$ for notational convenience. Suppose that $F_Y$ satisfies a domain-of-attraction condition of Type 1, 2, or 3, defined as follows.
where $a(y) \equiv \int_{s_Y}^y F_Y(v) d v / F_Y(y) \text { for } y>s_Y$. This requirement is mild, as most common distribution families fall into one of the three types: Type 1 (e.g., normal), Type 2 (e.g., Student's $t$), or Type 3 (e.g., uniform). The parameter $\xi$ is referred to as the tail index, which characterizes the heaviness of the left tail. Specifically, the tail is heavy if $\xi > 0$, light (or thin) if $\xi = 0$, and bounded if $\xi < 0$.
In addition, we make the following assumptions.
Assumption (ref) requires the existence and regular variation of the unconditional and conditional quantile density functions. The exponents may differ across $f_Y\left(F_Y^{-1}(\tau) \mid \alpha_i\right)$, $f_Y\left(F_Y^{-1}(\tau) \mid \gamma_t\right)$, and $f_Y\left(F_Y^{-1}(\tau) \mid \alpha_i, \gamma_t\right)$. Assumption (ref) places constraints on the conditional cumulative distribution function of $Y$. For instance, part (i) requires that, as $\tau \to 0$, the second moment of the conditional probability $P\left(Y_{it} < \beta(\tau) \mid \alpha_i\right)$ diminishes at a rate faster than or comparable to $\tau$ (or $\tau^2 N$ if $N$ is small). In other words, the conditional probability must decay sufficiently fast as $\tau$ becomes small. The second moment of the unconditional probability $P(Y < \beta(\tau))$ scales with $\tau^2$, so the condition imposed by Assumption (ref) is relatively mild and, in particular, weaker than requiring the conditional probability to decay at the same rate as the unconditional one. Assumption (ref) imposes regularity conditions on the relationship between the conditional and unconditional density functions of $Y$. For example, part (i) asserts that the conditional density $f_Y(y \mid \alpha_i)$ should not vanish at a rate faster than the unconditional density $f_Y(y)$ as $y \to s_Y$. In essence, the decline of the conditional density in the lower tail should be comparable to that of the unconditional density. Assumption (ref) imposes moment conditions necessary to satisfy the Lyapunov condition, ensuring asymptotic Gaussianity of the normalized sample quantile.
Under these four assumptions, we can establish
Formal derivations are provided in Lemmas (ref) and (ref) in Appendices (ref) and (ref), respectively. As a consequence, we obtain the following limiting distribution for the solution $\widehat{Z}_{NT}$ to $\min_{z \in \mathbb{R}} Q_{NT}(z,\tau)$, where $Q_{NT}(z, \tau) = \sigma_{NT}^{-1} W_{NT}(\tau) \cdot z + G_{NT}(z, \tau)$ as defined in equation (ref).
While the formal proof is deferred to Appendix (ref), we briefly highlight its key components here because it is of interest on its own. The statistic $W_{NT}(\tau)$ is decomposed as $$ W_{NT}(\tau) = A_{1,\tau} + A_{2,\tau} + A_{3,\tau} + A_{4,\tau}, $$ where $A_{1,\tau}$, $A_{2,\tau}$, and $A_{4,\tau}$ are the projections onto $\alpha_i$, $\gamma_t$, and $(\alpha_i, \gamma_t)$, respectively, and $A_{3,\tau}$ denotes the remainder term. Under two-way clustered sampling, the remainder term $A_{3,\tau}$ is, in general, asymptotically non-Gaussian. However, by considering an appropriate sequence of quantiles with $\tau \to 0$, we can ensure that $A_{3,\tau}$ is asymptotically negligible relative to the leading terms, which are Gaussian.\footnote{Assumption (ref) ensures that $\mathrm{Var}(A_{3,\tau}) = o(\tau)$ and $\mathrm{Var}(A_{4,\tau}) = 1$, which together imply that the term $A_{3,\tau}$ is asymptotically negligible relative to the Gaussian component $A_{4,\tau}$.} (See the proof of Lemma (ref) in Appendix (ref) for further details.) As a result, we obtain the asymptotic Gaussianity stated in equation (ref), despite the potential complications introduced by two-way clustering. This robust Gaussianity is a special feature of extreme quantiles, not shared by intermediate quantiles.
The corresponding convergence rate in cross-sectional settings is $\sqrt{\tau N}$ Chernozhukov2005. By analogy, one might expect a convergence rate of $\sqrt{\min\{\tau N, \tau T\}}$ under two-way clustered sampling. However, the actual convergence rate turns out to be faster than this benchmark. The intuition behind this apparently counterintuitive result is as follows. The convergence rate of the estimator $\widehat{Z}_{NT}$ is given by $\displaystyle\frac{a_{NT}}{\sigma_{NT}}$. In the non-degenerate case (i.e., when $\mathrm{Var}(W_{NT} \mid \alpha_i) > 0$ or $\mathrm{Var}(W_{NT} \mid \gamma_t) > 0$), we have $ \sigma_{NT} = O\left( \sqrt{\mathrm{Var}(A_{1,\tau})} \vee \sqrt{\mathrm{Var}(A_{2,\tau})} \right). $ By Assumption (ref), we obtain $\sqrt{\mathrm{Var}(A_{1,\tau})} = o(T)$ and $\sqrt{\mathrm{Var}(A_{2,\tau})} = o(N)$, which together imply $ \sigma_{NT} = o\left( \sqrt{N} \vee \sqrt{T} \right). $ As a result, we obtain $ \frac{a_{NT}}{\sigma_{NT}} = o\left( \sqrt{\tau N} \vee \sqrt{\tau T} \right), $ indicating that the rate of $\widehat{Z}_{NT}$ is faster than $\sqrt{\min\{\tau N, \tau T\}}$. In the degenerate case (i.e., $\mathrm{Var}(W_{NT} \mid \alpha_i) = 0$ and $\mathrm{Var}(W_{NT} \mid \gamma_t) = 0$), we have $\sigma_{NT} = O\left( \sqrt{\mathrm{Var}(A_{4,\tau})} \right) = O(1)$, and hence the convergence rate simplifies to $a_{NT} = \sqrt{\tau NT}$. This corresponds to the rate obtained in Chernozhukov2005 under $NT$ i.i.d.\ observations.
As implied by the discussion after Theorem (ref) from the previous subsection, we have the decomposition \[ \sigma^2_{NT} = \mathrm{Var}\left(W_{NT}(\tau)\right) = \mathrm{Var}(A_{1,\tau} + A_{2,\tau} + A_{3,\tau} + A_{4,\tau}), \] where the components are
By construction, the four components are orthogonal, and we have $$ \mathrm{Var}(A_{1,\tau} + A_{2,\tau} + A_{3,\tau} + A_{4,\tau}) = \mathrm{Var}(A_{1,\tau}) + \mathrm{Var}(A_{2,\tau}) + \mathrm{Var}(A_{3,\tau}) + \mathrm{Var}(A_{4,\tau}). $$
From the previous subsection, we know that the asymptotic variance of $A_{3,\tau} + A_{4,\tau}$ converges to 1. On the other hand, the variances of $A_{1,\tau}$ and $A_{2,\tau}$ are given by
where $D_{it} = \tau - \mathbbm{1}\left\{ Y_{it} < \beta(\tau) \right\}$. Thus, the asymptotic variance of $W_{NT}(\tau)$ depends on the behavior of $\tau$, $N$, $T$, and the second moments of the conditional probabilities.
For example, if $\tau T = o(1)$ and $\mathbb{E} \left[ P\left(Y_{it} < \beta(\tau) \mid \alpha_i \right)^2 \right] = o(\tau^2)$, then $\mathrm{Var}(A_{1,\tau}) \to 0$. Similarly, if $\tau N = o(1)$ and $\mathbb{E} \left[ P\left(Y_{it} < \beta(\tau) \mid \gamma_t \right)^2 \right] = o(\tau^2)$, then $\mathrm{Var}(A_{2,\tau}) \to 0$. In this case, the entire asymptotic variance of $W_{NT}(\tau)$ is driven by $A_{3,\tau} + A_{4,\tau}$ and converges to 1. Since the rates of these components are unknown in practice, we propose using the sample counterpart of $\sigma_{NT}$ to construct confidence intervals, rather than relying on its asymptotic approximation.
Note that $ \left| \mathbb{E}\left[ \tau - \mathbbm{1}\left\{ Y_{it} < \beta(\tau) \right\} \mid \alpha_i \right] - \mathbb{E}\left[ \tau - \mathbbm{1}\left\{ Y_{it} < \widehat{\beta}(\tau) \right\} \mid \alpha_i \right] \right| = o_p(1) $ uniformly over $i$, due to the smoothness of the conditional CDF and the consistency of $\widehat{\beta}(\tau)$. We are going to claim that $ \frac{1}{T} \sum_{t=1}^T \left( \tau - \mathbbm{1}\left\{ Y_{it} < \widehat{\beta}(\tau) \right\} \right) $ is a consistent estimator of $\mathbb{E}\left[ \tau - \mathbbm{1}\left\{ Y_{it} < \widehat{\beta}(\tau) \right\} \mid \alpha_i \right]$. To see this, observe that, for each fixed $i$, the observations $(Y_{it})_{t=1}^T$ are conditionally independent given $\alpha_i$. Therefore, for each fixed $i$ and any $y \in \mathbb{R}$, $ \frac{1}{T} \sum_{t=1}^T \mathbbm{1}\left\{ Y_{it} < y \right\} \stackrel{p}{\to} \mathbb{E}\left[ \mathbbm{1}\left\{ Y_{it} < y \right\} \mid \alpha_i \right] $ uniformly over $y$, by the uniform law of large numbers for empirical CDFs. Following the same arguments as in Equation (E.1) in the proof of Theorem 5 in chiang2024standard, this convergence holds uniformly over $(y, i) \in \mathbb{R} \times \mathbb{N}$.
Analogously, the conditional expectation $\mathbb{E} \left[ \tau - \mathbbm{1}\left\{ Y_{it} < \beta(\tau) \right\} \mid \gamma_t \right]$ can be consistently estimated by $ \frac{1}{N} \sum_{i=1}^N \left( \tau - \mathbbm{1}\left\{ Y_{it} < \widehat{\beta}(\tau) \right\} \right). $
Although the conditional expectation $\mathbb{E} \left[ \tau - \mathbbm{1}\left\{ Y_{it} < \beta(\tau) \right\} \mid \alpha_i, \gamma_t \right]$ cannot be directly estimated consistently, it is not necessary to estimate it. Under Assumptions (ref)--(ref), we have $\mathrm{Var}(A_{3,\tau} + A_{4,\tau}) \to 1$.
Putting all these pieces together, we define the variance estimator by \[ \widehat{\sigma}^2_{NT} = \frac{1}{\tau NT} \sum_{i=1}^N \left( \sum_{t=1}^T \left( \tau - \mathbbm{1}\left\{ Y_{it} < \widehat{\beta}(\tau) \right\} \right) \right)^2 + \frac{1}{\tau NT} \sum_{t=1}^T \left( \sum_{i=1}^N \left( \tau - \mathbbm{1}\left\{ Y_{it} < \widehat{\beta}(\tau) \right\} \right) \right)^2 + 1, \] and are going to show that it satisfies the consistency property: \[ \frac{\widehat{\sigma}_{NT}}{\sigma_{NT}} \stackrel{p}{\to} 1. \]
We now turn to the estimation of the tail index $\xi$ that appears in equation (ref). Define \[
\] Of these, the convergence $\varphi \xrightarrow{p} 1$ follows from the decomposition: \[
\] since the first two terms on the right-hand side are $o_p(\sqrt{\tau N} \vee \sqrt{\tau T})$ by Theorem 1.
Given that $\varphi \xrightarrow{p} 1$, it follows that \[ - \frac{1}{\log l} \log \widehat{\rho}_l \to \xi, \] by an argument analogous to the proof of Theorem 6.1, using properties (v) and (vi) from Lemma 9.1 in Chernozhukov2005.\footnote{ The only part of Theorem 6.1 in Chernozhukov2005 where the dependence structure plays a role is in establishing the convergence $\varphi \to 1$. For all other parts, the presence of two-way clustering does not affect the argument. Therefore, once $\varphi \xrightarrow{p} 1$ is established under two-way clustering, it follows that $ - \frac{1}{\log \ell} \log \widehat{\rho}_{\ell} \to \xi. $ }
Consequently, we have $ \widehat{\xi} \equiv - \frac{\log \widehat{\rho}_l}{\log l} \xrightarrow{p} \xi. $ Also, if we define $ \widehat{a}_{NT} \equiv \frac{\sqrt{\tau NT}}{\widehat{\beta}(m \tau) - \widehat{\beta}(\tau)}, $ then a consistent estimator of the variance of $\widehat{\beta}(\tau)$ is given by $$ \frac{\widehat{\xi}^2}{\left(m^{-\widehat{\xi}} - 1\right)^2} \cdot \frac{\widehat{\sigma}^2_{NT}}{\widehat{a}_{NT}^2}. $$ We summarize the above result in the following proposition.
This section presents simulation studies to evaluate the finite-sample performance of the inference procedure based on the asymptotic Gaussianity results. Two designs are employed.
First, consider the following data-generating process: $$ Y_{it} = \sigma_\alpha \alpha_i + \sigma_\gamma \gamma_t + \sigma_\epsilon \epsilon_{it}, $$ where $\alpha_i \sim \mathcal{N}(0,1)$, $\gamma_t \sim \mathcal{N}(0,1)$, and $\epsilon_{it} \sim \mathcal{N}(0,1)$ are mutually independent. Throughout the simulations, we fix $\sigma_\epsilon = 2.0$, and vary $\sigma_\alpha \in \{1.0, 1.5, 2.0, 2.5, 3.0\}$ and $\sigma_\gamma \in \{1.0, 1.5, 2.0, 2.5, 3.0\}$. The sample size is set to $N = 200$ and $T = 200$, and we consider two extreme quantile levels: $\tau \in \{0.05, 0.01\}$ that are relevant to major empirical applications.
For each Monte Carlo iteration, we compute the point estimate and the corresponding 95% confidence interval based on the Gaussian approximation as described in Section (ref) and the variance estimator introduced in Section (ref).
Table (ref) reports the coverage frequencies of the 95% confidence intervals over 1{,}000 Monte Carlo iterations. The left panel reports the results for $\tau = 0.05$, while the right panel reports those for $\tau = 0.01$. Within each panel, the rows correspond to different values of $\sigma_\alpha \in \{1.0, 1.5, 2.0, 2.5, 3.0\}$, while the columns correspond to different values of $\sigma_\gamma \in \{1.0, 1.5, 2.0, 2.5, 3.0\}$. The results show that the coverage frequencies are fairly close to the nominal level of 95%, thereby supporting the validity of the proposed inference method based on asymptotic Gaussianity.
We next consider an alternative data-generating process: $$ Y_{it} = \sigma_\alpha \alpha_i \cdot \sigma_\gamma \gamma_t + \sigma_\epsilon \epsilon_{it}, $$ where $\alpha_i \sim \mathcal{N}(0,1)$, $\gamma_t \sim \mathcal{N}(0,1)$, and $\epsilon_{it} \sim \mathcal{N}(0,1)$ are mutually independent. In this design, the model incorporates a non-additive interaction between the latent factors $\alpha_i$ and $\gamma_t$. All other aspects of the simulation design remain the same as in the previous setup.
For each Monte Carlo iteration, we compute the point estimate and the corresponding 95% confidence interval using the asymptotic Gaussian approximation (Section (ref)) and the variance estimator described in Section (ref).
Table (ref) reports the coverage frequencies of the 95% confidence intervals across 1{,}000 Monte Carlo iterations. The left panel reports the results for $\tau = 0.05$, while the right panel reports those for $\tau = 0.01$. Within each panel, the rows correspond to different values of $\sigma_\alpha \in \{1.0, 1.5, 2.0, 2.5, 3.0\}$, while the columns correspond to different values of $\sigma_\gamma \in \{1.0, 1.5, 2.0, 2.5, 3.0\}$. Again, the coverage frequencies are close to the nominal level of 95%. These results further support the robustness of our inference procedure based on asymptotic Gaussianity.
The baseline analysis thus far focuses on unconditional distributions. In this section, we extend it to the setting of quantile regression. Before proceeding, we reset and redefine the variables and notations. Let $Y$ denote a response variable taking values in $\mathbb{R}$, and let $X = \left(1, X_{-1}^\prime\right)^\prime$ be a $d \times 1$ vector of regressors. Here, $X_{-1}$ denotes the vector $X$ with its first coordinate (corresponding to the constant term) removed. Denote the conditional distribution function of $Y$ given $X = x$ by $F_Y(\cdot \mid x)$.
Suppose that we observe data $\{(Y_{it}, X_{it}^\prime) : 1 \le i \le N,\ 1 \le t \le T\}$ generated by $ (Y_{it}, X_{it}^\prime) = g(\alpha_i, \gamma_t, \varepsilon_{it}) $ for some Borel measurable function $g$, where $\alpha_1, \ldots, \alpha_N$, $\gamma_1, \ldots, \gamma_T$, and $\varepsilon_{11}, \ldots, \varepsilon_{NT}$ are identically distributed and mutually independent random variables. As in the unconditional quantile case, this formulates two-way cluster dependence through the latent factors $\alpha_i$ and $\gamma_t$.
Our objective is to conduct inference on features of the conditional quantile function \[ F_Y^{-1}(\tau \mid x) = \inf \left\{ y \in \mathbb{R} : F_Y(y \mid x) > \tau \right\}, \] where $\tau$ is close to zero. Assume that the conditional quantile function takes the linear form:
With this formulation, the coefficient vector $\beta(\tau)$ is the parameter of interest, and it is estimated by
where $\rho_tau$ defined by $\rho_\tau(u) = u(\tau - \mathbbm{1}\{u < 0\})$ denotes the check function.
As in the unconditional quantile case, we are interested in the limiting distribution of the statistic \[ \widehat{Z}_{NT} = r_{NT} \left( \widehat{\beta}(\tau) - \beta(\tau) \right), \] where $\{r_{NT}\}$ is a sequence that governs the rate of convergence. Our focus is on the asymptotic behavior of $\widehat{Z}_{NT}$, and thus of $\widehat{\beta}(\tau)$, in the left tail of the conditional distribution of $Y$ given $X$, in the limit as $\tau \to 0$.
For a random variable $u$ with its distribution function $F_u$, let the lower bound of its support be denoted by $s_u = 0$ or $-\infty$, for notational convenience and without loss of generality by location normalization. We reintroduce the domain-of-attraction condition in three types:
where $a(z) \equiv \int_{s_u}^y F_u(z) d v / F_u(z) \text { for } z>s_u$ With this definition in place, we now state the three assumptions inherited from Chernozhukov2005.
Assumptions (ref)--(ref) correspond to Conditions R1--R3 in Chernozhukov2005, respectively. Under these assumptions, Theorem 3.1 of Chernozhukov2005 implies that \[ K(x) =
\] for some $\mathbf{c} \in \mathbb{R}^d$.
With the above tools, we now turn to the asymptotic theory in our framework. Define $$ W_{NT}(\tau) \equiv \frac{-1}{\sqrt{\tau NT}} \sum_{i=1}^N \sum_{t=1}^T \left( \tau - \mathbbm{1}\left\{ Y_{it} < X_{it}^\prime \beta(\tau) \right\} \right) X_{it}, $$ and denote its variance by $$ \Sigma_{NT} = \mathrm{Var} \left( W_{NT}(\tau) \right). $$ Also, define the scaling factor $$ a_{NT} \equiv \frac{\sqrt{\tau NT}}{\mu_X^\prime \left( \beta(m \tau) - \beta(\tau) \right)}, \qquad \mu_X^\prime = \mathbb{E}[X^\prime], $$ and the normalized statistic $$ \widehat{Z}_{NT} \equiv a_{NT} \cdot \Sigma_{NT}^{-1/2} \cdot \mathcal{Q}_H \cdot \left( \widehat{\beta}(\tau) - \beta(\tau) \right), $$ where $$ \mathcal{Q}_H \equiv \mathbb{E} \left[ H(X)^{-1} X X^\prime \right] \quad with \quad H(x) \equiv
$$
To establish the asymptotic Gaussianity of $\widehat{Z}_{NT}$, we first consider the convergence of an intermediate statistic: $$ \widetilde{Z}_{NT} \equiv \frac{a_{NT}}{c_{NT}} \left( \widehat{\beta}(\tau) - \beta(\tau) \right), $$ where $c_{NT} \in \mathbb{R}_+$ is a scaling factor satisfying $$ \frac{c_{NT}}{a_{NT}} = o(1) \qquad \text{and} \qquad \frac{\Sigma_{NT}}{c_{NT}^2} = O(1). $$ Assumption (ref), introduced later, ensures the existence of such a scalar $c_{NT}$, as it guarantees that $\Sigma_{NT}$ converges faster than $\tau NT$. We will show that $\widetilde{Z}_{NT}$ is asymptotically Gaussian, and it then follows that $\widehat{Z}_{NT}$ also inherits asymptotic Gaussianity.
Similarly to the unconditional quantile case, it follows from equation (ref) that the transformation $\widetilde{Z}_{NT}$ of $\widehat{\beta}(\tau)$ minimizes the objective function
Using Knight's identity, the objective (ref) can be rewritten as
where \[ W_{NT}(\tau) \equiv \frac{-1}{\sqrt{\tau NT}} \sum_{i=1}^N \sum_{t=1}^T \left( \tau - \mathbbm{1}\left\{ Y_{it} < X_{it}^\prime \beta(\tau) \right\} \right) X_{it} \] and \[ G_{NT}(z, \tau) \equiv \frac{a_{NT}}{\sqrt{\tau NT}} \sum_{i=1}^N \sum_{t=1}^T \int_0^{X_{it}^\prime z / a_{NT}} \left( \mathbbm{1}\left\{ Y_{it} - X_{it}^\prime \beta(\tau) \le s \right\} - \mathbbm{1}\left\{ Y_{it} - X_{it}^\prime \beta(\tau) \le 0 \right\} \right) ds. \]
We have already stated Assumptions (ref)--(ref), which are inherited from Chernozhukov2005. In addition, we introduce a new assumption concerning the tail behavior of the conditional distributions given $\alpha_i$, $\gamma_t$, and both $\alpha_i$ and $\gamma_t$, which is specific to the two-way clustered data structure considered in our setting.
Assumption (ref) requires that Assumption (ref) continues to hold when conditioning on each of the variables $\alpha_i$, $\gamma_t$, and both $\alpha_i$ and $\gamma_t$, thereby accommodating the dependence structure inherent in two-way clustered data. The corresponding scaling functions $K_\alpha$, $K_\gamma$, and $K_{\alpha,\gamma}$ are allowed to differ from one another.
Finally, we introduce the following set of assumptions, which are analogous to those stated for the analysis of unconditional quantiles in Section (ref).
Assumptions (ref)--(ref) are analogous to Assumptions (ref)--(ref), respectively, in Section (ref) concerning unconditional quantiles. Assumption (ref) (i) requires that Assumption (ref) holds conditional on each of the latent variables $\alpha_i$, $\gamma_t$, and both $\alpha_i$ and $\gamma_t$. Assumption (ref) (ii) requires the existence of the quantile density function and its regular variation. It also imposes regular variation on the conditional quantile density functions. The exponents of the regular variation may differ across the conditional densities $f_u(F_u^{-1}(\tau) \mid \alpha_i)$, $f_u(F_u^{-1}(\tau) \mid \gamma_t)$, and $f_u(F_u^{-1}(\tau) \mid \alpha_i, \gamma_t)$. Similar to Assumption (ref), Assumption (ref) here imposes constraints on the conditional cumulative distribution function of the variable $u$ given $X$. Assumption (ref) governs the relationship between the unconditional and conditional density functions of $u$. For example, Assumption (ref) (i) requires that the density of $u$ does not decay significantly faster than its conditional counterpart given $\alpha_i$. In other words, the decline in the density of $u$ should not be much more rapid than the decline in the conditional density of $u$ given $\alpha_i$. Lastly, Assumption (ref) imposes certain moment conditions needed to verify the Lyapunov-type central limit theorem, which underpins the asymptotic Gaussianity of our estimator.
Under Assumptions (ref)--(ref), and analogously to the unconditional quantile case, we can show that $c_{NT}^{-1} W_{NT}(\tau)$ converges in distribution to a multivariate Gaussian distribution as $$ c_{NT}^{-1} W_{NT}(\tau) \xrightarrow{d} \mathcal{N}(0, \Gamma), $$ where $$ \Gamma = \lim_{N, T \to \infty,\ \tau \to 0} \frac{\Sigma_{NT}}{c_{NT}^2}. $$ In addition, we can establish the convergence of the second component: $$ G_{NT}(z, \tau) \xrightarrow{p} \mathbb{E}[G_{NT}(z, \tau)] = \frac{1}{2} \cdot \frac{m^{-\xi} - 1}{-\xi} \cdot z^{\prime} \mathcal{Q}_H z. $$ See Lemmas (ref) and (ref) in Appendices (ref) and (ref), respectively, for formal and detailed derivations of these results.
Consequently, we obtain the following limiting distribution for the solution $\widetilde{Z}_{NT}$ to the minimization problem $\min_{z \in \mathbb{R}^d} Q_{NT}(z, \tau)$, where $Q_{NT}(z, \tau) = c_{NT}^{-1} W_{NT}(\tau)^\prime z + G_{NT}(z, \tau)$, as defined in equation (ref). Since the limiting distribution of $\widetilde{Z}_{NT}$ is Gaussian, it follows immediately that $\widehat{Z}_{NT}$ also converges in distribution to a Gaussian distribution.
See Appendix (ref) for a detailed proof. Despite some additional complications, the key features of the argument remain essentially the same as those presented in the unconditional quantile case in Section (ref).
Similarly to the unconditional quantile case presented in Section (ref), we use the decomposition \[ \Sigma_{NT} = \mathrm{Var}\left(W_{NT}(\tau)\right) = \mathrm{Var}(A_{1,\tau} + A_{2,\tau} + A_{3,\tau} + A_{4,\tau}), \] where
By construction, the four components are orthogonal, and we have $$ \mathrm{Var}(A_{1,\tau} + A_{2,\tau} + A_{3,\tau} + A_{4,\tau}) = \mathrm{Var}(A_{1,\tau}) + \mathrm{Var}(A_{2,\tau}) + \mathrm{Var}(A_{3,\tau}) + \mathrm{Var}(A_{4,\tau}). $$ From the previous subsubsection, we know that the asymptotic variance of $A_{3,\tau} + A_{4,\tau}$ is $E[X_{it} X^\prime_{it}]$. On the other hand, the variances of $A_{1,\tau}$ and $A_{2,\tau}$ are given by
where $D_{it} = \tau - \mathbbm{1}\left\{ Y_{it} < X^\prime_{it}\beta(\tau) \right\}$. As in the case of unconditional quantiles, the asymptotic variance of $W_{NT}$ depends on the behavior of $\tau$, $N$, and $T$. Since the rates of these components are unknown in practice, we propose using the sample counterpart of $\Sigma_{NT}$ to construct confidence intervals, rather than relying on its asymptotic approximation.
Since $E[E[D_{it}|\alpha_i] E[D^\prime_{it}|\alpha_i]]$ and $E[E[D_{it}|\gamma_t] E[D^\prime_{it}|\gamma_t]]$ can be estimated with \\
respectivly, the variances of $A_{1,\tau}$ and $A_{2,\tau}$ can be estimated by
respectively. Although $E[E[D_{it}|\alpha_i, \gamma_t] E[D^\prime_{it}|\alpha_i, \gamma_t]]$ cannot be directly estimated consistently, it is not necessary to estimate. Therefore, we can construct $\widehat\Sigma_{NT}$ as
and show its consistency property \[ \Sigma_{NT}^{-1} \widehat\Sigma_{NT} \xrightarrow{p} 1. \]
To estimate the tail index $\xi$, for any $x,\dot{x} \in \mathcal{X}$, define \[ \varphi \equiv \frac{x^{\prime} \left( \widehat{\beta}(m \tau) - \widehat{\beta}(\tau) \right)}{x^{\prime} \left( \beta(m \tau) - \beta(\tau) \right)} \qquad\text{and}\qquad \widehat{\rho}_{x, \dot{x}, l} \equiv \frac{x^{\prime} \left( \widehat{\beta}(m l \tau) - \widehat{\beta}(l \tau) \right)}{\dot{x}^{\prime} \left( \widehat{\beta}(m \tau) - \widehat{\beta}(\tau) \right)}. \] As in the case of unconditional quantiles, we have $\varphi \xrightarrow{p} 1$ by Theorem 2. Consequently, \[ - \frac{1}{\log l} \log \widehat{\rho}_{\overline{X}, \overline{X}, l} \to \xi \] follows by applying the same argument as in the unconditional quantile case. Also, define \[ \widehat a_{NT} = \frac{\sqrt{\tau NT}}{(\widehat\beta(m \tau)-\widehat\beta(\tau))} \] and \[ \widehat{\mathcal{Q}}_H = \frac{1}{NT}\sum_{i=1}^N\sum_{t=1}^T \left(\widehat{\rho}^{-1}_{X_{it}, \overline X, 1} X_{it} X^\prime_{it}\right). \] Thus, we obtain the following result.
In this section, we demonstrate an application of our proposed method to explore the growth-at-risk (GaR) framework of \citet*{adrian2022term}, and to analyze the role of financial conditions in shaping the distribution of forecast GDP growth across 11 advanced economies (AEs). These 11 AEs are identified by the IMF as having systemically important financial sectors.
Following adrian2022term, we estimate the dynamics of the GDP growth distribution over horizons ranging from 1 to 12 quarters using local projections in a panel setting. Specifically, we examine the conditional distributions of GDP growth for both near-term (1--4 quarters ahead) and medium-term (5--12 quarters ahead) forecast horizons. In conducting the inference, we account for the two-way clustering structure inherent in the panel data, where two ways of strong dependence arise from country-specific factors $\alpha_i$ and time-specific factors $\gamma_t$.
Let $\Delta y_{i, t+h}$ denote the annualized average GDP growth rate for country $i$ from quarter $t$ to quarter $t+h$. Let $x_{i, t}$ denote a vector of conditioning variables, which includes current GDP growth, the inflation rate, and the financial conditions index adrian2022term.
Our estimation and inference procedure is based on Equation (1) in adrian2022term. Specifically, we estimate the conditional quantile of $\Delta y_{i, t+h}$ via panel quantile regression on $x_{i, t}$. The quantile-specific regression slope $\delta_\tau^{(h)}$ is estimated by \[ \widehat{\delta}_\tau^{(h)} = \rm argmin_{\delta \in \mathbb{R}^d} \sum_{i=1}^{N} \sum_{t=1}^{T-h} \rho_\tau \left( \Delta y_{i, t+h} - x_{i, t}^\prime \delta \right), \] where $\rho_\tau$ defined by $\rho_\tau(u) = u(\tau - \mathbbm{1}\{u < 0\})$ is the check function. For inference, we apply our Gaussian limit distribution together with our proposed two-way cluster-robust standard error estimator specifically developed for extremal quantiles.
Table (ref) presents the estimates of the coefficients for various horizons ($h=2$, $h=4$, $h=6$, $h=8$, and $h=10$) focusing on the 5% conditional quantile, together with their two-way cluster-robust standard errors. The table includes the coefficients for the FCI, current GDP growth, and inflation across the horizons $h \in \{2,4,8,10,12\}$ and quantiles $\tau \in \{0.04,0.02,0.01\}$.
Observe that the FCI has a significantly negative impact on GDP growth at near-term horizons ($h = 2$ and $h = 4$), while its effects turn positive at medium-term horizons ($h = 8$ and $h = 10$). This pattern suggests that looser financial conditions are associated with higher short-run growth but may be linked to increased downside risk over the medium term. We emphasize that these significant results are obtained even after accounting for the two-way clustering structure in the panel data as well as the extreme quantiles.
The coefficients on current GDP growth (dlgdp) are positive and statistically significant across all forecast horizons, indicating that recent economic momentum is a strong predictor of future growth. Inflation (infl) generally exerts a negative effect on GDP growth, although the magnitude of this impact tends to diminish at longer horizons.
These empirical results are broadly consistent with adrian2022term, and our robust inference results reconfirm the earlier findings. The value added by our analysis is thus the guaranteed robustness of these empirical conclusions with statistical inference accounting for the two-way clustering and extreme quantiles.
This study investigates extreme quantiles in the context of data sampled under two-way clustering. Our objective is to characterize the asymptotic properties of extreme quantiles in the presence of such dependence structures.
We begin by analyzing the behavior of unconditional quantiles in the tails when the data exhibit two-way clustering, and subsequently extend our analysis to the case of quantile regression. As noted by Menzel2021, the presence of two-way clustering can lead to asymptotic distributions of sample means and related statistics that deviate from Gaussianity. In this study, we demonstrate that extreme quantiles retain asymptotic Gaussianity even under two-way clustering for extremal quantiles of intermediate orders. This result enables valid and robust inference for tail behavior despite the complex dependence structure. Simulation studies corroborate our theoretical findings.
Accounting for both two-way clustering and tail estimation is essential in applied research. By incorporating these considerations, our framework provides a more robust and reliable basis for inference on the extremal behavior of distributions when data are subject to multi-dimensional clustering.
We apply our theoretical framework to the empirical analysis of growth-at-risk. While our contribution does not concern identification or estimation per se, but rather the construction of valid standard errors, it is nevertheless important for empirical research to assess whether existing findings are robust to appropriate inference procedures. We show that, even after accounting for two-way clustering and extreme quantiles, the main empirical conclusions in the existing literature remain qualitatively unchanged.