# Extreme Conditional Quantile Estimation in High Dimensions: A Comparative Study

Laurent Gardes, Alex Podgorny

<figure><img src="image_1.png" /></figure>

## To cite this version:

Laurent Gardes, Alex Podgorny. Extreme Conditional Quantile Estimation in High Dimensions: A Comparative Study. 2025. hal-05059524v2

HAL Id: hal-05059524

https://hal.science/hal-05059524v2

Preprint submitted on 6 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Extreme Conditional Quantile Estimation in High Dimensions: A Comparative Study

L. Gardes(1) and A. Podgorny(1)

(1) Université de Strasbourg & CNRS, IRMA
7 rue René Descartes, 67000 Strasbourg, France.

## Abstract

This paper addresses the problem of estimating extreme conditional quantiles in high-dimensional settings. We mainly focus on the case where the conditional distribution is heavy-tailed. We consider recent estimation procedures based on machine learning techniques, including neural networks and random forests, as well as dimension reduction approaches specifically designed for extreme values. A comprehensive simulation study evaluates their performance across various scenarios, investigating the influence of the covariate dimension, the second-order tail behavior, the complexity of the link between covariate and tail index, the correlation structure of the random covariate and the choice of the intermediate quantile level. Our findings provide practical guidelines for selecting appropriate methods and highlight the strengths and limitations of each approach.

## 1 Introduction

Let $Y$ be a real-valued response variable and $X$ its associated random covariate of support $\mathcal{X} \subset \mathbb{R}^p$. Given a sample $(X_1, Y_1), \dots, (X_n, Y_n)$ of $n$ independent replications of $(X, Y)$, we are interested in the estimation of the conditional quantile of level $\alpha \in (0, 1)$ at a point $x \in \mathbb{R}^p$ given by

$$Q(\alpha | x) := \inf \{t \in \mathbb{R}; \bar{F}(t | x) \le \alpha\},$$

where $\bar{F}(t | x) := \mathbb{P}(Y > t | X = x)$ is the conditional survival function of $Y$ given $X = x$.

Numerous methods exist for conditional quantile estimation, see for instance Koenker and Bassett [18], Li and Racine [20] and Yu and Moyeed [26], but most of them are not reliable in a high-dimensional setting i.e., for a large dimension $p$ of the covariate. For instance kernel regression (Berlinet et al. [8] and Samanta [23]) or nearest-neighbor methods (Bhattacharya and Gangopadhyay [6]) suffer severely from the curse of dimensionality (Bellman [7]), leading
to poor estimation performance and computational inefficiency as the number of covariates increases.

To tackle the challenges of high-dimensional quantile estimation, several strategies have been proposed. In this paper, we focus on two broad families of methods. The first includes machine learning approaches such as neural networks, random forests, gradient boosting or penalized regression models. These methods aim to estimate the mapping $x \mapsto Q(\alpha | x)$ directly, leveraging the structure of the data through flexible, data-driven models. The second family of methods relies on dimension reduction techniques, such as Sliced Inverse Regression (SIR) (Li [19]) or Sliced Average Variance Estimation (SAVE) (Cook and Weisberg [11]). These approaches seek to identify a q-dimensional linear subspace $S \subset \mathbb{R}^p$ with $q \ll p$ that retains all the information in $X$ relevant for predicting $Y$. In this framework, the conditional quantile $Q(\alpha | x)$ is assumed to depend only on a projection $B^\top x$, for some basis $B \in \mathbb{R}^{p \times q}$ of $S$, reducing the effective dimensionality of the problem and making the classical nonparametric techniques usable in the reduced space.

In the remainder of this paper, we shift our attention from classical conditional quantiles to the more delicate problem of estimating an extreme conditional quantile $Q(\beta_n | x)$, where the level $\beta_n$ converges to 0 as the sample size increases. These quantiles play a crucial role in various domains, including finance, environmental science, and risk management, and their estimation presents additional challenges due to the scarcity of observations in the tail of the distribution.

To ensure reliable estimation in such cases, we assume that for all $x \in \mathcal{X}$, the conditional distribution of $Y$ given $X = x$ belongs to a maximum domain of attraction providing a theoretical framework for extrapolating beyond observed data. Specifically, we consider the following model:

(M) For all $x \in \mathcal{X} \subset \mathbb{R}^p$, the cumulative distribution function $F(\cdot | x)$ of $Y$ given $X = x$ is continuous, strictly increasing and there exist functions $\gamma: \mathcal{X} \to \mathbb{R}$ and $a: [1, \infty) \times \mathcal{X} \to (0, \infty)$ such that, for any $t > 0$,

$$ \lim_{\alpha \to 0^+} \frac{Q(t\alpha | x) - Q(\alpha | x)}{a(\alpha^{-1} | x)} = \frac{t^{-\gamma(x)} - 1}{\gamma(x)} =: K_{\gamma(x)}(t^{-1}), \quad (1) $$

where the right-hand side is interpreted as $\log(t^{-1})$ when $\gamma(x) = 0$.

The function $\gamma$ is the conditional extreme value index function. It drives the decay of the conditional tail distribution. In order to estimate extreme conditional quantiles, we exploit the asymptotic behavior given by model (M). More specifically, for a probability level $\beta_n$ close to zero, we apply the limit in (1) with $t = \beta_n/\alpha_n$ to get

$$ Q(\beta_n | x) \approx Q(\alpha_n | x) + \frac{a(\alpha_n^{-1} | x)}{\gamma(x)} \left[ \left( \frac{\alpha_n}{\beta_n} \right)^{\gamma(x)} - 1 \right], \quad (2) $$
where $\alpha_n$ is a sufficiently 'large' level in order that the conditional quantile $Q(\alpha_n | x)$ can be estimated without using the convergence (1) in model (M). In particular, in Gnecco et al. [15], the authors assume that $\alpha_n > \beta_n$ with $\alpha_n \to 0$ and $n\alpha_n \to \infty$. Such a quantile is referred to as a classical conditional quantile. Naturally, the faster the convergence in (1), the better the approximation (2). The rate of convergence in (1) is usually controlled by the following condition.

(SOC) For all $x \in \mathcal{X}$, there exist functions $\gamma : \mathcal{X} \to \mathbb{R}$, $\rho : \mathcal{X} \to (-\infty, 0]$, $a : [1, \infty) \times \mathcal{X} \to (0, \infty)$ and $A : [1, \infty) \times \mathcal{X} \to (0, \infty)$ with $A(y | x) \to 0$ as $y \to \infty$ such that, for any $t > 0$,

$$ \lim_{\alpha \to 0^+} \frac{1}{A(\alpha^{-1} | x)} \left( \frac{Q(t\alpha | x) - Q(\alpha | x)}{a(\alpha^{-1} | x)} - K_{\gamma(x)}(t^{-1}) \right) = H_{\gamma(x), \rho(x)}(t^{-1}), \quad (3) $$

with

$$ H_{\gamma(x), \rho(x)}(y) := \frac{1}{\rho(x)} (K_{\gamma(x)+\rho(x)}(y) - K_{\gamma(x)}(y)). $$

This condition is referred to as the second-order condition (see, for example, de Haan and Ferreira [16, Section 2.3]). The function $|A(\cdot | x)|$ is regularly varying of index $\rho(x) \le 0$ that is to say that for all $t > 0$,

$$ \lim_{y \to \infty} \left| \frac{A(ty | x)}{A(y | x)} \right| = t^{\rho(x)}. $$

As a consequence, the so-called second-order parameter $\rho(x)$ drives the rate of convergence in (1); a value of $\rho(x)$ close to 0 is associated to a slow convergence and thus to a poor approximation (2).

To sum up, estimating $Q(\beta_n | x)$ requires:

* the estimation of a classical conditional quantile $Q(\alpha_n | x)$;

* the estimation of the scale parameter $a_n(x) := a(\alpha_n^{-1} | x)$ and of the extreme value index $\gamma(x)$. These parameters govern the tail behavior and enable the extrapolation beyond classical conditional quantiles.

In the specific case where $\gamma(x) > 0$ for all $x \in \mathcal{X}$, we have that $a(\alpha_n^{-1} | x) \sim \gamma(x)Q(\alpha | x)$ as $\alpha \to 0$ (see for instance de Haan and Ferreira [16, Lemma 1.2.9]). In this case, convergence (1) rewrites

$$ \lim_{\alpha \to 0^+} \frac{Q(t\alpha | x)}{Q(\alpha | x)} = t^{-\gamma(x)}. \quad (4) $$

When $\gamma(x) > 0$, the conditional distribution is said to be heavy-tailed and $\gamma$ is referred to as the conditional tail index. This is equivalent to say that the survival function $\bar{F}(\cdot | x)$ is regularly varying of index $-1/\gamma(x)$ that is to say, for all $t > 0$,

$$ \lim_{y \to \infty} \frac{\bar{F}(ty | x)}{\bar{F}(y | x)} = t^{-1/\gamma(x)}. $$
Equation (4) leads to a simplified estimation strategy for the estimation of the extreme conditional quantile since only the estimation of the conditional tail index $\gamma$ is required. The approximated extreme quantile is then given by:

$$Q(\beta_n | x) \approx Q(\alpha_n | x) \left(\frac{\alpha_n}{\beta_n}\right)^{\gamma(x)}. \qquad (5)$$

We pay particular attention to the estimation of extreme conditional quantiles in the heavy-tailed case $\gamma(x) > 0$, as it is the most relevant and challenging setting in many real-world applications. Heavy-tailed distributions frequently arise in finance, insurance, environmental sciences, and various risk assessment problems, where accurate tail estimation is crucial.

In this paper, we focus on four recently proposed methods adapted to this setting. The first two are machine learning-based: the extremal random forest method developed by Gnecco et al. [15] and the extremal neural networks approach introduced by Pasche and Engelke [21]. These methods can deal with the general case $\gamma(x) \in \mathbb{R}$ in model (M). In Gnecco et al. [15], the authors propose a variant that is specifically dedicated to heavy-tailed distributions ($\gamma(x) > 0$). In this paper, we also propose an adaption of the extremal neural networks approach of Pasche and Engelke [21] to the case $\gamma(x) > 0$. The other two methods are based on dimension reduction techniques specifically designed for extremes: the TIREX method of Aghbalou et al. [2], and the CTI method proposed by Gardes and Podgorny [13]. As these methods have only been introduced recently, few numerical studies have evaluated their relative performance. In particular, the use of dimension reduction techniques developed specifically for extreme value analysis has not yet been systematically assessed for conditional extreme quantile estimation. The main contribution of this paper is to conduct a comprehensive numerical study, evaluating the performance of these techniques across a variety of models and settings. This analysis provides new insights into the strengths and limitations of each method, assessing their robustness to key factors such as the dimension of the covariate space, the choice of the classical quantile level $\alpha_n$, and the underlying model structure.

The structure of the paper is as follows. Section 2 is dedicated to the estimation of classical conditional quantiles using the generalized random forest algorithm. As explained above, this estimation is the first step in the estimation process of an extreme conditional quantile. In Section 3, the following methods for the estimation of an extreme quantile are presented:

1. extremal random forest (Gnecco et al. [15]) and its variant for heavy-tailed distribution,

2. extremal neural networks (Pasche and Engelke [21]) and its adaptation to heavy-tailed distribution,

3. dimension reduction with TIREX (Aghbalou et al. [2]),
4. dimension reduction with the CTI method (Gardes and Podgorny [13]).

Finally, these methods are compared through a simulation study in Section 4. Of course, many other methods can be considered. For instance, a gradient boosting approach has been proposed by Velthoen et al. [24] and Youngman [25]. To avoid overloading this paper, we chose not to include this approach in our study; however, a comparison with neural networks and random forests can be found in Pasche and Engelke [21] and Gnecco et al. [15].

## 2 Classical conditional quantile estimation

As mentioned in the introduction, the first step in the estimation of an extreme conditional quantile is the estimation of a classical conditional quantile $Q(\alpha_n | x)$ where the order $\alpha_n$ converges to 0 not too fast. Estimation of classical conditional quantiles has been the subject of extensive research. Various effective methods exist, but in this work, we systematically estimate it using the generalized random forest method of Athey et al. [5], which has demonstrated strong performance for nonparametric conditional quantile estimation in high-dimensional settings.

The random forest algorithm was initially introduced by Breiman [9] for estimating the conditional mean $\mu(x) := \mathbb{E}(Y | X = x)$. This method was later extended to the use of other loss functions by the generalized random forest algorithm, see Athey et al. [5]. As a consequence, this method can be used to estimate various quantities beyond the conditional mean, such as conditional quantiles. Hereafter, we present the generalized random forest algorithm for the estimation of a conditional quantile.

For a level $\alpha > 0$, the quantile $Q(\alpha | x)$ minimizes the conditional expectation of the quantile loss function $\rho_{\alpha}(c) := c(\mathbb{I}\{c > 0\} - \alpha)$:

$$Q(\alpha | x) = \arg \min_{q \in \mathcal{Y}} \mathbb{E}(\rho_{\alpha}(Y - q) | X = x),$$

where $\mathcal{Y}$ is the support of the distribution of the response variable $Y$. A natural way to estimate $Q(\alpha | x)$ is to replace the conditional expectation by its empirical counterpart:

$$Q(\alpha | x) := \arg \min_{q \in \mathcal{Y}} \sum_{i=1}^{n} W(x, X_i) \rho_{\alpha}(Y_i - q),$$

where $W(x, X_i)$ are positive weights summing to 1 and measuring the similarity between $Q(\alpha | x)$ and $Q(\alpha | X_i)$. The value of $W(x, X_i)$ is expected to be close to 0 if $Q(\alpha | x)$ and $Q(\alpha | X_i)$ are very different and close to 1 if $Q(\alpha | x) \approx Q(\alpha | X_i)$. One method among many others to estimate these similarity weights is the generalized random forest algorithm. More specifically, as for the classical random forest algorithm, generalized random forest aggregate B decision trees where each tree is drawn with replacement from the
original dataset. The algorithm is described below. Let $s_n(\alpha) \in \{1, \dots, n\}$ and $m_{RF}(\alpha) \in \{1, \dots, s_n(\alpha)\}$.

1. For $b \in \{1, \dots, B\}$, select a random subsample with replacement of size $s_n(\alpha)$ of the training set. Next, randomly select $k := \min(p, \lfloor \sqrt{p} \rfloor + 20)$ components of the covariate. We denote by $P^{(1,b)} := \{X_i^{(b)}; i = 1, \dots, s_n(\alpha)\}$ the obtained random vectors of dimension $k$ and by $\{Y_i^{(b)}; i = 1, \dots, s_n(\alpha)\}$ the associated responses variables.

2. For a finite set $E$, let $|E|$ be its number of elements. For $j = 1, 2, \dots$, as long as $|P^{(j,b)}| \ge m_{RF}(\alpha)$, split the node $P^{(j,b)}$ into two axis-aligned child nodes $C_1$ and $C_2 = P^{(j,b)} \setminus C_1$ and compute the criterion

$$ \Delta^{(j,b)}(C_1, C_2) = \sum_{\ell=1}^{2} \frac{1}{|C_{\ell}|} \left( \sum_{X_i^{(b)} \in C_{\ell}} \mathbb{I}\{Y_i^{(b)} > \hat{Q}(\alpha | P^{(j,b)})\} \right)^2, \quad (6) $$

and

$$ \hat{Q}(\alpha | P^{(j,b)}) := \arg \min_{q \in \mathcal{Y}} \sum_{i=1}^{n} \rho_{\alpha}(Y_i^{(b)} - q) \mathbb{I}\{X_i^{(b)} \in P^{(j,b)}\}. \quad (7) $$

3. Denote by $(P^{(2j,b)}, P^{(2j+1,b)})$ the couple which maximizes (6).

4. Repeat steps 2. and 3. until it is no longer possible. The obtained nodes are called the final leaves.

Let us emphasize that the couple that maximizes the criterion $\Delta^{(j,b)}(\cdot, \cdot)$ over all possible splits is the couple of sets where one contains all the $Y_i^{(b)}$'s exceeding the estimated conditional quantile $\hat{Q}(\alpha | P^{(j,b)})$ and the other all the $Y_i^{(b)}$'s located below $\hat{Q}(\alpha | P^{(j,b)})$. As a consequence, for observations of the covariate belonging to the same final leaf, the corresponding conditional quantiles are expected to be close.

The last step of the generalized random forest algorithm is the computation of the estimated similarity weights. More specifically, let $L_b(x)$ denote the terminal leaf of tree b containing x and let

$$ w^{(b)}(x, X_i) := \begin{cases} 1/|L_b(x)| & \text{if } X_i \in L_b(x), \\ 0 & \text{otherwise.} \end{cases} $$

The estimated similarity weights are given by

$$ W_{\alpha}^{(RF)}(x, X_i) := \frac{1}{B} \sum_{b=1}^{B} w^{(b)}(x, X_i). $$
Averaging over the B trees smooths out variability across individual trees, ensuring a capacity to capture complex conditional structures. The classical conditional quantile estimator is finally given by

$$ \tilde{Q}^{(\text{RF})}(\alpha | x) := \arg \min_{q \in \mathcal{Y}} \sum_{i=1}^{n} W_{\alpha}^{(\text{RF})}(x, X_i) \rho_{\alpha}(Y_i - q). \quad (8) $$

This method can be extended to estimate multiple quantiles simultaneously at different levels $\alpha_1 > \dots > \alpha_K > \alpha_{K+1} := 0$. For a given node P, instead of partitioning the data into two groups based on a single conditional quantile estimate, observations are divided into K intervals

$$ I_k(P) := [\hat{Q}(\alpha_k | P), \hat{Q}(\alpha_{k+1} | P)), \quad k = 1, \dots, K. $$

For the b-th tree, the optimal split $(C_1, C_2 = P \setminus C_1)$ of node P maximizes the multiple quantiles criterion:

$$ \sum_{k=1}^{K} \sum_{\ell=1}^{2} \frac{1}{|C_{\ell}|} \left( \sum_{X_i^{(b)} \in C_{\ell}} \mathbb{I}\{Y_i^{(b)} \in I_k(P)\} \right)^2. \quad (9) $$

This generalized random forest algorithm is used for two different purposes. First, it is employed to estimate the classical conditional quantile of order $\alpha_n$ which is required to get the empirical excesses in Sections 3.1 and 3.2 and appears in the expression of the extreme conditional quantile estimators. Next, the generalized random forest algorithm with a multiple quantile criterion is applied to get the similarity weights in the negative log-likelihood, see Section 3.1. We end this section by giving in Table 1 the hyper-parameters required to perform this estimation procedure of $Q(\alpha_n | x)$ together with a proposition of selection.

<table><caption>Table 1: Hyper-parameters for the estimation of a classical conditional quantile.</caption><thead><tr><td>Hyper-parameter</td><td>Description</td><td>Selection</td></tr></thead><tbody><tr><td><i>B</i></td><td>number of trees</td><td><i>B</i> = 2000</td></tr><tr><td><i>s</i><sub>n</sub>(<i>α</i>)</td><td>size of the subsamples</td><td><i>s</i><sub>n</sub>(<i>α</i>) = |<i>n</i>/2|</td></tr><tr><td><i>m</i><sub>RF</sub>(<i>α</i>)</td><td>minimal number of points if a node</td><td><i>m</i><sub>RF</sub>(<i>α</i>) = 5</td></tr></tbody></table>

## 3 Tail parameters estimation

After estimating the classical conditional quantile, the second step consists in estimating the tail parameters. We present here three main approaches: extreme random forest, extreme neural networks and dimension reduction techniques. Given these findings, we focus on the most competitive methods, which have demonstrated superior performance in high-dimensional settings.
## 3.1 Extreme random forest

The extreme random forest estimator was proposed by Gneco et al. [15]. It is based on a result of Pickands [22], which states that under model (M), the excess distribution above a high threshold converges to a Generalized Pareto Distribution (GPD): for all $z > 0$,

$$ \lim_{\alpha \to 0^+} \mathbb{P} \left( \frac{Y - Q(\alpha | x)}{a(\alpha^{-1} | x)} < z \mid Y > Q(\alpha | x), X = x \right) = G(z, 1, \gamma(x)), \quad (10) $$

where $G(z, a, \gamma) := 1 - (1 + \gamma z/a)^{-1/\gamma}$ is the cumulative distribution function of a GPD with parameters $(a, \gamma)$. The corresponding negative log-density function is given for $z > 0$ by

$$ g(z, a, \gamma) := -\log \frac{\partial G(z, a, \gamma)}{\partial z} = \log(a) + \left(1 + \frac{1}{\gamma}\right) \log \left(1 + \frac{\gamma}{a}z\right). \quad (11) $$

The estimation procedure is presented hereafter. For a given classical conditional quantile level $\alpha_n$, let

$$ \{\hat{Y}_i(\alpha_n); i = 1, \dots, M_n(\alpha_n)\}, $$

be the random variables $Y_i$'s larger than the estimated classical conditional quantile $\check{Q}^{(\text{RF})}(\alpha_n | X_i)$ which is defined in (8). The random variables $\hat{Y}_i(\alpha_n) - \check{Q}^{(\text{RF})}(\alpha_n | X_i)$ are referred to as the empirical excesses. If a large weight is assigned to the empirical excesses for which the associated covariate $X_i$ is such that the conditional distribution of Y given X = $X_i$ is similar to the one of Y given X = $x$, we get a sample that, according to (10), follows approximatively a GPD with parameters $(a_n(x), \gamma(x))$. The difference between the two conditional distributions is measured by the similarity weights

$$ \{W_{\vec{\alpha}}^{(\text{RF})}(x, X_i); i = 1, \dots, n\}, $$

which are obtained with the generalized random forest algorithm described in Section 2 by using the multiple quantiles criterion (9) with $\vec{\alpha} = (0.9, 0.5, 0.1)$. As in Gneco et al. [15], this choice of $\vec{\alpha}$ permits to cover a large part of the conditional distribution of $Y | X$. This leads us to the definition of the estimator of an extreme conditional quantile as introduced in Gneco et al. [15]. Let us first give some notations.

For $i \in \{1, \dots, n\}$, let $U_i := (Y_i - \check{Q}^{(\text{RF})}(\alpha_n | X_i))_+$ where $(\cdot)_+$ is the positive part function. Let also introduced the estimators

$$ (\hat{a}_n, \hat{\gamma}_n) = \arg \min_{(a, \gamma)} \sum_{i=1}^{n} g(Y_i - Q_Y(\alpha_n), a, \gamma) \mathbb{I}\{Y_i > Q_Y(\alpha_n)\}, \quad (12) $$

where

$$ Q_Y(\alpha_n) := \arg \min_{q \in \mathcal{Y}} \sum_{i=1}^{n} \rho_{\alpha_n}(Y_i - q). \quad (13) $$
Note that $\hat{a}_n$ and $\hat{\gamma}_n$ are the estimators of the tail-parameters obtained in the unconditional case that is to say when no covariate $X$ is recorded with the response variable $Y$.

**Definition 1.** Let $x \in \mathcal{X}$ and let $\alpha_n$ be a classical conditional quantile level. For $\beta_n < \alpha_n$, the random forest (RF) estimator of the extreme conditional quantile $Q(\beta_n | x)$ is given by

$$ \hat{Q}^{(RF)}(\beta_n | x) := \tilde{Q}^{(RF)}(\alpha_n | x) + \frac{\hat{a}_n^{(RF)}(x)}{\hat{\gamma}_n^{(RF)}(x)} \left[ \left( \frac{\alpha_n}{\beta_n} \right)^{\hat{\gamma}_n^{(RF)}(x)} - 1 \right], $$

where the estimated tail-parameters $\hat{a}_n^{(RF)}(x)$ and $\hat{\gamma}_n^{(RF)}(x)$ minimize the penalized negative log-likelihood

$$ \sum_{i=1}^{n} W_{\tilde{\alpha}}^{(RF)}(x, X_i)g(U_i, a, \gamma)\mathbb{I}\{U_i > 0\} + \lambda(\gamma - \hat{\gamma}_n)^2, $$

for some regularization parameter $\lambda > 0$ and with $\hat{\gamma}_n$ defined in (12).

The penalization is introduced to prevent excessive fluctuations of the conditional tail index estimator across different values of $x$. Large variations can lead to erratic behavior in extreme quantile estimates, particularly in regions where data is sparse. Since the fluctuation of the scale parameter estimator is, to a large extent, due to that of the extreme value index estimator, we do not found useful to introduce an additional term of penalization.

A heavy-tailed variant with the log-exceedances is also proposed by Gneco et al. [15]. It is based on the following asymptotic result (see de Haan and Ferreira [16, Theorem 3.2.2]):

$$ \lim_{\alpha \to 0_+} \mathbb{P} \left( \log \frac{Y}{\tilde{Q}(\alpha | x)} < z \mid Y > Q(\alpha | x), X = x \right) = \mathcal{E}(1/\gamma(x)), $$

where $\mathcal{E}(1/\gamma(x))$ denotes the cumulative distribution function of an exponential distribution with mean $\gamma(x)$. This suggests an alternative estimation strategy by approximating the weighted expectation of the log-exceedances $\log \hat{Y}_i(\alpha_n) - \log \hat{Q}^{(RF)}(\alpha_n | X_i)$ by $\gamma(x)$. More specifically, introducing the random variables

$$ V_i := \left( \log \frac{Y_i}{\tilde{Q}^{(RF)}(\alpha_n | X_i)} \right)_+ $$

the following estimator was introduced in Gneco et al. [15].

**Definition 2.** Let $x \in \mathcal{X}$ and let $\alpha_n$ be a classical conditional quantile level. For $\beta_n < \alpha_n$, the heavy-tailed version of the random forest estimator (HRF) of the extreme conditional quantile $Q(\beta_n | x)$ is given by

$$ \hat{Q}^{(HRF)}(\beta_n | x) = \tilde{Q}^{(RF)}(\alpha_n | x) \left( \frac{\alpha_n}{\beta_n} \right)^{\hat{\gamma}_n^{(HRF)}(x)}, $$
where the estimated tail index parameter is given by:

$$ \hat{\gamma}_n^{(HRF)}(x) \in \arg \min_{\gamma>0} \sum_{i=1}^{n} W_{\hat{\alpha}}^{(RF)}(x, X_i)(V_i - \gamma)^2 \mathbb{I}\{V_i > 0\} + \lambda(\gamma - \hat{\gamma}_n)^2, $$

for some regularization parameter $\lambda > 0$ and with $\hat{\gamma}_n$ defined in (12).

In Table 2, we give all the parameters required for the extreme random forest algorithm.

<table><caption>Table 2: Hyper-parameters for the extreme random forest algorithm.</caption><thead><tr><th colspan="3">For the estimator Q<sup>(RF)</sup>(α<sub>n</sub> | x)</th></tr><tr><th>Hyper-parameter</th><th>Description</th><th>Selection</th></tr></thead><tbody><tr><td>α<sub>n</sub></td><td>classical quantile level</td><td>Choice left to the user</td></tr><tr><td>B</td><td>number of trees</td><td>B = 2000</td></tr><tr><td>s<sub>n</sub>(α<sub>n</sub>)</td><td>size of the subsamples</td><td>s<sub>n</sub>(α<sub>n</sub>) = [n/2]</td></tr><tr><td>m<sub>RF</sub>(α<sub>n</sub>)</td><td>minimal number of points if a node</td><td>m<sub>RF</sub>(α<sub>n</sub>) = 5</td></tr></tbody><thead><tr><th colspan="3">For the similarity weights W<sub>α</sub><sup>(RF)</sup>(x, X<sub>i</sub>)</th></tr><tr><th>Hyper-parameter</th><th>Description</th><th>Selection</th></tr></thead><tbody><tr><td>B</td><td>number of trees</td><td>B = 2000</td></tr><tr><td>s<sub>n</sub>(α̂)</td><td>size of the subsamples</td><td>s<sub>n</sub>(α̂) = [n/2]</td></tr><tr><td>m<sub>RF</sub>(α̂)</td><td>minimal number of points if a node</td><td>cross-validation</td></tr></tbody><thead><tr><th colspan="3">Estimation of the tail parameter(s)</th></tr><tr><th>Hyper-parameter</th><th>Description</th><th>Selection</th></tr></thead><tbody><tr><td>λ</td><td>regularization parameter</td><td>cross-validation</td></tr></tbody></table>

The random forest approach is advantageous for extreme quantile estimation due to its theoretical guarantees and computational efficiency. Gnecco et al. [15] established the consistency under reasonable assumptions, ensuring reliable estimation as the sample size grows. Additionally, the method is computationally efficient, making it well-suited for large datasets.

## 3.2 Extreme neural networks

Neural networks are powerful models for capturing complex, non-linear relationships in data, and their use has become increasingly popular across a range of fields, including extreme value theory in recent years (Allouche et al. [3, 4]). In this section, we focus on the Pasche and Engelke [21] approach's. The field of neural networks is vast, and many additional methods exist that are beyond the scope of this work. For a more comprehensive introduction to neural networks, we refer the reader to classic references such as Aggarwal [1].
<figure><img src="image_3.png" /><figcaption>Figure 1: Schematic representation of an FNN with two hidden layers. The operation between each layer is denoted above the arrows.</figcaption></figure>

### 3.2.1 Introduction to feedforward neural networks

A feedforward neural networks (FNN) maps an input $e \in \mathbb{R}^p$ to an output $\hat{o}(e) \in \mathbb{R}^{p_o}$. It consists of $L-1$ hidden layers and 1 output layer, each performing a linear transformation followed by a non-linear activation function. Formally, for each hidden layer $l \in \{1, \dots, L-1\}$, the transformation can be written as:

$$h^{(l)} = \sigma \left( \theta^{(l)} h^{(l-1)} + b^{(l)} \right) \in \mathbb{R}^{d_l},$$

where $\theta^{(l)} \in \mathbb{R}^{d_l \times d_{l-1}}$ is the weight matrix (with $d_0 = p$), $b^{(l)} \in \mathbb{R}^{d_l}$ is the bias vector, and $\sigma$ is the activation function. Note that $\sigma$ applies to a vector element-wise i.e., $\sigma((x_1, \dots, x_d)^\top) = (\sigma(x_1), \dots, \sigma(x_d))^\top$. The input for the first hidden layer is $h^{(0)} = e$. The output layer produces the final prediction of dimension $p_o$:

$$\hat{o}(e) = f \left( \theta^{(L)} h^{(L-1)} + b^{(L)} \right) \in \mathbb{R}^{p_o}.$$

where $f$ is an arbitrary function, $\theta^{(L)} \in \mathbb{R}^{p_o \times d_{L-1}}$ and $b^{(L)} \in \mathbb{R}^{p_o}$. Figure 1 illustrates the architecture of a simple FNN with two hidden layers ($L=3$).

The output function $\hat{o}(\cdot) = \hat{o}(\cdot | W)$ is thus a quite complex function depending on a vector $W$ of parameters. This vector contains the coefficients of the matrices $\{\theta^{(l)}, l = 1, \dots, L\}$ and the components of the vectors $\{b^{(l)}, l = 1, \dots, L\}$ so that the dimension of $W$ is

$$D_{FNN} := \sum_{l=1}^{L} (d_{l-1} + 1)d_l,$$
with $d_L = p_o$. Given a training set $\{(e_1, o_1), \dots, (e_n, o_n)\}$, the parameter vector of the output function is denoted by $\hat{\mathcal{W}}_n$ and is obtained by maximizing the average loss function:

$$ \mathcal{AL}(\mathcal{W}) := \frac{1}{n} \sum_{i=1}^{n} \mathcal{L}(\hat{o}(e_i | \mathcal{W}), o_i), $$

where $\mathcal{L}(\cdot, \cdot)$ is a given loss function that depends on the specific context in which the FNN is used (e.g., regression, classification, etc.).

This optimization problem is solved by using a gradient descent algorithm: for a given learning rate $\eta$, the parameter vector is updated by using the rule

$$ \mathcal{W} \leftarrow \mathcal{W} - \eta \nabla_{\mathcal{W}}(\mathcal{AL}(\mathcal{W})). \qquad (14) $$

To be computationally efficient, the gradient needs to be computed iteratively, layer by layer, starting from the last layer and propagating back to the first (backpropagation algorithm). This update rule is repeated until convergence (with respect to a criterion choose by the user) or until a given number of iterations is reached. In practice, instead of computing the gradient for each point of the training set, which is computationally demanding, we randomly divide the training set into $n/s_B$ mini-batches $\{B_j, j = 1, \dots, n/s_B\}$ of size $s_B \in \{1, \dots, n\}$. We then perform iteratively from $B_1$ to $B_{n/s_B}$ the update rule (14) by replacing the gradient $\nabla_{\mathcal{W}}(\mathcal{AL}(\mathcal{W}))$ by

$$ \frac{1}{s_B} \sum_{i \in B_j} \nabla_{\mathcal{W}} (\mathcal{L}(\hat{o}(e_i | \mathcal{W}), o_i)). $$

A complete pass through all mini-batches $\{B_j, j = 1, \dots, n/s_B\}$ is referred to as an epoch. In general, several epochs (say $n_e > 1$) are required to ensure convergence. Let us emphasize that the choice of the mini-batch size $s_B$ plays an important role in training dynamics. Typically, values of $s_B$ range from 32 to 512, depending on the dataset size $n$ and computational constraints.

### 3.2.2 Estimation of an extreme conditional quantile

When observations $\{(x_1, y_1), \dots, (x_n, y_n)\}$ are drawn from an independent sample $\{(X_1, Y_1), \dots, (X_n, Y_n)\}$, Pasche and Engelke [21] propose to estimate the tail parameters $a_n(x)$ and $\gamma(x)$ by the output of an FNN. Their method, referred to as the Extreme Quantile Regression Networks, was initially tailored for flood risk forecasting in the Swiss Aare river catchment. For $\zeta = (\nu, \gamma) \in (0, \infty) \times (-1/2, \infty)$, the FNN is trained by using the loss function

$$ \mathcal{L}^{(\text{GPD})}(\zeta, u) := \left[ \left(1 + \frac{1}{\gamma}\right) \log\left(1 + \gamma \frac{u}{\nu}\right) + \log(\nu) - \log(1 + \gamma) \right] \mathbb{I}\{u > 0\}. $$

Note that $\mathcal{L}^{(\text{GPD})}(\zeta, u)$ is a reparameterization of the negative log-density of a GPD with parameters $a$ and $\gamma$ (see (11)) where $\nu = a(1 + \gamma)$. The authors have
observed that this new parameterization leads to a more stable estimation of the GPD parameters. The training set is $\{(e_1, o_1), \dots, (e_n, o_n)\}$ with $e_i = x_i$ and $o_i = u_i$ where

$$u_i := (y_i - \check{Q}^{(\text{RF})}(\alpha_n | x_i))_+,$$

with $\check{Q}^{(\text{RF})}(\alpha_n | x_i)$ is defined in (8) for some classical conditional quantile level $\alpha_n$. For an input $x$, denoting by $\mathcal{W}$ the parameter vector of the FNN, the output is the two-dimensional vector

$$\zeta(x | \mathcal{W}) = (\hat{\nu}(x | \mathcal{W}), \hat{\gamma}(x | \mathcal{W})).$$

Following the recommendation of Pasche and Engelke [21], this output is obtained by applying the activation function $f$ given for $x \in \mathbb{R}$ by

$$f(x) = (\exp(x), 0.6 \tanh(x) + 0.1),$$

so that the scale-related parameter $\hat{\nu}(x) \in (0, +\infty)$, and the tail index estimate $\hat{\gamma}(x) \in (-0.5, 0.7)$. Using the procedure described in Section 3.2, the parameter vector of the trained FNN is the vector $\hat{W}_n^{(\text{GPD})}$ minimizing the penalized average loss

$$\mathcal{A}\mathcal{L}^{(\text{GPD})}(\mathcal{W}, \lambda) := \frac{1}{n} \sum_{i=1}^{n} \mathcal{L}^{(\text{GPD})} \left( \hat{\zeta}(x_i | \mathcal{W}), u_i \right) + \lambda \|\mathcal{W}\|_2, \quad (15)$$

for some regularization parameter $\lambda > 0$ and where $\| \cdot \|_2$ is the $L^2$-norm. The penalty term has been added to improve the generalization power of the FNN. We are now in a position to give the definition of the estimator of the extreme conditional quantile $Q(\beta_n | x)$.

**Definition 3.** Let $x \in \mathcal{X}$ and let $\alpha_n$ be a classical conditional quantile level. For $\beta_n < \alpha_n$, the neural networks (NN) estimator of the extreme conditional quantile $Q(\beta_n | x)$ is given by

$$\hat{Q}^{(NN)}(\beta_n | x) = \check{Q}^{(RF)}(\alpha_n | x) + \frac{\hat{a}_n^{(NN)}(x)}{\hat{\gamma}_n^{(NN)}(x)} \left[ \left( \frac{\alpha_n}{\beta_n} \right)^{\hat{\gamma}_n^{(NN)}(x)} - 1 \right],$$

where the estimated parameters are given by

$$\hat{a}_n^{(NN)}(x) = \frac{\hat{\nu}(x | \hat{W}_n^{(GPD)})}{1 + \hat{\gamma}(x | \hat{W}_n^{(GPD)})}, \quad \text{and} \quad \hat{\gamma}_n^{(NN)}(x) = \hat{\gamma}(x | \hat{W}_n^{(GPD)}),$$

where $\hat{W}_n^{(GPD)}$ minimizes (15) for some regularization parameter $\lambda > 0$.

We finally propose to adapt the above estimator to the situation where the conditional tail index is positive, that is, for heavy-tailed conditional distributions. Note that this adaptation was not done in the paper by Pasche and Engelke [21]. In this case, the neural networks is trained using the loss function

$$\mathcal{L}^{(\text{HT})}(\gamma, v) = (v - \gamma)^2 \mathbb{I}\{v > 0\},$$
with the training set $\{(e_1, o_1), \dots, (e_n, o_n)\}$ where $e_i = x_i$,

$$o_i := v_i := \log \left( \frac{y_i}{\hat{Q}^{(\text{RF})}(\alpha_n | x_i)} \right),$$

and where $\alpha_n$ is a classical conditional quantile level. For an input $x$, we denote by $\hat{\gamma}(x | \mathcal{W})$ the output of the FNN calculated with the parameter vector $\mathcal{W}$. Here, the activation function for the output layer is $f(x) = x$. As before, the trained parameter vector is $\hat{W}_n^{(\text{HT})}$ which minimizes the penalized average loss

$$\mathcal{A}\mathcal{L}^{(\text{HT})}(\mathcal{W}, \lambda) := \frac{1}{n} \sum_{i=1}^{n} \mathcal{L}^{(\text{HT})} (\hat{\gamma}(x_i | \mathcal{W}), v_i) + \lambda \|\mathcal{W}\|_2, \quad (16)$$

for some positive regularization parameter $\lambda$. The estimator of the extreme conditional quantile for heavy-tailed distributions is given hereafter.

**Definition 4.** Let $x \in \mathcal{X}$ and let $\alpha_n$ be a classical conditional quantile level. For $\beta_n < \alpha_n$, the heavy-tailed version of the NN estimator (HNN) is given by:

$$\hat{Q}^{(\text{HNN})}(\beta_n | x) = \hat{Q}^{(\text{RF})}(\alpha_n | x) \left(\frac{\alpha_n}{\beta_n}\right)^{\hat{\gamma}_n^{(\text{HNN})}(x)},$$

where the estimated tail index parameter is given by:

$$\hat{\gamma}_n^{(\text{HNN})}(x) = \hat{\gamma}(x | \hat{W}_n^{(\text{HT})}),$$

where $\hat{W}_n^{(\text{HT})}$ minimizes (16) for some regularization parameter $\lambda > 0$.

Neural networks-based methods offer high flexibility and do not require strong distributional assumptions. However, the choice of architecture and regularization is crucial to avoid overfitting; see Table 3 where all the required hyper-parameters are listed. Furthermore, they lack theoretical guarantees regarding consistency and convergence, which can make their results harder to interpret compared to classical statistical methods.

## 3.3 Dimension reduction procedures

While machine learning methods generally perform well in high-dimensional settings, reducing the dimensionality of the covariate space can sometimes give better results. The basic idea is to add a constraint to the model for the couple $(X, Y)$ by assuming that all the information contained in the covariate $X$ about the tail behavior of $Y$ is captured by a lower-dimensional projection of the predictor vector, $B^\top X \in \mathbb{R}^{p \times q}$, with $q < p$. Hence, if $B$ is known (or can be estimated), inference can be performed by replacing $X$ with the lower-dimensional covariate $B^\top X$. When $q \ll p$, dimension reduction methods allows us to circumvent the curse of dimensionality.
<table><caption>Table 3: Hyper-parameters for the extreme neural networks algorithm.</caption><thead><tr><th colspan="3">For the estimator Q<sup>(RF)</sup>(α<sub>n</sub> | x)</th></tr><tr><th>Hyper-parameter</th><th>Description</th><th>Selection</th></tr></thead><tbody><tr><td>α<sub>n</sub></td><td>classical quantile level</td><td>Choice left to the user</td></tr><tr><td>B</td><td>number of trees</td><td>B = 2000</td></tr><tr><td>s<sub>n</sub>(α<sub>n</sub>)</td><td>size of the subsamples</td><td>s<sub>n</sub>(α<sub>n</sub>) = [n/2]</td></tr><tr><td>m<sub>RF</sub>(α<sub>n</sub>)</td><td>minimal number of points if a node</td><td>m<sub>RF</sub>(α<sub>n</sub>) = 5</td></tr><tr><th colspan="3">For the stochastic gradient descent algorithm</th></tr><tr><th>Hyper-parameter</th><th>Description</th><th>Selection</th></tr><tr><td>η</td><td>learning rate</td><td>η = 10<sup>-4</sup></td></tr><tr><td>n<sub>e</sub></td><td>number of epochs</td><td>Early stopping*</td></tr><tr><td>s<sub>B</sub></td><td>size of mini-batches</td><td>128</td></tr><tr><th colspan="3">Architecture of the neural networks</th></tr><tr><th>Hyper-parameter</th><th>Description</th><th>Selection</th></tr><tr><td>L</td><td>number of layers</td><td>cross-validation</td></tr><tr><td>σ</td><td>activation function for the hidden layers</td><td>σ = tanh</td></tr><tr><td>d<sub>1</sub>, ..., d<sub>L-1</sub></td><td>dimensions of the weight matrices</td><td>cross-validation</td></tr><tr><th colspan="3">Estimation of the tail parameter(s)</th></tr><tr><th>Hyper-parameter</th><th>Description</th><th>Selection</th></tr><tr><td>λ</td><td>regularization parameter</td><td>cross-validation</td></tr></tbody></table>

To our knowledge, three dimension reduction methods have been developed for extreme values. The earliest, proposed by Gardes [12], is computationally demanding and its implementation is not easy when $q > 1$. Two more recent approaches will be presented in this section: the tail inverse regression proposed by Aghbalou et al. [2] and the method based on the central tail index subspace introduced in Gardes and Podgorny [13].

### 3.3.1 Tail inverse regression

Aghbalou et al. [2] propose two dimension reduction methods, TIREX1 and TIREX2, which can be seen as adaptations of Sliced Inverse Regression (SIR) (Li [19]) and SAVE (Cook and Weisberg [11]), respectively. However, TIREX1 is theoretically limited to estimating subspaces of dimension at most 1. Hereafter, we chose to focus on TIREX2 method.

The TIREX2 method estimates a basis for the Extreme Sufficient Direction Reduction (Extreme SDR) subspace, defined as follows:

**Definition 5.** Let $y^*$ be the right endpoint of the distribution of Y. A linear

<footnote>* The early stopping procedure is described in Section 4.</footnote>
subspace $S_E$ spanned by a basis $B \in \mathbb{R}^{p \times q}$ is an Extreme SDR if

$$ \lim_{y \to y^*} \mathbb{E} \left| \frac{\mathbb{P}(Y > y | B^\top X) - \mathbb{P}(Y > y | X)}{\mathbb{P}(Y > y)} \right| = 0. \quad (17) $$

Moreover, $S_E$ is called a central Extreme SDR subspace if it is an Extreme SDR subspace containing all the Extreme SDR subspaces.

If $B$ spans an Extreme SDR, the tails of the conditional distributions $Y | X = x$ and $Y | B^\top X = B^\top x$ are, roughly speaking, close to each other. However, Equation (17) does not imply that $Q(\alpha | x)/Q(\alpha | B^\top x) \to 1$ as $\alpha \to 0$. The authors' primary focus was not on extreme quantile estimation, but our simulations suggest that the method can be effectively used for this purpose in practice. Note also that the central Extreme SDR subspace does not necessarily exist but if it is the case, it is unique and given by the intersection of all Extreme SDR subspaces. In all what follows, we assume that the central Extreme SDR $S_E$ exists.

Let us now give the main idea leading to the estimator of the linear subspace $S_E$ proposed by Aghbalou et al. [2]. Let $Z := \Sigma^{-1/2}(X - m)$ be the normalized random vector of covariates with $m = \mathbb{E}(X)$ and $\Sigma = \mathbb{V}(X)$ and let $P_E$ be the orthogonal projection matrix on the subspace $S_E$. Under the linearity condition $\mathbb{E}(Z | P_E Z) = P_E Z$ and if $\mathbb{V}(Z | P_E Z)$ is almost surely constant, the authors establish, under additional technical conditions (see Aghbalou et al. [2, Theorem 1 and 2]), that the $p \times p$ matrix

$$ M_T := \int_0^1 D_n(u)[D_n(u)]^\top du, $$

with, denoting by $I_p$ the identity matrix and by $Q_Y$ the quantile function of $Y$ (see (13)),

$$ D_n(u) := \frac{1}{\alpha_n} \mathbb{E} \left[ (ZZ^\top - I_p) \mathbb{I}\{Y > Q_Y(\alpha_n)\} \right], $$

is such that $\text{span}(M_T) \subset S_E$. The TIREX estimator of the central Extreme SDR subspace is the subspace spanned by the $q$ largest eigenvalues of the estimated matrix:

$$ \hat{M}_n := \int_0^1 \hat{D}_n(u) [\hat{D}_n(u)]^\top du, $$

where

$$ \hat{D}_n(u) := \Sigma_n^{-1/2} \left( \frac{1}{[n\alpha_n]} \sum_{i=1}^n \overline{X_i} \overline{X_i}^\top \mathbb{I}\{Y_i > Y_{n-[n\alpha_n]+1,n}\} \right) \Sigma_n^{-1/2}, $$

with $Y_{1,n} \le \dots \le Y_{n,n}$ the order statistics associated to the sample $\{Y_1, \dots, Y_n\}$,

$$ \overline{X}_i := X_i - \frac{1}{n} \sum_{j=1}^n X_j \text{ and } \Sigma_n = \frac{1}{n} \sum_{i=1}^n \left( X_i - \frac{1}{n} \sum_{j=1}^n X_j \right) \left( X_i - \frac{1}{n} \sum_{j=1}^n X_j \right)^\top. $$
Denoting by $\hat{B}_n^{(T)}$ the $p \times q$ matrix of the $q$ largest eigenvalues of $\hat{M}_n$, we proposed the following kernel-based estimator of $Q(\beta_n | x)$ when the conditional distribution of $Y$ given $X = x$ is heavy-tailed.

**Definition 6.** Let $x \in \mathcal{X}$ and let $\alpha_n$ be a classical conditional quantile level. For $\beta_n < \alpha_n$ and when $\gamma(x) > 0$, the TIREX estimator is given by:

$$ \hat{Q}^{(T)}(\beta_n | x) = \check{Q}^{(RF)}(\alpha_n | x) \left(\frac{\alpha_n}{\beta_n}\right)^{\hat{\gamma}_n^{(T)}(x)}, \quad (18) $$

where the estimated tail index is given for a non-negative function $K$ on $\mathbb{R}^q$ and a bandwidth $h_n$ by:

$$ \hat{\gamma}_n^{(T)}(x) \in \arg \min_{\gamma \in \mathbb{R}} \sum_{i=1}^{n} K \left( \frac{[\hat{B}_n^{(T)}]^\top (x - X_i)}{h_n} \right) (V_i - \gamma)^2 \mathbb{I}_{V_i > 0}. $$

The required hyperparameters to perform the TIREX approach are listed in Table 4.

<table><caption>Table 4: Hyper-parameters for the TIREX method.</caption><thead><tr><th colspan="3">For the estimator Q<sup>(RF)</sup>(α<sub>n</sub> | x)</th></tr><tr><th>Hyper-parameter</th><th>Description</th><th>Selection</th></tr></thead><tbody><tr><td>α<sub>n</sub></td><td>classical quantile level</td><td>user choice</td></tr><tr><td>B</td><td>number of trees</td><td>B = 2000</td></tr><tr><td>s<sub>n</sub>(α<sub>n</sub>)</td><td>size of the subsamples</td><td>s<sub>n</sub>(α<sub>n</sub>) = [n/2]</td></tr><tr><td>m<sub>RF</sub>(α<sub>n</sub>)</td><td>minimal number of points if a node</td><td>m<sub>RF</sub>(α<sub>n</sub>) = 5</td></tr></tbody><thead><tr><th colspan="3">Estimation of the basis B̂<sub>n</sub><sup>(T)</sup></th></tr><tr><th>Hyper-parameter</th><th>Description</th><th>Selection</th></tr></thead><tbody><tr><td>q</td><td>dimension of the subspace</td><td>cross-validation</td></tr></tbody><thead><tr><th colspan="3">For the Kernel method</th></tr><tr><th>Hyper-parameter</th><th>Description</th><th>Selection</th></tr></thead><tbody><tr><td>K</td><td>kernel density</td><td>K(·) = I{||·||<sub>∞</sub> &lt; 1}</td></tr><tr><td>h<sub>n</sub></td><td>bandwidth parameter</td><td>h<sub>n</sub> = n<sup>-0.2/q</sup></td></tr></tbody></table>

### 3.3.2 Central Tail Index Subspace

Unlike TIREX, Gardes and Podgorny [13] propose a dimension reduction approach which is specifically designed for the estimation of the conditional tail index in the heavy-tailed case. In this section, we assume that for all $x \in \mathcal{X}$, the conditional distribution function of $Y$ given $X = x$ is heavy-tailed with tail index $\gamma(x) > 0$. Under additional assumptions, the authors prove that for
all $B \in \mathbb{R}^{p \times q}$, the conditional distribution of $Y$ given $B^\top X = B^\top x$ remains heavy-tailed with tail index

$$ \xi_B(B^\top x) := \max_{z: B^\top z = B^\top x} \gamma(z). \qquad (19) $$

We then define the Tail Index Dimension Reduction (TIDR) subspace as follows.

**Definition 7.** A linear subspace $S$ spanned by a basis $B \in \mathbb{R}^{p \times q}$ is a Tail Index Dimension Reduction (TIDR) subspace if the conditional distributions of $Y | X$ and $Y | B^\top X$ share the same conditional tail index, that is to say if $\xi_B(B^\top X) = \gamma(X)$ almost surely. Moreover, $S$ is called a central Tail Index (CTI) subspace if it is an TIDR subspace containing all the TIDR subspaces.

In what follows, we assume that the CTI subspace exists. It is thus unique as the intersection of all the TIDR subspaces. We denote by $B_0 \in \mathcal{B}_q$ its canonical basis (see Gardes and Podgorny [13, Definition 2]). Here, $\mathcal{B}_q$ denotes the set that gathers the canonical bases of all linear subspaces of dimension $q$. Note that for all $x \in \mathcal{X}$ and $B \in \mathcal{B}_q$,

$$ \gamma(x) = \xi_{B_0}(B_0^\top x) \le \xi_B(B^\top x) = \max_{z: B^\top z = B^\top x} \gamma(z), $$

since $\gamma(x)$ belongs to the support of the conditional distribution of $\gamma(X)$ given $B^\top X = B^\top x$. Hence, for a compact subset $\mathcal{X}_0$ in the interior of $\mathcal{X}$ satisfying the condition

$$ (C) \text{ For } B \in \mathcal{B}_q, \text{ if } \xi_B(B^\top x) = \gamma(x) \text{ for all } x \in \mathcal{X}_0 \text{ then } B = B_0, $$

we get that

$$ \arg\min_{B \in \mathcal{B}_q} \mathbb{E} \left[ \xi_B(B^\top X) \mathbb{I}\{X \in \mathcal{X}_0\} \right] =: \arg\min_{B \in \mathcal{B}_q} \Psi(B, \mathcal{X}_0) = \{B_0\}. $$

Under condition (C), we propose to estimate $B_0$ by the canonical basis $\hat{B}_n^{(\text{CTI})}$ which minimizes the empirical version of $\Psi(B, \mathcal{X}_0)$ given by:

$$ \hat{\Psi}_n(B, \mathcal{X}_0) := \frac{1}{n} \sum_{i=1}^{n} \hat{\xi}_{n,B}(B^\top X_i) \mathbb{I}\{X_i \in \mathcal{X}_0\}. $$

Similarly to TIREX, we define the following kernel-based estimator of an extreme conditional quantile.

**Definition 8 (CTI).** Let $x \in \mathcal{X}$ and let $\alpha_n$ be a classical conditional quantile level. For $\beta_n < \alpha_n$, the CTI estimator of the conditional extreme quantile is given by:

$$ \hat{Q}^{(\text{CTI})}(\beta_n | x) = \check{Q}^{(\text{RF})}(\alpha_n | x) \left( \frac{\alpha_n}{\beta_n} \right)^{\hat{\gamma}_n^{(\text{CTI})}(x)}, \qquad (20) $$

where the estimated tail index is given for a non-negative $\mathcal{K}$ on $\mathbb{R}^q$ and a band-width $h_n$ by:

$$ \hat{\gamma}_n^{(\text{CTI})}(x) \in \arg\min_{\gamma \in \mathbb{R}} \sum_{i=1}^{n} \mathcal{K} \left( \frac{[\hat{B}_n^{(\text{CTI})}]^\top (x - X_i)}{h_n} \right) (V_i - \gamma)^2 \mathbb{I}_{V_i > 0}. $$
We gather in Table 5 all the hyperparameters that are required for the CTI method.

<table><caption>Table 5: Hyper-parameters for the CTI method.</caption><thead><tr><th colspan="3">For the estimator ̃Q<sup>(RF)</sup>(α<sub>n</sub> | x)</th></tr><tr><th>Hyper-parameter</th><th>Description</th><th>Selection</th></tr></thead><tbody><tr><td>α<sub>n</sub></td><td>classical quantile level</td><td>Choice left to the user</td></tr><tr><td>B</td><td>number of trees</td><td>B = 2000</td></tr><tr><td>s<sub>n</sub>(α<sub>n</sub>)</td><td>size of the subsamples</td><td>s<sub>n</sub>(α<sub>n</sub>) = [n/2]</td></tr><tr><td>m<sub>RF</sub>(α<sub>n</sub>)</td><td>minimal number of points if a node</td><td>m<sub>RF</sub>(α<sub>n</sub>) = 5</td></tr></tbody></table><table><thead><tr><th colspan="3">Estimation of the basis ̂B<sub>n</sub><sup>(CTI)</sup></th></tr><tr><th>Hyper-parameter</th><th>Description</th><th>Selection</th></tr></thead><tbody><tr><td>q</td><td>dimension of the subspace</td><td>Procedure in Gardes and Podgorny [13]</td></tr></tbody></table><table><thead><tr><th colspan="3">For the Kernel method</th></tr><tr><th>Hyper-parameter</th><th>Description</th><th>Selection</th></tr></thead><tbody><tr><td>K</td><td>kernel density</td><td>K(·) = I{||·||<sub>∞</sub> &lt; 1}</td></tr><tr><td>h<sub>n</sub></td><td>bandwidth parameter</td><td>h<sub>n</sub> = n<sup>-0.2/q</sup></td></tr></tbody></table>

Dimension reduction methods provide an efficient way to handle high-dimensional covariate spaces by projecting the data onto a lower-dimensional subspace while preserving relevant information for tail index estimation. Additionally, it can enhance interpretability by identifying key directions that influence extreme events. However, they rely on the correct identification of the relevant subspace, which may not always be theoretically justified.

# 4 Simulation Study

In this section, we evaluate the performance of the various methods presented in this paper for estimating extreme conditional quantiles in a high-dimensional setting. The goal is to assess their accuracy and robustness under controlled conditions where the true quantiles are known. The study is structured as follows. Section 4.1 describes the general setup for the data-generating process and the evaluation protocol. Section 4.2 discusses the selection of tuning parameters for each method. Section 4.3 presents and analyzes the results of the experiments across different scenarios.

## 4.1 General Setup

We generate $R = 100$ independent replications of samples $(Y_1, X_1), \dots, (Y_n, X_n)$, each of size $n = 2000$. The covariates $X \in [0, 1]^p$ are generated from a Gaussian copula with correlation matrix $\Sigma_{ij} = \nu^{|i-j|}$ where $\nu \in [0, 1]$, ensuring an
AR(1)-type correlation structure between coordinates with uniform marginal distributions. Note that when $\nu = 0$, the components of $X$ are independent. The conditional distribution of $Y | X = x$ is assumed to be heavy-tailed, with a conditional quantile function given for $x = (x_1, \dots, x_p)^\top$ by

$$Q(u | x) = q(u, \gamma(x)) \times (1 + 6\phi(x_1, x_2)),$$

where $q(\cdot, \gamma(x))$ is a heavy-tailed quantile function of a distribution with tail index $\gamma(x) > 0$ and $\phi$ denotes the density of a standard bivariate normal distribution centered at (0,0) with correlation coefficient 0.9. Note that this conditional quantile function is also used by Pasche and Engelke [21] in their simulation study.

On each replication, we apply all the estimation methods presented in this paper to estimate the extreme conditional quantile of level $\beta_n = 2/n$. The performance of each method is evaluated by computing the estimation error on the set $\{x_i, i = 1, \dots, 1000\}$ of the first 1000 points of a Halton sequence in $[-1, 1]^p$. For each replication $r$, we evaluate the estimation accuracy using the symmetric relative error at each evaluation point $x_i$, defined as

$$E^{(r, \bullet)}(x_i) := \max \left( \left| \frac{\hat{Q}^{(\bullet)}(\beta_n | x_i)}{Q(\beta_n | x_i)} - 1 \right|, \left| \frac{Q(\beta_n | x_i)}{\hat{Q}^{(\bullet)}(\beta_n | x_i)} - 1 \right| \right),$$

with $\bullet \in \{\text{RF, HRF, NN, HNN, T, CTI}\}$. This criterion is scale-invariant and penalizes both underestimation and overestimation in a balanced way. The global error for the $r$-th replication is then defined as the median of the pointwise errors over the evaluation set:

$$\bar{E}^{(r, \bullet)} = \text{Med} \left( \{E^{(r, \bullet)}(x_i) : i = 1, \dots, 2000\} \right).$$

## 4.2 Choice of Tuning Parameters

All methods presented in this paper involve several hyperparameters that can significantly affect their performance. In this section, we describe the tuning procedures adopted for each approach. The hyper-parameter values used for each method are summarized in the tables provided in their respective sections. Unless otherwise stated, these values follow the recommendations of the authors who introduced the corresponding methods. When parameters require tuning, we adopt a standard three-fold cross-validation procedure. More precisely, the training data is randomly partitioned into three approximatively equally sized subsets. For each fold, two subsets are used to train the model, and the remaining subset is used to evaluate the performance. This process is repeated three times, rotating the validation subset at each iteration. The configuration minimizing the average of the three validation errors is selected.

For the random forest-based methods RF and HRF, we perform the cross-validation procedure over the grid of hyperparameters considered in Gnecco et al. [15].
Specifically, we are searching the regularization parameter $\lambda$ in the set $\{0, 0.01, 0.1\}$ and the minimum leaf size $m_{RF}(\vec{\alpha})$ in the set $\{5, 40, 100\}$, resulting in a total of nine combinations. For a subset of indices $F \subset \{1, \dots, n\}$, the validation error for the method RF with the fold $F$ is given by

$$CV^{(RF)}(\lambda, m_{RF}(\vec{\alpha})) := \sum_{i \in F} g(U_i, \hat{a}_n^{(RF)}(X_i), \hat{\gamma}_n^{(RF)}(X_i)) \mathbb{I}\{U_i > 0\},$$

where $\hat{a}_n^{(RF)}$ and $\hat{\gamma}_n^{(RF)}$ are obtained by the procedure described in Definition 1 with the observations $\{(X_i, Y_i); i \notin F\}$. For the method HRF, the validation error for the fold $F$ is

$$CV^{(HRF)}(\lambda, m_{RF}(\vec{\alpha})) := \sum_{i \in F} (V_i - \hat{\gamma}_n^{(HRF)}(X_i))^2 \mathbb{I}\{V_i > 0\},$$

where the estimator $\hat{\gamma}_n^{(HRF)}$ is obtained as in Definition 2 by using the subsample $\{(X_i, Y_i); i \notin F\}$.

For the neural networks based approaches, we use a cross-validation grid over combinations of network architecture and $L^2$ regularization. As Pasche and Engelke [21], we consider three architectures: shallow networks with three hidden layers of width 10, intermediate networks with two layers of width 64, and deeper networks with three layers of width 128. For each architecture, we test three values of penalty: $\lambda \in \{0, 10^{-5}, 10^{-3}\}$, leading again to nine configurations. Here, denoting by $L$ the number of layers, the validation error for a fold $F$ is

$$CV^{(NN)}(\lambda, L) := \sum_{i \in F} \mathcal{L}^{(GPD)}(\hat{\zeta}(X_i | \hat{W}_n^{(GPD)}), U_i),$$

for the NN approach and

$$CV^{(HNN)}(\lambda, L) := \sum_{i \in F} \mathcal{L}^{(HT)}(\hat{\gamma}(X_i | \hat{W}_n^{(HT)}), V_i),$$

where the parameters $\hat{W}_n^{(GPD)}$ and $\hat{W}_n^{(HT)}$ of the neural network are calculated by using the subsample $\{(X_i, Y_i); i \notin F\}$.

We also use an early stopping strategy to determine the number of epochs. More specifically, recall that in the three-fold cross-validation procedure, three validation errors are computed by rotating the fold dedicated to validation. The number of epochs used to obtain each validation error is selected as follows: the data belonging to the two training folds are randomly divided into a set $E_T$ (80% of the data) and a set $E_V$. The model is trained on $E_T$, and the validation error is evaluated on $E_V$. This evaluation is performed over several epochs, and the training is stopped when no improvement in the validation error is observed after a fixed number of consecutive epochs (referred to as the patience threshold). The model weights corresponding to the best validation performance are
retained. In our experiments, the patience threshold is set to 20 epochs. Once the best hyperparameters are selected through cross-validation, a final training is performed on the entire dataset, again applying early stopping with the same strategy.

For dimension reduction methods, the choice of the reduced dimension $q$ plays a central role. In the case of TIREX, $q$ is selected by cross-validation. For the CTI approach, $q$ is estimated using the specific procedure proposed by Gardes and Podgorny [13]. For the compact $\mathcal{X}_0$, we take $\mathcal{X}_0 := [\varepsilon_p, 1 - \varepsilon_p]^p$ with $\varepsilon_p = (1 - 0.9^p)/2$ so that approximately 90% of observations belong to $\mathcal{X}_0$. Finally, the choice of the classical quantile level $\alpha_n$ cannot be tuned via cross-validation, since it directly affects the number of training examples (i.e., the number of observed excesses). As a consequence, $\alpha_n$ must be chosen by the user based on prior knowledge or practical constraints. We investigate the impact of this choice in Experiment 1.

During the preparation of this work, we also investigated whether alternative hyperparameter configurations could significantly improve the performance of the compared methods. For the neural networks based approaches in particular, we experimented with deeper architectures, additional regularization strategies such as dropout, and different gradient-based optimization algorithms. These explorations did not lead to substantial performance gains, and we therefore retained the hyperparameter values recommended by the original authors. Further research in this direction would nevertheless be valuable, as the choice of hyperparameters can have a considerable influence on the accuracy and stability of these methods.

## 4.3 Experimental Results

We now present the results of four experiments designed to evaluate the sensitivity and robustness of the different methods under various conditions. Each experiment investigates a specific factor influencing estimation accuracy:

1) the threshold level $\alpha_n$;

2) the dimension $p$ of the covariate;

3) the second order parameter function $\rho$;

4) complexity of the conditional tail index function $\gamma$;

5) the correlation parameter $\nu$ between covariates.

For each experiment, we also look at the respective contributions of conditional extreme value index estimation and SDR subspace (for TIREX and CTI methods) to the overall performance of the extreme conditional estimators. More specifically, for each replication $r$ and method $\bullet$, let

$$ \bar{E}_{\gamma}^{(r, \bullet)} = \text{Med} \left( \left\{ \left| \hat{\gamma}^{(\bullet)}(x_i) - \gamma(x_i) \right| : i = 1, \dots, 2000 \right\} \right), $$
be a measure of the error in the estimation of the conditional extreme value index. Let also

$$ \bar{E}_{B}^{(r, \bullet)} = \| P_{\hat{B}_{n}^{(\bullet)}} - P_{B} \|_{F}, $$

be the Frobenius norm of the difference between the projection operators $P_M = MM^{\top}$ with $M \in \{B, \hat{B}_{n}^{(\bullet)}\}$.

In order to illustrate the influence of the estimation of the extreme value index and the SDR subspace on the accuracy of the conditional extreme quantile estimators, we represent the scatter plots

$$ \{ (\bar{E}_{\gamma}^{(r, \bullet)}, \bar{E}^{(r, \bullet)}) ; r = 1, \dots, R \} \text{ and } \{ (\bar{E}_{B}^{(r, \bullet)}, \bar{E}^{(r, \bullet)}) ; r = 1, \dots, R \}, $$

for each experiment and each method •. To avoid overloading the main part of the paper, these plots are gathered in the appendix. For each experiment, the heavy-tailed methods (HNN, HRF, TIREX, CTI) display a non-negligible relationship between the quality of the $\gamma(x)$ estimation and the overall quantile estimation performance. Surprisingly, the NN estimator of the extreme quantile does not seem to be affected at all by the accuracy of the extreme value index estimation, and only a weak correlation is observed for the RF method when applied to Fréchet or Burr distributions. This behavior may be explained by a compensation effect in the joint estimation of the scale and shape parameters. Regarding the estimation of the SDR subspace, it clearly appears – as expected – that the better the estimation of B, the better the conditional extreme quantile estimation obtained by the TIREX and CTI methods. It should also be noted that, in most cases, the estimation of B provided by the CTI approach is more accurate.

**Experiment 1** We assess the sensitivity of the estimators to the choice of the classical quantile level $\alpha_n$, which determines the number of excesses used in the tail parameter estimation. We consider three different values for this level, namely $\alpha_n \in \{0.05, 0.1, 0.2\}$. The dimension of the covariate is fixed to $p = 10$ and we set the correlation parameter $\nu = 0$, meaning that the covariates are independent. Two different conditional distributions are considered:

* **Fréchet:** $q(u, \gamma(x)) = (-\log(1 - u))^{-\gamma(x)}$,

* **Absolute student:** $q(u, \gamma(x)) = t_{1/\gamma(x)}(1 - u/2)$, where $t_\nu$ denotes the quantile function of a Student distribution with $\nu$ degrees of freedom.

As in Pasche and Engelke [21], the conditional tail index is given for $x = (x_1, \dots, x_p)^{\top}$ by

$$ \gamma(x) = \left( \frac{7}{1 + \exp(4x_1 + 1.2)} + 3 \right)^{-1} \in [0.1, 0.3]. \quad (21) $$

The results are reported in Figure 2. For the Fréchet conditional distribution, all methods are relatively robust to the choice of $\alpha_n$, with the exception of
HNN, which exhibits increased variability at small values of $\alpha_n$. In this setting, heavy-tailed variants, HRF and HNN, generally outperform GPD-based methods RF and NN, except HNN for $\alpha_n = 0.05$. In particular, the random forest method RF is particularly unstable.

For the Student conditional distribution, the GPD-based machine learning approaches RF and NN outperform their heavy-tailed counterparts. The dimension reduction approaches (TIREX and CTI), are highly sensitive to the choice of $\alpha_n$. For these methods to perform well, a small value of $\alpha_n$ seems essential.

<figure><img src="image_4.png" /><figcaption>Figure 2: **Experiment 1** – Boxplots of the symmetric relative error for the Fréchet and Absolute student distributions for different values of $\alpha_n$.</figcaption></figure>

**Experiment 2** We now evaluate the impact of the covariate dimension $p$. We take the setup of the first experiment but fix the intermediate conditional quantile level at $\alpha_n = 0.1$, and consider $p \in \{5, 10, 30, 200\}$. Note that the effective number of observations used for tail estimation is of order $n\alpha_n = 200$. Hence, when $p = 200$, the experiment corresponds to a high-dimensional setting where the number of covariates is comparable to the effective sample size. In this case, the CTI method, which relies on solving an optimization problem, becomes computationally prohibitive and is therefore not included in this configuration. The results, depicted in Figure 3, indicate that, except for the HNN method, the dimension of the covariate does not strongly influence the performance of the conditional quantile estimator. The neural networks approach HNN, which is designed for the heavy-tailed case, performs better than method NN when $p = 5$
but strongly deteriorates when $p$ increases. For the Fréchet distribution, the best results are clearly obtained by the approaches HRF and CTI for moderate dimensions ($p \le 30$).

<figure><img src="image_5.png" /><figcaption>Figure 3: **Experiment 2** – Boxplots of the symmetric relative error for the Fréchet and Absolute student distributions for different values of $p$.</figcaption></figure>

**Experiment 3** As mentioned in the introduction, the second-order parameter function $\rho < 0$ quantifies the rate of convergence in the limit result of model (M). We examine its influence by considering, in addition to Fréchet and student distributions, the Burr conditional distribution:

$$q(u, \gamma(x)) = (u^{\rho} - 1)^{-\gamma(x)/\rho}, \text{ with } \rho \in \{-0.2, -0.5, -1, -\infty\}.$$

Here, the second-order parameter function is equal to $\rho$ and does not depend on $x$. Note also that

$$\lim_{\rho \to -\infty} (u^{\rho} - 1)^{-\gamma(x)/\rho} = u^{-\gamma(x)}.$$

Hence, when $\rho = -\infty$, the Burr distribution coincides with a standard Pareto distribution. For the Fréchet distribution, $\rho = -1$ for all $x$, whereas for the student distribution, $\rho$ depends on $x$ and ranges approximately between $-0.2$ and $-2/3$. In this experiment, we take $\alpha_n = 0.1$, $\nu = 0$, $p = 10$.
The boxplots in Figure 4 highlight the effect of $\rho(x)$ on estimation accuracy. These results suggest that the second-order parameter plays a decisive role in determining the relative performance of the estimators. In general, when $\rho(x)$ is far from zero (i.e., fast convergence), heavy-tailed approaches tend to outperform GPD-based ones. Conversely, when $\rho(x)$ is close to zero (i.e., convergence is slower), the GPD methods yield more accurate estimates.

<figure><img src="image_6.png" /><figcaption>Figure 4: **Experiment 3** – Boxplots of the symmetric relative error for the Fréchet, Absolute student and Burr distributions.</figcaption></figure>

**Experiment 4** In this experiment, we investigate how the complexity of the conditional tail index $\gamma$ affects the performance of the estimators. To this end, we take $\alpha_n = 0.1$, $\nu = 0$, $p = 10$ and consider a Fréchet distribution with three different configurations for $\gamma$, each reflecting increasing levels of structural complexity:

$$ 
\begin{align*}
\gamma_1(x) &= 0.3 + 0.5 \mathbb{I}\{|x_1| > 1/2\}, \\
\gamma_2(x) &= 0.3 + 0.5 \mathbb{I}\{|x_1| + |x_2| > 1\}, \\
\gamma_3(x) &= 0.3 + 0.5 \mathbb{I}\{\|x\|_1 > p/2\}.
\end{align*}
$$

The conditional tail index $\gamma_1$ only depends on the first component of the covariate. In this case, it can be shown that $S = \text{span}(B)$ with $B = (1, 0, \dots, 0)^\top$ is the CTI subspace as introduced in Definition 7. The function $\gamma_2$ depends on the first two components of the covariate so that the CTI subspace is spanned
by the vectors $e_1 := (1, 0, \dots, 0)^\top$ and $e_2 := (0, 1, 0, \dots, 0)^\top$. Note that on the non linear set $\{(x, y); |x| + |y| > 1\}$, the function $\gamma_2$ is equal to 0.8 and to 0.3 on $\{(x, y); |x| + |y| \le 1\}$. Finally, for the function $\gamma_3$, the CTI is the set $\mathbb{R}^p$ and the function depends on the whole covariate.

The results are depicted in Figure 5. First remark that all methods perform poorly when the conditional tail index is $\gamma_3$. This can be explained by the complexity of the pattern that cannot be effectively captured, even by flexible machine learning-based approaches. For the function $\gamma_1$, the best results are obtained by the HRF and CTI methods. Finally, for $\gamma_2$, HRF appears to benefit from its ability to detect nonlinear structures, in contrast to CTI, which relies on linear projections.

<figure><img src="image_7.png" /><figcaption>Figure 5: **Experiment 4** – Boxplots of the symmetric relative error for the Fréchet distribution with three different tail index functions $\gamma_1$, $\gamma_2$, and $\gamma_3$. The experiment illustrates the impact of structural complexity in the conditional tail index.</figcaption></figure>

**Experiment 5** In this last experiment, we evaluate the robustness of the methods to correlation between covariates. Recall that the covariate vectors $X \in [0, 1]^p$ are generated from a Gaussian copula with correlation matrix $\Sigma_{ij} = \nu^{|i-j|}$, which induces an AR(1)-type dependence structure across coordinates while preserving uniform marginal distributions. We consider three values of the correlation parameter, $\nu \in \{0, 0.5, 0.8\}$, corresponding respectively to independent, moderately correlated, and strongly correlated covariates. This setting allows us to assess how the performance of the different estimators is affected by increasing multicollinearity in the design. For the other settings, they are identical to those in Experiment 1, and we fixe $\alpha_n = 0.1$ and $p = 10$.
<figure><img src="image_8.png" /><figcaption>Figure 6: **Experiment 5** — Boxplots of the symmetric relative error for the Fréchet and Absolute student distributions for different values of $\nu$.</figcaption></figure>

The results are shown in Figure 6. For the Fréchet distribution, all estimators exhibit good robustness to the introduction of correlation with only a slight degradation in performance as the correlation parameter $\nu$ increases. By contrast, for the Student distribution, the presence of correlation has a more pronounced effect. While the neural networks approach (NN) method remains remarkably stable across all correlation levels, the performance of the other approaches deteriorates as $\nu$ increases.

**Conclusion** Across the five experiments conducted in this study, we evaluated the performance of several recently proposed methods for extreme conditional quantile estimation under a variety of conditions. Our findings provide several practical insights for users seeking to apply these methods.

First, HNN systematically underperforms and should generally be avoided. In contrast, NN exhibits moderate performance but proves robust across all scenarios, making it a reliable and safe choice when no prior information is available about the tail structure.

Among all the approaches, RF is the better choice when the second-order parameter function $\rho$ is close to zero, corresponding to low convergence in the extreme value limit. However, it becomes highly unstable as $\rho$ moves away from zero. On the other hand, HRF is more stable and outperforms RF in scenarios where $\rho(x) < -0.5$.

Regarding dimension reduction methods, CTI is better or at least comparable to TIREX in most settings. When the tail index is defined on a lower dimensional linear subspace, CTI also outperforms all the tested methods as soon as $\rho(x) < -0.5$.
The experiment on correlated covariates reveals that correlation can lead to heterogeneous behavior across estimators and settings. Method NN exhibiting remarkable robustness while high correlations tend to distort the estimation of the other methods.

More broadly, our results show that GPD-based methods are generally less sensitive to the choice of the intermediate level $\alpha_n$ than their heavy-tailed counterparts. For heavy-tailed approaches, using a relatively small value such as $\alpha_n = 0.05$ or $\alpha_n = 0.1$ is recommended to ensure reliable estimation.

These findings offer practical guidance for selecting appropriate methods based on the underlying model characteristics and emphasize the importance of the unknown second-order function and of the structural complexity of the conditional tail index when estimating extreme quantiles.

## References

[1] Aggarwal, C. C. (2018). *Neural networks and deep learning: a textbook*. Springer.

[2] Aghbalou, A., Portier, F., Sabourin, A., & Zhou, C. (2024). Tail Inverse Regression: dimension reduction for prediction of extremes, *Bernoulli*, **30**(1), 503–533.

[3] Allouche, M., Girard, S., & Gobet, E. (2022). EV-GAN: Simulation of extreme events with ReLU neural networks, *Journal of Machine Learning Research*, **23**(150), 1–39.

[4] Allouche, M., Girard, S., & Gobet, E. (2024). Estimation of extreme quantiles from heavy-tailed distributions with neural networks, *Statistics and Computing*, **34**(1), 12.

[5] Athey, S., Tibshirani, J., & Wager, S. (2019). Generalized random forests, *The Annals of Statistics*, **47**(2), 1148–1178.

[6] Bhattacharya, P. K., & Gangopadhyay, A. K. (1990). Kernel and nearest-neighbor estimation of a conditional quantile, *The Annals of Statistics*, **18**(3)1400–1415.

[7] Bellman, R. (1961). *Adaptive Control Processes: A Guided Tour*. Princeton University Press.

[8] Berlinet, A., Gannoun, A., & Atzner-Løber, E. (2001). Asymptotic normality of convergent estimates of conditional quantiles, *Statistics*, **35**(2), 139–169.

[9] Breiman, L. (2001). Random forests, *Machine Learning*, **45**, 5–32.
[10] Carreau, J., & Girard, S. (2011). Spatial extreme quantile estimation using a weighted log-likelihood approach, *Journal de la Société Française de Statistique*, **152**(3), 66–82.

[11] Cook, R. D., & Weisberg, S. (1991). Sliced inverse regression for dimension reduction: Comment, *Journal of the American Statistical Association*, **86**(414), 328–332.

[12] Gardes, L. (2018). Tail Dimension Reduction for extreme quantile estimation, *Extremes*, **21**(1), 57–95.

[13] Gardes, L., & Podgorny, A. (2025). Dimension reduction for the estimation of the conditional tail-index, *Scandinavian Journal of Statistics*, **52**(3), 1444–1476.

[14] Gardes, L., & Stupfler, G. (2014). Estimation of the conditional tail index using a smoothed local Hill estimator, *Extremes*, **17**(1), 45–75.

[15] Gneco, N., Terefe, E. M., & Engelke, S. (2024). Extremal random forests, *Journal of the American Statistical Association*, **119**(548), 3059–3072.

[16] de Haan, L., & Ferreira, A. (2006). *Extreme Value Theory: An Introduction*, Springer, New York.

[17] Hill, B. M. (1975). A simple general approach to inference about the tail of a distribution, *The Annals of Statistics*, **3**(5), 1163–1174.

[18] Koenker, R. & Bassett Jr, G. (1978). Regression quantiles, *Econometrica*, **46**(1), 33–50.

[19] Li, K-C. (1991). Sliced inverse regression for dimension reduction, *Journal of the American Statistical Association*, **86**(414), 316–327.

[20] Li, Q., & Racine, J. S. (2008). Nonparametric estimation of conditional CDF and quantile functions with mixed categorical and continuous data, *Journal of Business and Economic Statistics*, **26**(4), 423–434.

[21] Pasche, O. C., & Engelke, S. (2024). Neural networks for extreme quantile regression with an application to forecasting of flood risk, *The Annals of Applied Statistics*, **18**(4), 2818–2839.

[22] Pickands III, J. (1975). Statistical inference using extreme order statistics, *The Annals of Statistics*, **3**(1), 119–131.

[23] Samanta, M. (1989). Non-parametric estimation of conditional quantiles, *Statistics and Probability Letters*, **7**(5), 407–412.

[24] Velthoen, J., Dombry, C., Cai, J. J., & Engelke, S. (2023). Gradient boosting for extreme quantile regression, *Extremes*, **26**(4), 639–667.
[25] Youngman, B. D. (2019). Generalized additive models for exceedances of high thresholds with an application to return level estimation for US wind gusts, *Journal of the American Statistical Association*, **114**(528), 1865–1879.

[26] Yu, K., & Moyeed, R. A. (2001). Bayesian quantile regression, *Statistics and Probability Letters*, **54**(4), 437–447.

## Appendix – Additional plots

This appendix gathers the scatter plots presented in Section 4.3. These plots help disentangle the respective contributions of extreme value index estimation and subspace recovery to the overall performance of the extreme conditional quantile estimators. Figures 7–11 display the scatter plots

$$ \{ (\bar{E}_{\gamma}^{(r, \bullet)}, \bar{E}^{(r, \bullet)}) ; r = 1, \dots, R \}, $$

for each experiment and method. The scatter plots

$$ \{ (\bar{E}_{B}^{(r, \bullet)}, \bar{E}^{(r, \bullet)}) ; r = 1, \dots, R \}, $$

are depicted in Figures 12–16.
<figure><img src="image_9.png" /><figcaption>Figure 7: Scatter plots of the errors of the conditional extreme quantile estimators versus the errors of the extreme value index estimators for Experiment 1 (effect of $\alpha_n$).</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>Figure 8: Scatter plots of the errors of the conditional extreme quantile estimators versus the errors of the extreme value index estimators for Experiment 2 (effect of $p$).</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 9: Scatter plots of the errors of the conditional extreme quantile estimators versus the errors of the extreme value index estimators for Experiment 3 (effect of $\rho$).</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>Figure 10: Scatter plots of the errors of the conditional extreme quantile estimators versus the errors of the extreme value index estimators for Experiment 4 (complexity of $\gamma(x)$).</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>Figure 11: Scatter plots of the errors of the conditional extreme quantile estimators versus the errors of the extreme value index estimators for Experiment 5 (correlated covariates).</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 12: Scatter plots of the errors of the conditional extreme quantile estimators versus the errors of the SDR subspace estimators for Experiment 1 (effect of $\alpha_n$).</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13: Scatter plots of the errors of the conditional extreme quantile estimators versus the errors of the SDR subspace estimators for Experiment 2 (effect of the dimension $p$).</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Figure 14: Scatter plots of the errors of the conditional extreme quantile estimators versus the errors of the SDR subspace estimators for Experiment 3 (effect of $\rho$).</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 15: Scatter plots of the errors of the conditional extreme quantile estimators versus the errors of the SDR subspace estimators for Experiment 4 (complexity of $\gamma(x)$).</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 16: Scatter plots of the errors of the conditional extreme quantile estimators versus the errors of the SDR subspace estimators for Experiment 5 (correlated covariates).</figcaption></figure>