High-Dimensional Extreme Quantile Regression

Quantile regression has become a widely recognized and useful alternative to classical least-squares regression for analyzing heterogeneous data. Since its introduction by Koenker and Bassett (1978), quantile regression has gradually been extended to a wide variety of data analytic settings; for a comprehensive review, see Koenker (2005); Koenker et al. (2017).

While traditional quantile regression allows exploration of a wide range of conditional quantiles for $\tau\in[\tau_{l},\tau_{u}]$ with $\tau_{l},\tau_{u}\in(0,1)$, there is often interest in the extreme tails where $\tau$ is close to 0 or 1. Without loss of generality, we focus the discussion on $\tau$ close to 1. The inherent challenge in estimating tail quantiles lies in the fact that the number of observations in the tails, that is, above the $\tau$th quantile, is often small. In various contexts, methods have been developed to study extreme quantile regression by leveraging extreme value theory. Chernozhukov (2005), Wang, Li and He (2012), and Wang and Li (2013) studied extreme conditional quantiles under the linear quantile regression framework, and Li and Wang (2019) focused on the linear quantile autoregressive model. Daouia et al. (2013) and Gardes and Stupfler (2019) considered nonparametric smoothing methods. Wang and Tsai (2009) and Youngman (2019) applied generalized Pareto models to analyze the exceedance at a high threshold. Velthoen et al. (2023) and Gnecco et al. (2024) explored tree-based methods such as gradient boosting and random forest.

With advancements in data collection techniques such as genomics, economics, finance, and imaging studies, the dimension of covariates $p$ is becoming larger and can grow with the sample size $n$. Existing methodology and theory in quantile regression for high-dimensional covariates have primarily focused on a central quantile level or compact quantile sets in $(0,1)$. In situations where $p$ grows with $n$, various studies have investigated asymptotic behaviors of quantile regression estimators, e.g., Welsh (1989), He and Shao (2000), Belloni et al. (2019), Pan and Zhou (2021), and He et al. (2023). Assuming the sparsity in regression coefficients, researchers have explored the $\ell_{1}$-penalized quantile regression estimator and its generalizations (Belloni and Chernozhukov, 2011; Fan et al., 2014), concave penalties (Wang, Wu and Li, 2012), and smoothed quantile regression with concave penalties (Tan et al., 2022). Additionally, the debiased estimator of high-dimensional quantile regression has been studied for inference at central quantiles (Zhao et al., 2014; Bradic and Kolar, 2017).

In this article, we address the challenge of estimating extreme conditional quantiles of $Y\in\mathbb{R}$ given a set of predictors $\mathbf{X}\in\mathbb{R}^{p}$ in a high-dimensional context. The challenge lies in the scarcity of efficient samples at the tail and the sparsity of samples in high-dimensional space. As noted by Gnecco et al. (2024), a lack of samples exceeding the corresponding conditional $\tau$th quantile leads to empirical estimators with large bias and variance. The high dimensionality of the predictor space introduces additional bias, as noise covariates can obscure true signals.

Current research on extreme quantiles primarily focuses on estimation within fixed-dimensional settings and is either not applicable or less competitive in high-dimensional settings. Although tree-based methods (Velthoen et al., 2023; Gnecco et al., 2024) demonstrate effectiveness in relatively large and fixed covariate dimensions, addressing scenarios where the covariate dimension grows with the sample size, comparable to the sample size itself, remains inadequately explored. Recently, Sasaki et al. (2024) proposed a high-dimensional tail index regression model, assuming that the dimension grows with $n$ and is comparable to $n(1-\tau)$. However, this approach imposes constraints such as a Pareto tail and specific link functions between the extreme value index and covariates.

In this article, we propose novel estimators for extreme quantiles in high dimensions by integrating concepts from extreme value theory with regularized estimators. Building upon a linear conditional quantile model, we employ regularized quantile regression to estimate intermediate conditional quantiles and extrapolate them to extreme quantiles. Here, “intermediate quantile” refers to a quantile level $\tau_{n}$ that approaches one at a moderate rate such that $n(1-\tau_{n})\rightarrow\infty$, whereas “extreme quantile” denotes a quantile level $\tau_{n}$ that approaches one more rapidly.

The theoretical analysis of extreme quantiles in high dimensions is challenging. The limiting distribution of a quantile estimator could be intractable in high dimensions, even at central quantile levels. Existing work on central quantiles within $(0,1)$ often requires a common assumption, namely $f_{Y}(Q_{Y}(\tau|\mathbf{X})|\mathbf{X})>C>0$, where $f_{Y}(\cdot|\mathbf{X})$ is the conditional density of $Y$ given $\mathbf{X}$, and $Q_{Y}(\tau|\mathbf{X})$ is the $\tau$th conditional quantile of $Y$; see for instance Belloni and Chernozhukov (2011), and Wang, Wu and Li (2012). However, this assumption is violated in the tails as $\tau$ approaches $1$. It is even more challenging to achieve a uniform result of the regularized estimator over the entire tail, which is crucial for constructing an effective extreme value index estimator. This paper addresses the gap and establishes the uniform error rate for quantile regression estimator in high dimensions and at intermediate quantiles. This rate includes terms analogous to those in high-dimensional quantile regression at a central quantile level, augmented by an inflationary component attributed to escalating quantile levels or diminishing effective sample sizes. In addition, we propose a refined Hill approach for estimating the extreme value index, which is based on a fixed number of intermediate high-dimensional quantile estimates. This framework enables us to analyze the uniform behavior of regularized quantile regression in the tail region and assess the rate at which $f_{Y}(Q_{Y}(\tau|\mathbf{X})|\mathbf{X})$ tends to zero. The error rates of the proposed refined Hill and extreme quantile estimators offer theoretical insights into threshold selection in high-dimensional contexts.

The rest of the article is organized as follows. In Section 2, we present the proposed estimators for the extreme value index and extreme conditional quantiles, and derive theoretical results for the proposed estimators using techniques from both extreme value theory and high-dimensional statistics theory. We assess the finite sample performance of the proposed methods through a simulation study in Section 3 and the analysis of auto insurance data in Section 4. Technical details and additional information for Sections 3 and 4 are provided in the online Supplementary Material.

Suppose we observe a random sample $\{(\mathbf{X}_{i},Y_{i}),i=1,\ldots,n\}$ of the random vector $(\mathbf{X},Y)$, where $Y_{i}$ is the univariate response variable and $\mathbf{X}_{i}=(X_{i1},\ldots,X_{ip})^{T}$ is the $p$-dimensional centralized design vector. This article considers the high dimensional case: $p:=p(n)\rightarrow\infty$ as $n\rightarrow\infty$. Let $Q_{Y}(\tau|\mathbf{X})=\inf\{y:F_{Y}(Y|\mathbf{X})\geq\tau\}$ denote the $\tau$th conditional quantile of $Y$ given $\mathbf{X}$, where $F_{Y}(\cdot|\mathbf{X})$ is the cumulative distribution function (CDF) of $Y$ given $\mathbf{X}$. Denote $\mathcal{X}\subset\mathbb{R}^{p}$ as the support of $\mathbf{X}$.

Throughout the article, we assume that $F_{Y}(\cdot|\mathbf{X})$ is in the maximum domain of attraction of an extreme value distribution $G_{\gamma}(\cdot)$ with the extreme value index (EVI) $\gamma>0$, denoted by $F_{Y}(\cdot|\mathbf{X})\in D(G_{\gamma})$. That is, for a given random sample $Y_{1},\ldots,Y_{n}$ from $F_{Y}(\cdot|\mathbf{X})$, there exist constants $a_{n}>0$ and $b_{n}\in\mathbb{R}$ such that

as $n\rightarrow\infty$, for $1+\gamma y\geq 0$. In this paper, we assume $\gamma>0$, which means that $Y|\mathbf{X}$ has heavy-tailed distributions as commonly seen in many applications, such as stock market returns, insurance claims, earthquake magnitudes, river flows during floods.

The main objective of this paper is to estimate the conditional quantile $Q_{Y}(\tau_{n}|\mathbf{X})$ when $\tau_{n}\rightarrow 1$ and $p\rightarrow\infty$ as $n\rightarrow\infty$. We address two distinct regimes for the quantile level $\tau_{n}$ as it approaches 1: intermediate order quantile levels such that $n(1-\tau_{n})\rightarrow\infty$, and extreme order quantile levels with $n(1-\tau_{n})\rightarrow C$, where $C\geq 0$ is some constant.

We focus on the following tail linear quantile regression model:

where $\mathbf{Z}=(1,\mathbf{X}^{T})^{T}$, $\tau_{ln}\rightarrow 1$ as $n\rightarrow\infty$, and the quantile slope coefficients $\boldsymbol{\beta}(\tau)$ may vary across $\tau\in$ $[\tau_{ln},1)$. The linear quantile model in (1) is specifically assumed at the upper tail. Similar tail model specifications have also been considered in Wang, Li and He (2012), Li and Wang (2019), Xu et al. (2022). In high-dimensional settings, we extend this assumption and additionally impose a sparsity condition to ensure model identifiability. We assume there are $s$ tail-relevant variables, that is, $s:=\#\{j=1,2,\ldots,p:\beta_{j}\left(\tau\right)\neq 0,$ $\exists\tau\in[\tau_{ln},1)\}$, and $s=o(n)$. A tail-relevant variable may influence the tail quantiles of $Y|\mathbf{X}$, even if it does not affect the central quantiles. Without loss of generality, we assume that the first $s$ slope coefficients are nonzero, i.e., $\boldsymbol{\beta}(\tau)=(\beta_{0}(\tau),\beta_{1}(\tau),\ldots,\beta_{s}(\tau),\mathbf{0}_{p-s})^{T}$.

We propose a three-step procedure to estimate the extreme conditional quantile. In the first step, we obtain $\ell_{1}$-penalized quantile estimators at a sequence of intermediate quantile levels. For heavy-tailed distributions, the Hill estimator (Hill, 1975) is the most commonly used method to estimate the EVI. However, based on the average of log excesses, the Hill estimator requires analysis of the tail quantile process (de Haan and Ferreira, 2006), which poses both computational and theoretical challenges in high dimensions. To address these challenges, we develop a refined Hill estimator in the second step, which relies on quantile estimates at a fixed number of intermediate quantile levels that approach to one at the same rate. In the third step, we develop an extrapolation estimator for the extreme conditional quantile $Q_{Y}(\tau_{n}|\mathbf{X})$ by leveraging extreme value theory and results from the first two steps.

Before presenting the proposed estimators and their theoretical properties, we first introduce the notations used throughout the paper. For a sequence of random variables $E_{1},E_{2},\ldots,E_{n}$, we denote its order statistics as $E_{(1)}\leq E_{(2)}\leq\ldots\leq E_{(n)}$. Given a random sample $\mathbf{Z}_{1},\ldots,\mathbf{Z}_{n}$, let $\mathbb{G}_{n}(f)=\mathbb{G}_{n}\{f(\mathbf{Z}_{i})\}:=n^{-1/2}\sum_{i=1}^{n}[f(\mathbf{Z}_{i})-\mathrm{E}\{f(\mathbf{Z}_{i})\}]$ and $\mathbb{E}_{n}f=\mathbb{E}_{n}f\left(\mathbf{Z}_{i}\right):=n^{-1}\sum_{i=1}^{n}f(\mathbf{Z}_{i})$. We denote the $\ell_{2}$, $\ell_{1}$, $\ell_{\infty}$ and $\ell_{0}$ norms by $\|\cdot\|$, $\|\cdot\|_{1}$, $\|\cdot\|_{\infty}$ and $\|\cdot\|_{0}$, respectively. Let $\|\boldsymbol{\beta}\|_{1,n}=\sum_{j=1}^{p}\widehat{\sigma}_{j}|\beta_{j}|$ denote the weighted $\ell_{1}$-norm with $\widehat{\sigma}_{j}^{2}:=\mathbb{E}_{n}(X_{ij}^{2})$. Given a vector $\boldsymbol{\delta}\in\mathbb{R}^{p}$, and a set of indices $T\subset\{1,\ldots,p\}$, we denote by $\boldsymbol{\delta}_{T}$ the vector where $\delta_{Tj}=\delta_{j}$ if $j\in T,\delta_{Tj}=0$ if $j\notin T$. We use $a\lesssim b$ to denote $a=O(b)$, meaning $a\leq cb$ for some constant $c>0$ that does not depend on $n$; and $a\asymp b$ to denote $a=O(b)$ and $b=O(a)$. We use $a\lesssim_{P}b$ to denote $a=O_{P}(b)$. Additionally, we use $a\vee b=\max\{a,b\}$ and $a\wedge b=\min\{a,b\}$. For notational simplicity, let $F_{i}(\cdot)=F_{Y}(\cdot|\mathbf{X}_{i})$ be the conditional distribution function of $Y$ given $\mathbf{X}_{i}$, and denote $f_{i}(\cdot)=f_{Y}(\cdot|\mathbf{X}_{i})$. Finally, let $\mathbb{X}=(\mathbf{X}_{1},\ldots,\mathbf{X}_{n})^{\top},\mathbb{Z}=(\mathbf{Z}_{1},\ldots,\mathbf{Z}_{n})^{\top}$, and $\mathbf{Y}=(Y_{1},Y_{2},\ldots,Y_{n})^{T}$.

We conduct a simulation study to investigate the performance of the proposed method for estimating extreme quantiles in high-dimensional settings. The data are generated from

where $X_{ij}\sim U(0,1)$ for $j=1,2,\ldots,p$ and $\varepsilon_{i}$ are independent and identically distributed student t random variables, and $w$ is a constant controlling the degree of heteroscedasticity. Under this model, the $\tau$th conditional quantile of $Y$ is $Q_{Y}(\tau|\mathbf{X}_{i})=\alpha(\tau)+\mathbf{X}_{i}^{T}\boldsymbol{\beta}(\tau),$ where $\mathbf{X}_{i}=\left(X_{i1},X_{i2},\ldots,X_{ip}\right)^{T}$, $\alpha(\tau)=Q_{\varepsilon}(\tau)$, $\boldsymbol{\beta}(\tau)=\left(1+wQ_{\varepsilon}(\tau),1,0\ldots,0\right)^{T}$, and $Q_{\varepsilon}(\tau)$ is the $\tau$th quantile of $\varepsilon_{i}$. We consider six different cases. In Cases 1, 3, and 5, $\varepsilon_{i}\sim t(5)$ and in Cases 2, 4 and 6, $\varepsilon_{i}\sim t(2)$. Cases 1-2 have homoscedastic errors with $w=0$, while $w=0.5$ in Cases 3-4, and $w=0.9$ in Cases 5-6. Note that, for $t(v)$, the EVI is $\gamma=1/v$. We consider two settings for the sample size $n$ and dimension $p$: (1) $p=p(n)$ grows with $n$ but $p$ is smaller than $n$, with $p=\lceil n^{0.5}\rceil$ and $n=1000,5000$; (2) $p$ is larger than or equal to $n$, with $(n,p)=(1000,1000)$ and $(1000,1200)$. For each scenario, the simulation is repeated $500$ times.

For each simulation, we estimate $Q_{Y}(\tau_{n}|\mathbf{x})$ at $\tau_{n}=0.995,0.999$ using three methods: the extreme quantile regression method (EQR) without assuming common slope in Wang, Li and He (2012); the high-dimensional quantile regression (HQR) of Belloni and Chernozhukov (2011), and the proposed high-dimensional extreme quantile regression (HEQR). For EQR, the number of upper-order statistics $k$ is set to $[4.5n^{1/3}]$, as recommended in Wang, Li and He (2012).

For each simulation, we calculate integrated squared error (ISE), defined as

where $\mathbf{X}_{1}^{*},\ldots,\mathbf{X}_{L}^{*}$ are evaluation points of the covariates, and we define the mean integrated squared error (MISE) as the average ISE across $500$ simulations. In our simulation, we set $L=100$, and let $\mathbf{X}_{l}^{*}$ be random replicates of $\mathbf{X}$.