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

HAL
open science

# Extreme-PLS with missing data under weak dependence

Stéphane Girard, Cambyse Pakzad

► To cite this version:

Stéphane Girard, Cambyse Pakzad. Extreme-PLS with missing data under weak dependence. 2025. [https://inria.hal.science/hal-05361239v1](https://inria.hal.science/hal-05361239v1).

HAL Id: hal-05361239

[https://inria.hal.science/hal-05361239v1](https://inria.hal.science/hal-05361239v1)

Preprint submitted on 12 Nov 2025

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 Attribution 4.0 International License
EXTREME-PLS WITH MISSING DATA UNDER WEAK DEPENDENCE

STÉPHANE GIRARD¹ AND CAMBYSE PAKZAD²

¹Univ. Grenoble Alpes, Inria, CNRS, Grenoble INP, LJK, 38000 Grenoble, France

²Modal'X, Université Paris Nanterre, 200 Avenue de la République, 92001 Nanterre, France

**ABSTRACT.** This paper develops a theoretical framework for *Extreme Partial Least Squares* (EPLS) dimension reduction in the presence of missing data and weak temporal dependence. Building upon the recent EPLS methodology for modeling extremal dependence between a response variable and high-dimensional covariates, we extend the approach to more realistic data settings where both serial correlation and missing-ness occur. Specifically, we consider a single-index inverse regression model under heavy-tailed conditions and introduce a *Missing-at-Random* (MAR) mechanism acting on the covariates, whose probability depends on the extremeness of the response. The asymptotic behavior of the proposed estimator is established within an $\alpha$-mixing framework, leading to consistency results under regularly varying tails. Extensive Monte-Carlo experiments covering eleven dependence schemes (including ARMA, GARCH, and nonlinear ESTAR processes) demonstrate that the method performs robustly across a wide range of heavy-tailed and dependent scenarios, even when substantial portions of data are missing. A real-world application to environmental data further confirms the method's capacity to recover meaningful tail directions.

**KEYWORDS:** Extreme Value Theory, Partial Least Squares, Missing Data, Weak Dependence, Dimension Reduction, Heavy Tails.

# 1. INTRODUCTION

In modern statistical applications, analysts increasingly face interconnected challenges that defy classical assumptions. Indeed, the data collected are often complex, sometimes arising as time series exhibiting serial dependence, with numerous predictors/covariates, which weights on statistical procedures and their computational cost. Furthermore and realistically, datasets often suffer from missing data, particularly in the extreme values of the response variable – the very observations often most critical to study. Together, these issues build up a compounded sparsity with exacerbated curse of dimensionality by the rarity of extreme events, which in turn may be under-recorded because of the missing-ness mechanism. This convergence undermines standard statistical methods and creates a pressing need for tailored robust techniques.

Addressing the challenge of raw sparsity due to high-dimensions, the field of sufficient dimension reduction (SDR) seeks a low-dimensional projection $P_S X$ of the high-dimensional covariate $X$ preserving the information about the response variable $Y$. More precisely, it is such that $Y$ conditionally to $X$ has same distribution as $Y$ conditionally to $P_S X$, or equivalently $Y$ and $X$ are independent conditionally to $P_S X$ (Cook, 2007). A classic example is the single-index model, which assumes $\dim(S) = 1$ and therefore $Y$ only depends on a single linear combination of $X$ with a deterministic vector $\beta$. Next, one would substitute $X$ by $P_S X$ in further statistical methods as it carries the same amount of information on $Y|X$ without the dimensionality curse.

E-mail address: ¹stephane.girard@inria.fr, ²cambyse.pakzad@parisnanterre.fr.
Two prominent SDR methods are Partial Least Squares (PLS), initiated by Wold (1975), and Sliced Inverse Regression (SIR), introduced by Li (1991). PLS finds projections with high variance and high covariance with $Y$, meaning that the information in $X$ that impacts $Y$ is encompassed in $P_X X$, while being the most exhaustive. SIR, in contrast, leverages the simpler inverse regression 'X against Y', usually at the price of a certain linearity condition and constant variance. This inverse model, where the inference is rather conditioned on the predictors/covariates $X$, somewhat contradicts the philosophy of Fisher, as noted but also encouraged by Cook (2007).

A nascent research area aims to merge SDR with extreme-value theory since nonparametric estimators of extreme conditional features (Girard et al., 2021) may be impacted both by the scarcity of extremes and the high-dimensional setting. The challenging task is to capture the information relevant specifically to the response tails. Pioneering work by Gardes (2018) introduced a notion of tail conditional independence. Later, Xu et al. (2022) and Aghbalou et al. (2024) developed estimators for extreme quantiles and tail inverse regression, respectively. However, these approaches largely assume independent data and a fully observed response, conditions that are often violated in practice.

This paper bridges this gap by establishing the theoretical foundations for dimension reduction in the extremes under more realistic conditions: serial dependence and missing data. Our focus is on Extreme-PLS (EPLS, Bousebata et al. (2023)) relying on PLS principles for estimating the linear combination $\beta^\top X$ of $X$ that best explains the extreme values of $Y$. The asymptotic properties of the estimated direction $\hat{\beta}$ are established assuming a single-index inverse regression model and extreme-value type conditions. Here, the asymptotic properties of the estimator are rather derived under the more realistic context of dependent and missing data. Specifically, a Missing At Random (MAR, Rubin (1976)) mechanism is assumed on the covariate $X$, the missingness probability depending on the extreme-ness of the response variable $Y$. This, for instance, captures scenarios where the failure of a sensor (leading to a missing covariate) is more likely during an extreme weather event. We refer to Beirlant et al. (2023); Xu et al. (2022) for recent works on the estimation of tail features in the case where the missingness mechanism is assumed on the response variable $Y$. Finally, an $\alpha$-mixing framework is adopted both for the observations $(X, Y)$ and the masking process associated with the missing data, in order to account for the temporal dependence common in real-world data like financial or environmental time series.

The remainder of the paper is organized as follows. Some links between EPLS and the Marginal Expected Shortfall (MES, Cai et al. (2015)) risk measure are established in Section 2. The missing data and dependence frameworks are described in Section 3 and the EPLS estimator is derived. Its asymptotic properties are established in Section 4 under tail conditions while its finite sample behavior is illustrated on large scale experiments in Section 5 and Section 6. The first part of the experiments is conducted on simulated data involving eleven different dependence schemes. The second part illustrates the performance of the method compared to four competitor baselines on environmental data involving measures at nearly 5800 triplets of weather stations. Proofs are postponed to the Appendix.

## 2. EPLS AND MARGINAL EXPECTED SHORTFALL

Let us consider the vector $w(y)$ referred to as the EPLS direction in Bousebata et al. (2023). It is defined as the unit vector maximizing over $||w|| = 1$ the covariance between $Y$ and the projection
$w^\top X$ of $X$ on $w$ given that $Y$ exceeds a large threshold $y > 0$:

$$w(y) = \arg \max_{\|w\|=1} \operatorname{cov}(w^\top X, Y | Y \ge y), \qquad (1)$$

where $\| \cdot \|$ denotes the Euclidean norm. We also introduce the marginal expected shortfall of $X$ at level $y$ which plays a central role in this study:

$$\mathrm{MES}_X(y) = \mathbb{E}(X | Y > y),$$

see Cai et al. (2015). Note that $\mathrm{MES}_X(y) \in \mathbb{R}^p$ and the associated jth coordinate is given by $\mathrm{MES}_X^{(j)}(y) = \mathrm{MES}_{X^{(j)}}(y)$. This risk measure may alternatively be expressed thanks to the tail-moment defined as $m_X(y) = \mathbb{E}(X1_{\{Y \ge y\}})$, whenever it exists. This notation will be kept throughout our work. In fact, denoting $\bar{F}(y) = \mathbb{P}(Y \ge y)$ the survival function of $Y$, one may write $\mathrm{MES}_X(y) = m_X(y)/\bar{F}(y)$. Let us also remark that

$$\mathrm{MES}_Y(y) = \mathbb{E}(Y | Y > y) := \mathrm{ES}_Y(y)$$

is the so-called Expected Shortfall (ES). The linear optimization problem (1) under a quadratic constraint benefits from a closed-form solution obtained with Lagrange multipliers method and involving the cross-MES between $X$ and $Y$. The next proposition precises this point, it is a simple rewriting of (Bousebata et al., 2023, Proposition 1).

**Proposition 1.** Suppose that $\mathbb{E}(\|X\|1_{\{Y \ge y\}}) < \infty$, $\mathbb{E}(|Y|1_{\{Y \ge y\}}) < \infty$ and $\mathbb{E}(\|XY\|1_{\{Y \ge y\}}) < \infty$ for all $y > 0$. Then, the unique solution of the optimization problem (1) is given for all $y > 0$ by:

$$w(y) = \frac{v(y)}{\|v(y)\|}, \quad \text{where } v(y) := \mathrm{MES}_{XY}(y) - \mathrm{MES}_X(y)\mathrm{ES}_Y(y). \qquad (2)$$

In the following, we aim at investigating the behavior of $w(y)$ for large thresholds $y$. without resorting neither to a linear conditional expectation assumption as in Xu et al. (2022) nor to a conditional independence assumption as in Gardes (2018); Saracco (1997). In contrast, additional assumptions on the distribution tails are introduced below.

### 3. STATISTICAL FRAMEWORK

Here and throughout, $\odot$ and $/$ are respectively the componentwise product and ratio. For a generic random vector $Z \in \mathbb{R}^p$, let us recall that $Z^{(j)}$ denotes the jth coordinate of $Z$ and the jth coordinate of $m_Z$ is then given by $m_Z^{(j)} = m_{Z^{(j)}}$, $j \in \{1, \dots, p\}$.

**Model and tail conditions.** Let us consider the single-index inverse regression model introduced in Bousebata et al. (2023):

(M1) : $X = g(Y)\beta + \varepsilon$ where $X$ and $\varepsilon$ are $p$-dimensional random vectors, $g: \mathbb{R} \to \mathbb{R}$ is an unknown link function and $\beta \in \mathbb{R}^p$ is the unknown direction of interest.

Model (M1) is an inverse regression model since the covariates $X = (X^{(1)}, \dots, X^{(p)})^\top$ are written as functions of the response variable $Y$. A Bayesian version of EPLS is introduced in Arbel et al. (2024) basing on (M1) while similar models were used to establish the theoretical properties of SIR, see for instance Bernard-Michel et al. (2009); Cook (2007).

Three assumptions on the link function $g$ and the distribution tail of $Y$ and $\varepsilon$ are considered. They rely on the notion of regularly-varying functions. Recall that $\varphi \in RV_\theta$ with $\theta \in \mathbb{R}$ if and only if $\varphi$ is positive and

$$\lim_{y \to \infty} \frac{\varphi(ty)}{\varphi(y)} = t^\theta,$$
for all $t > 0$. This property is denoted for short by $\varphi \in RV_{\theta}$. We refer to Bingham et al. (1987) for a detailed account on regular variations.

(A₁) : Y is a real random variable with density function $f \in RV_{-1/\gamma-1}$ with $\gamma \in (0, 1)$.

(A₂) : $g \in RV_{\kappa}$ with $\kappa > 0$.

(A₃) : There exists $q \in (1, +\infty)$ such that $E(\|\varepsilon\|^q) < \infty$.

Assumption (A₁) implies that the distribution of Y is in the Fréchet maximum domain of attraction with positive tail-index $\gamma$, see (Bingham et al., 1987, Theorem 1.5.8) and (de Haan and Ferreira, 2007, Theorem 1.2.1). This domain of attraction consists of heavy-tailed distributions, such as Pareto, Student t or Burr distributions, see Beirlant et al. (2004) for other examples. The larger $\gamma$ is, the heavier the tail. The restriction to $\gamma < 1$ ensures that the first-order moment $E(|Y|1_{\{Y \ge y\}})$ exists for all $y > 0$. Assumption (A₂) states that $g$ ultimately behaves like a power function. Assumption (A₃) can be interpreted as an assumption on the tails of $\|\varepsilon\|$. It is satisfied, for instance, by distributions with exponential-like tails such as Gaussian, Gamma or Weibull distributions.

**Missig data framework.** We consider the situation where some coordinates of the covariate X may not be observed. A random mask $\Lambda \in \{0, 1\}^p$ is applied on the covariate, such that only the $\Lambda \odot X$ part of X is available for estimation. The missing at random (Rubin, 1976) mechanism is adopted here, i.e. $\Lambda$ is dependent of Y but not of X:

(A₄) : $\Lambda$ is a $\{0, 1\}^p$ random vector such that, for all $j \in \{1, \dots, p\}$, $\Lambda^{(j)} | Y$ are independent Bernoulli random variables with success parameter $\lambda_j(Y) = c_j\lambda(Y)$ with $\lambda \in RV_{\tau}$, $\tau < 0$ and $c_j > 0$.

Note that the rate of non-missing data may depend on the coordinate $j$ through the multiplicative constant $c_j > 0$ and decreases in the distribution tail as a power function with exponent $\tau < 0$. We refer to Mayer et al. (2022) for more information on missing data. Let $(X_i, Y_i, \Lambda_i, \varepsilon_i)_{1 \le i \le n}$ be a random sequence sampled from a random vector $(X, Y, \Lambda, \varepsilon)$ under model (M₁). Some vectors $X_i$ may not be actually observed because of the missing-ness and our true observations set is thus $(Y_i, \Lambda_i \odot X_i)_{1 \le i \le n}$.

**Serial dependence.** We adopt the $\alpha$-mixing framework for its generality, see (Doukhan, 1994, Proposition 1 and diagram page 20). For any two positive integers $a \le b \le +\infty$, let $\mathcal{F}_a^b$ be the $\sigma$-algebra generated by $\{(X_i, Y_i, \Lambda_i, \varepsilon_i)\}_{a \le i \le b}$. The sequence $\{(X_i, Y_i, \Lambda_i, \varepsilon_i)\}_{1 \le i \le n}$ is said $\alpha$-mixing if $\alpha(n) \to 0$ as $n \to \infty$ where

$$ \alpha(n) := \sup_{k \ge 1} \sup_{A \in \mathcal{F}_1^k} \sup_{B \in \mathcal{F}_{k+n}^\infty} |\mathbb{P}(A \cap B) - \mathbb{P}(A)\mathbb{P}(B)|. $$

To setup the standard 'small-block/large-block' argument (see Drees and Rootzén (2010)), the following assumptions are introduced:

(M₂) : $(X_i, Y_i, \Lambda_i, \varepsilon_i)_{1 \le i \le n}$ is a strictly-stationary $\alpha$-mixing time serie, sampled from the random vector $(X, Y, \Lambda, \varepsilon)$ under model (M₁), and there exist two sequences of integers $1 \ll \ell_n \ll r_n \ll n$ and $n\alpha(\ell_n) \ll r_n$, such that

$$ \sum_{j \ge 1} \alpha^{\frac{\delta}{2+\delta}}(j) < \infty \text{ for some } \delta > 0. \quad (3) $$

While conditions $1 \ll \ell_n \ll r_n \ll n$ and $n\alpha(\ell_n) \ll r_n$ are at the basis of the 'small-block/large-block' argument, assumption (3) is standard in the literature of central limit theorems for strongly mixing processes, see Ibragimov (1962, Theorem 1.7).
**Inference.** Let $y_n \to \infty$ as the sample size $n$ tends to infinity. The EPLS with missing data (2) is estimated by the random direction $\hat{\beta}_{\Lambda}(y_n) := \hat{v}_{\Lambda}(y_n)/\|\hat{v}_{\Lambda}(y_n)\| \in \mathbb{R}^p$, based on

$$ \hat{v}_{\Lambda}(y_n) := \frac{1}{\hat{m}_{\Lambda}(y_n)} \left( \hat{F}(y_n) \hat{m}_{Y\Lambda \odot X}(y_n) - \hat{m}_Y(y_n) \hat{m}_{\Lambda \odot X}(y_n) \right) \quad (4) $$

the ratio being understood component-wise with, for any random vector $Z \in \mathbb{R}^p$,

$$ \hat{m}_Z(y_n) := \frac{1}{n} \sum_{i=1}^{n} Z_i \mathbb{1}_{\{Y_i \ge y_n\}}, $$

and where $\hat{F} := \hat{m}_1$ is the empirical survival function.

## 4. ASYMPTOTIC PROPERTIES OF THE ESTIMATOR

We first provide a tool establishing the joint asymptotic properties of $\hat{m}_{\Lambda \odot X}(y_n)$, $\hat{m}_{\Lambda}(y_n)$, $\hat{m}_{Y\Lambda \odot X}(y_n)$ and $\hat{m}_Y(y_n)$ when $y_n \to \infty$. This latter condition ensures that the rate of convergence, which is driven by the number of effective observations $n(\lambda \bar{F})^{\theta}(y_n)$ used in the estimators, tends to infinity as the sample size increases.

From now on, we denote by $d$ the number of non-zero components of $\beta$, supposed positive $d \ge 1$, namely $d := |\{j \in \{1, \dots, p\} : \beta^{(j)} \neq 0\}|$. Without loss of generality and for simplicity, one may assume the support of $\beta$ to be exactly $\{1, \dots, d\}$.

**Proposition 2.** Suppose $(M_2)$, $(A_1)$, $(A_2)$, $(A_3)$ and $(A_4)$ hold. Assume also that

$$ 1 < \gamma \min \{2\kappa + 2q + \tau; q\kappa + 2\tau; q\kappa/2 + \tau\}, $$

$$ \gamma \max \left\{ \frac{q}{q-1}\kappa + \tau; \frac{q}{q-1}(\kappa+1) + 2\tau; \frac{q}{q-2-\delta} + \tau; (2+\delta)\kappa + \tau + 1; 2\kappa + \tau + 2; \kappa+1; 2+\delta \right\} < 1. $$

Introduce $\Xi_n$ the $\mathbb{R}^{3d+2}$-random vector defined as:

$$ \left\{ \left( \frac{\hat{m}_{\Lambda^{(j)}X^{(j)}}(y_n)}{m_{\Lambda^{(j)}X^{(j)}}(y_n)} - 1 \right)_{1 \le j \le d}, \left( \frac{\hat{m}_{\Lambda^{(j)}(y_n)}}{m_{\Lambda^{(j)}(y_n)}} - 1 \right)_{1 \le j \le d}, \left( \frac{\hat{m}_{Y\Lambda^{(j)}X^{(j)}}(y_n)}{m_{Y\Lambda^{(j)}X^{(j)}}(y_n)} - 1 \right)_{1 \le j \le d}, \left( \frac{\hat{m}_{Y^i(y_n)}}{m_{Y^i(y_n)}} - 1 \right)_{i \in \{0,1\}} \right\}. $$

For any $\theta > 1 + \frac{\delta}{2+\delta} \in (1, 2)$ such that $n(\lambda \bar{F})^{\theta}(y_n) \gg 1$ holds, the sequence $n^{1/2}(\lambda \bar{F})^{\theta/2}(y_n) \Xi_n$ converges in probability to zero in $\mathbb{R}^{3d+2}$. When $\theta = 1 + \frac{\delta}{2+\delta}$, if furthermore

$$ \gamma \max \left\{ (2 + \delta)(\kappa + 1) + \tau; \frac{q(2 + \delta)}{q - (2 + \delta)} + \tau \right\} < 1, $$

then the sequence $n^{1/2}(\lambda \bar{F})^{\theta/2}(y_n) \Xi_n$ converges in distribution, at least up to a subsequence, to a centered Gaussian random vector in $\mathbb{R}^{3d+2}$.

Proposition 2 holds in the weaker $\psi, \phi, \rho$-mixing frameworks, see Doukhan (1994), with a better rate that corresponds to $\delta = 0$. Yet in any case, we found ourselves limited to only provide a non-explicit bound on the limiting covariance. Our main result establishes that $\hat{v}_{\Lambda}(y_n)$ and $\beta$ are asymptotically aligned in $\mathbb{R}^p$, with respect to the convergence in probability and rate $n^{1/2}(\lambda \bar{F})^{\theta/2}(y_n)$.

**Theorem 1.** Under the assumptions of Proposition 2, for any $\theta > 1 + \frac{\delta}{2+\delta}$, it holds that:

$$ \min \left\{ n^{1/2}(\lambda \bar{F})^{\theta/2}(y_n), g(y_n) \bar{F}^{1/q}(y_n) \right\} (\hat{\beta}_{\Lambda}(y_n) - \beta) \xrightarrow{n \to +\infty} 0, \quad \text{in } \mathbb{R}^p. $$

When $\theta = 1 + \frac{\delta}{2+\delta}$, if moreover $n^{1/2}(\lambda \bar{F})^{\theta/2}(y_n) \ll g(y_n) \bar{F}^{1/q}(y_n)$ and $ng^{-2}(y_n)(\lambda \bar{F})(y_n)^{\theta-2/q} \ll 1$, then, the following convergence holds, at least up to a subsequence,

$$ n^{1/2}(\lambda \bar{F})^{\theta/2}(y_n) (\hat{\beta}_{\Lambda}(y_n) - \beta) \xrightarrow{n \to +\infty} \beta \odot G, \quad \text{in } \mathbb{R}^p, $$
where G is some centered Gaussian vector in $\mathbb{R}^p$.

Note that, when the noise admits moments at any order, we may consider $q = +\infty$, thus boiling the conditions down to $\gamma \max\{(2+\delta)\kappa + \tau + 1; 2\kappa + \tau + 2; \kappa + 1; 2+\delta\} < 1$. Regarding the threshold regime associated with a finite value of $q$, let us consider $y_n := n^c$ with $c \in (0, 1)$. The property $n(\lambda \bar{F})^\theta(y_n) \gg 1$ is satisfied for a generic $\theta > 0$ whenever $c < \frac{\gamma}{\theta(1-\gamma\tau)}$. It is then checked through numerical computations that the set of all parameters $(q, \delta, \gamma, \kappa, \tau)$ satisfying the conditions needed for Theorem 1 is non-empty.

## 5. VALIDATION ON SIMULATED DATA

The large scale experiments are presented in Section 5.1, some implementation details are provided in Section 5.2 and the obtained results are discussed in Section 5.3.

5.1. **Experimental design.** The performance of our method is assessed by a Monte Carlo simulation experiment with $N = 500$ independent repetitions. The experimental design involves four setups for drawing a sample $(X_i, Y_i)_{1\le i\le n}$ with $n = 500$ and $p = 101$, see Paragraphes 5.1.1–5.1.4. In any case, the deterministic $\beta$'s coordinates are set to $\beta^{(j)} = \sqrt{2}\sin(2\pi j/p)$ for $1 \le j \le p$ and $g: y \in \mathbb{R}^+ \mapsto y^\kappa$ with $\kappa = 0.5$ fixed for shortness. To balance the contributions of signal and noise, we propose to rescale the conditional noise variance by multiplying $\varepsilon$ with $g(y)/10$. Once both the response variable $Y$ and the noise $\varepsilon$ are simulated, and the deterministic functions $\beta$ and $g$ are chosen, the covariate sampling $\{X_1, \dots, X_n\}$ is readily derived from (M$_1$). The mask dynamic framework is also common to all setups and is described in Paragraph 5.1.5. Finally, some remarks on the tail heaviness and the strict stationarity of the processes are provided in Paragraph 5.1.6 and Paragraph 5.1.7 respectively.

5.1.1. *First setup: univariate ARMA-resp./multivariate GARCH-noise.* The output and noise random variables are distributed as follows:

* $(Y_i)$ is a ARMA(1, 1)-time serie, i.e. $Y_i = \phi_{\text{resp}}Y_{i-1} + \theta_{\text{resp}}\eta_{i-1} + \eta_i$, where the innovations $\eta_i$ are i.i.d. BurrXII-distributed with parameters $-\rho$ and $1/\gamma$, i.e., with survival c.d.f. $\bar{F}(y) = (1 + y^{-\rho})^{1/\gamma}$ with $y \ge 0$, $\gamma \in (0, 1)$ being the tail-index and $\rho < 0$ being the second-order parameter. Here, $\gamma \in \{0.1, 0.4, 0.7\}$. Throughout the whole section, $\rho = -1$ is fixed. Since ARMA(1, 1) is a linear process, the response $(Y_i)$ is also heavy-tailed with same tail-index $-1/\gamma$, akin to the innovations $(\eta_i)$. The ARMA(1, 1) parameters $(\phi_{\text{resp}}, \theta_{\text{resp}})$ are picked among $\{(0.8, -0.3), (0.99, -0.98)\}$. The standard case $(\phi_{\text{resp}}, \theta_{\text{resp}}) = (0.8, -0.3)$ induces a stable and moderate mean-reverting process. In contrast, $(\phi_{\text{resp}}, \theta_{\text{resp}}) = (0.99, -0.98)$ creates a pathological dynamic where a single shock/innovation can trigger a lengthy period of volatility and extreme values, making events highly persistent.

* $(\varepsilon_i)$ conditionally to $Y = y$ is a GARCH(1, 1)-multivariate time serie, i.e., $\varepsilon_i = \sigma_i\eta'_i$ and $\sigma_i = \omega_{\text{noise}} + \alpha_{\text{noise}}\varepsilon_{i-1}^2 + \beta_{\text{noise}}\eta_{i-1}^2$ where the innovations $(\eta'_i)$ are i.i.d. centered Gaussian vectors in $\mathbb{R}^p$ with Toeplitz correlation matrix, i.e., $\text{Corr}(\varepsilon_i, \varepsilon_j | Y = y) = \rho_c^{|i-j|}$ with $\rho_c = 0.8$. Three sets of parameters are considered for the GARCH(1, 1) model of the noise: $(\omega_{\text{noise}}, \alpha_{\text{noise}}, \beta_{\text{noise}}) \in \{(0.05, 0.1, 0.85), (0.05, 0.05, 0.94)\}$. These configurations are chosen to capture a range of volatility dynamics. The standard case $(\omega_{\text{noise}}, \alpha_{\text{noise}}, \beta_{\text{noise}}) = (0.05, 0.1, 0.85)$ yields persistent but ultimately transitory volatility clustering, meaning that a period of high variance is likely to be followed by another, but the memory of this volatility decays over time. The parameters $(\omega_{\text{noise}}, \alpha_{\text{noise}}, \beta_{\text{noise}}) = (0.05, 0.05, 0.94)$ define
a process with extreme heteroskedasticity. Here, $\alpha_{noise} + \beta_{noise} = 0.99$ is close to one, and the process almost behaves like an Integrated GARCH (IGARCH), see (Engle and Bollerslev, 1986; Nelson, 1990), where volatility shocks become permanent.

5.1.2. *Second setup: univariate GARCH-resp./multivariate ARMA-noise.* Next, we move onto the second experimental setup by reversing the serial dependence structure of response and noise in the previous setup.

* The response ($Y_i$) is an univariate GARCH(1, 1)-process with i.i.d. BurrXII(-$\rho$, 1/$\gamma$)-distributed innovations. Concerning its parameters, we opt here for ($\omega_{resp}$, $\alpha_{resp}$, $\beta_{resp}$) = ($\omega_{noise}$, $\alpha_{noise}$, $\beta_{noise}$), meaning that they are the same as when the noise was GARCH(1, 1), while the innovations are to be understood as univariate i.i.d. BurrXII(-$\rho$, 1/$\gamma$)-distributed.

* The multivariate ARMA(1, 1)-structure on the noise has parameters corresponding to the ARMA-response case, i.e., ($\phi_{noise}$, $\theta_{noise}$) = ($\phi_{resp}$, $\theta_{resp}$). The innovations are similar to the GARCH-noise case, namely i.i.d. centered Gaussian vector in $\mathbb{R}^p$ with Toeplitz correlation matrix and $\rho_c = 0.8$.

5.1.3. *Third setup: ESTAR-resp./GARCH-noise.* We consider the Exponential Smooth Transition Autoregressive (ESTAR) model (van Dijk et al., 2002) for the response variable to model nonlinear dynamics. Here, the ESTAR(1) process is defined as:

$$Y_i = \phi_{\text{low}} Y_{i-1} + \phi_{\text{high}} Y_{i-1} [1 - \exp(-Y_{i-1}^2)] + \eta_i,$$

where $\eta_i \sim \text{BurrXII}(-\rho, 1/\gamma)$ are i.i.d. innovations and ($\phi_{\text{low}}$, $\phi_{\text{high}}$) = (0.2, 0.95). Endowed with these parameters, the dynamics smoothly shift from a quickly mean-reverting or low-persistence regime ($\phi_{\text{low}} = 0.2$) near equilibrium to a regime ($\phi_{\text{high}} = 0.95$) where large deviations tend to persist during extreme events. Concerning the noise, the process follows the multivariate GARCH(1, 1) specification from the first setup.

5.1.4. *Fourth setup: Independence.* The independent case is considered with the following combination of parameters:

$$(\phi_{\text{resp}}, \theta_{\text{resp}}) = (0, 0), \quad (\omega_{\text{noise}}, \alpha_{\text{noise}}, \beta_{\text{noise}}) = (1, 0, 0), \quad \alpha_{\text{bar}} = 0.$$

In this configuration, extremes occur sporadically and independently. The volatility is constant/clusterless, namely, an extreme value provides no information about the likelihood of another extreme value at a future time.

5.1.5. *Missing data mechanism.* The mask ($\Lambda_i$) is a multivariate Binary Autoregressive (or BAR) time serie, with marginal correction such that ($A_4$) is satisfied with $c_j = 1$ and $\lambda(y) = y^\tau$, $\tau \in \{-0.1, -0.5, -0.9\}$. Precisely, the BAR dynamics with cross-sectional independence writes as, for any $1 \le j \le p$ and $2 \le i \le n$,

$$\Lambda_1^{(j)} \sim \text{Bernoulli}(p_1), \quad \pi_{i,j} = \alpha_{\text{bar}}\Lambda_{i-1}^{(j)} + (1 - \alpha_{\text{bar}})p_i, \quad \text{and} \quad \Lambda_i^{(j)} \sim \text{Bernoulli}(\pi_{i,j}).$$

The parameters of the BAR process are $\alpha_{\text{bar}} \in [0, 1)$ which tunes the dependence and $p_i \in [0, 1]$ the target Bernoulli parameter for $\Lambda_i^{(j)}$. In our case, we fix $\alpha_{\text{bar}} = 0.5$ and $p_i = \lambda(Y_i)$ which corresponds to $c_j = 1$ in ($A_4$). Nonetheless, each $\Lambda_i^{(j)}$ though Bernoulli-distributed has its mean different from $p_i$. To remedy this hindrance and fit into our theoretical framework, we use a correction mechanism based on a rejection sampling technique. It starts with an initialization $\Lambda_1^{(j)} \sim \text{Bernoulli}(p_1)$. Next,
for any $i \ge 2$, we compute $\pi_{i,j}$ and independently sample $\tilde{\Lambda}_i^{(j)} \sim \text{Bernoulli}(\pi_{i,j})$ as candidate value and a uniform $U$. Then, one applies the asymmetric correction:

$$ \Lambda_i^{(j)} = \begin{cases} 1 & \text{if } \pi_{i,j} < p_i \text{ and } U < \gamma_{i,j}^{+} \text{ and } \tilde{\Lambda}_i^{(j)} = 0, \\ 0 & \text{if } \pi_{i,j} > p_i \text{ and } U < \gamma_{i,j}^{-} \text{ and } \tilde{\Lambda}_i^{(j)} = 1, \\ \tilde{\Lambda}_i^{(j)} & \text{otherwise.} \end{cases} $$

where $U \sim \text{Uniform}(0, 1)$ is independent of any source of randomness and the correction probabilities are $\gamma_{i,j}^{+} = \frac{p_i - \pi_{i,j}}{1 - \pi_{i,j}}$, $\gamma_{i,j}^{-} = \frac{\pi_{i,j} - p_i}{\pi_{i,j}}$. The idea is that if $\pi_{i,j} < p_i$ (under-generation of success), we flip the candidate value from 0 to 1 with probability $\gamma_{i,j}^{+}$ and keep it 1 if it is already the candidate value. We symmetrically proceed in the case of success over-generation. This ensures marginal consistency, i.e. $\mathbb{E}(\Lambda_i^{(j)}) = p_i$ for any $1 \le i \le n$ and $1 \le j \le p$. Indeed, one may condition with respect to $\pi_{i,j}$ and do a case analysis. Write $\mathbb{E}(\Lambda_i^{(j)}) = \mathbb{E}(\Lambda_i^{(j)} | \pi_{i,j})$. On the event $A = \{\pi_{i,j} < p_i\}$, we have by independence,

$$ \mathbb{E}(\Lambda_i^{(j)} | \pi_{i,j}) = \mathbb{P}(\tilde{\Lambda}_i^{(j)} = 1 | \pi_{i,j}) + \mathbb{P}(\tilde{\Lambda}_i^{(j)} = 0, U_i^{(j)} < \gamma_{i,j}^{+} | \pi_{i,j}) = \pi_{i,j} + (1 - \pi_{i,j}) \cdot \gamma_{i,j}^{+} = p_i. $$

On the event $B = \{\pi_{i,j} > p_i\}$, similarly $\mathbb{E}(\Lambda_i^{(j)} | \pi_{i,j}) = \mathbb{P}(\tilde{\Lambda}_i^{(j)} = 1, U_i^{(j)} \ge \gamma_{i,j}^{-} | \pi_{i,j}) = p_i$. On the event $C = \{\pi_{i,j} = p_i\}$, one has $\gamma_{i,j}^{+} = \gamma_{i,j}^{-} = 0$ which implies $\Lambda_i^{(j)} = \tilde{\Lambda}_i^{(j)}$ and thus, $\mathbb{E}(\Lambda_i^{(j)} | \pi_{i,j}) = \mathbb{E}(\tilde{\Lambda}_i^{(j)} | \pi_{i,j}) = \pi_{i,j} = p_i$. The claim follows by taking the expectation, since $\{A, B, C\}$ is a partition of $\Omega$.

### 5.1.6. Tail heaviness.
The GARCH(1, 1)-structure tends to increase the tail heaviness, even with Gaussian innovations. Thus, we consider in the second setup that $\gamma \in \{0.1, 0.4, 0.5\}$ with 0.5 instead of 0.7 in the first setup. The case $\gamma = 0.5$ is somewhat the limiting possibility if one wants to reach $\tau \le -0.5$. In fact, (Basrak et al., 2002, Remark 3.4) provides a way to compute the resulting tail-index $\gamma_g$ (with $g$ as GARCH) of the time series, by solving $\mathbb{E}((\alpha_{\text{resp}}\tilde{\eta}_1^2 + \beta_{\text{resp}})^{1/(2\gamma_g)}) = 1$ where $\tilde{\eta}_1$ is a standardized BurrXII(-$\rho$, $1/\gamma$)-random variable. Some Monte-Carlo computations yield, when $\gamma = 0.1$, in the standard GARCH(1, 1)-case a new tail-index $\gamma_g \approx 0.28$ and in the IGARCH-like case, $\gamma_g \approx 0.31$. When $\gamma = 0.4$, it respectively provides $\gamma_g \approx 0.43$ and $\gamma_g \approx 0.45$. In any case, the Monte-Carlo method fails when $\gamma = 0.5$.

### 5.1.7. Strict-stationarity and $\alpha$-mixing.
Denote momentarily $x \in \{\text{resp, noise}\}$. From Mokkadem (1988), the ARMA(1, 1) process with i.i.d. finite mean innovations (Gaussian or heavy-tailed) and $|\phi_x| < 1$ is strictly-stationary and $\alpha$-mixing. For the GARCH(1, 1)-process to be strictly-stationary, it is enough to verify by Monte-Carlo computations that $\mathbb{E}(\log(\alpha_x Z + \beta_x)) < 0$ when $Z \sim Y$ or $Z \sim \mathcal{N}(0, 1)$ (noise case). Thus, by Carrasco and Chen (2002) and since $\alpha_x + \beta_x < 1$, it is also $\alpha$-mixing. Nonetheless, the log-moment condition is not numerically true when $\gamma = 0.5$ (and even less when $\gamma = 0.7$), in the case of GARCH(1, 1)-response, indicating that the process fails to be strictly-stationary. Still, we decided to include this particularly pathological case as our method performs well. At last, the mask $(\Lambda_i)$ is a finite-state, irreducible, aperiodic Markov chains hence is $\alpha$-mixing, see (Bradley, 2005, Theorem 3.2).

### 5.2. Selection of the threshold.
The choice of the diverging threshold is a recurrent problem in extreme-value statistics, see Scarrott and MacDonald (2012) for a review. In the considered setup, a natural method is to select the value $y_n$ which maximizes the empirical covariance between $Y$ and $\hat{\beta}_{\Lambda}(y_n)^{\top}X$ using a grid search on $(0, \infty)$. This optimization problem can be simplified by remarking that the considered empirical covariance is a piecewise constant function of $y_n$ with $n$
jumps at $Y_1, \dots, Y_n$. It is thus sufficient to consider a discrete grid; to this end let us introduce $\bar{r}(k)$ defined as

$$ \frac{1}{k} \sum_{i=1}^{k} Y_{n-i+1,n} \hat{\beta}_{\Lambda} (Y_{n-k+1,n})^{\top} X_{(n-i+1,n)} - \frac{1}{k} \sum_{i=1}^{k} Y_{n-i+1,n} \frac{1}{k} \sum_{i=1}^{k} \hat{\beta}_{\Lambda} (Y_{n-k+1,n})^{\top} X_{(n-i+1,n)}, $$

and where $X_{(n-i+1,n)}$ denotes the concomitant of $Y_{n-i+1,n}$, i.e., the random variable $X_s$ with $s \in \{1, \dots, n\}$ being the unique index such that $Y_s = Y_{n-i+1,n}$. In practice, we let $\hat{k} := \operatorname{argmax}_{5 \le k \le n/5} \bar{r}(k)$ in order to prevent instabilities of estimates built on too few data points.

### 5.3. Results.
The performance of our method is assessed through the comparison of the two graphs $(j/p, \hat{\beta}_{\Lambda}^{(j)}(Y_{n-\hat{k}+1,n}))$ and $(j/p, \beta_j)$ indexed by $1 \le j \le p$. Figures 1–11 display the results on each of the 11 considered dependence frameworks. Each of them includes six or nine panels associated with different set of parameters $(\gamma, \tau, \kappa)$.

In all $9 \times 9 + 2 \times 6 = 93$ situations, the orange curve represents the graph associated with $\beta$, while the blue one is associated with the averaged value of $\hat{\beta}(Y_{n-\hat{k}+1,n})$ over the $N = 500$ Monte-Carlo replications, where $\hat{k}$ is selected following the above described procedure. The light blue area corresponds to the confidence region comprising the top 5–95% values of the Monte Carlo replications.

Figure 1 displays the no-serial dependence case. Adversely, Figure 2 corresponds to a standard ARMA(1, 1)-response $(\phi_{\text{resp}}, \theta_{\text{resp}}) = (0.8, -0.3)$ together with a standard GARCH(1, 1)-noise $(\omega_{\text{noise}}, \alpha_{\text{noise}}, \beta_{\text{noise}}) = (0.05, 0.1, 0.85)$. The results for the IGARCH-like noise are presented in Figure 3. Similarly, Figure 4 is for the pathological ARMA(1, 1)-response with $(\phi_{\text{resp}}, \theta_{\text{resp}}) = (0.99, -0.98)$ with standard GARCH(1, 1)-noise, and Figure 5 when instead the GARCH(1, 1)-noise is close to IGARCH, i.e., $(\omega_{\text{noise}}, \alpha_{\text{noise}}, \beta_{\text{noise}}) = (0.05, 0.05, 0.94)$. Next, the second experimental setup is gathered in Figures 6–9 with mirrored sets of parameters. This is synthesized in Table 1. Finally, Figures 10 and Figure 11 display the case of an ESTAR(1)-response with $\phi_{\text{low}} = 0.2$ and

<table><caption>TABLE 1. Organization of the results on simulated data.</caption><thead><tr><th>Response / Noise</th><th>Indep.</th><th>Stand. ARMA</th><th>Stand. GARCH</th><th>Patho. ARMA</th><th>Patho. GARCH</th></tr></thead><tbody><tr><td>Indep.</td><td>Figure 1</td><td></td><td></td><td></td><td></td></tr><tr><td>Stand. ARMA</td><td></td><td></td><td>Figure 2</td><td></td><td>Figure 3</td></tr><tr><td>Stand. GARCH</td><td></td><td>Figure 6</td><td></td><td>Figure 8</td><td></td></tr><tr><td>Patho. ARMA</td><td></td><td></td><td>Figure 4</td><td></td><td>Figure 5</td></tr><tr><td>Patho. GARCH</td><td></td><td>Figure 7</td><td></td><td>Figure 9</td><td></td></tr><tr><td>ESTAR</td><td></td><td></td><td>Figure 10</td><td></td><td>Figure 11</td></tr></tbody></table>

$\phi_{\text{high}} = 0.95$ with both GARCH(1, 1)-noise models. The first two rows are devoted to a standard noise with $(\omega_{\text{resp}}, \alpha_{\text{resp}}, \beta_{\text{resp}}) = (0.05, 0.1, 0.85)$ for $\gamma \in \{0.1, 0.5\}$ while the last two rows concern the pathological noise with $(\omega_{\text{resp}}, \alpha_{\text{resp}}, \beta_{\text{resp}}) = (0.05, 0.05, 0.94)$.

In all considered situations, the EPLS method provides reliable estimates in terms of bias and variance. It appears to be robust with respect to heavy tails, serial dependence, and missing data, by showing good performance in pathological and nonlinear settings. Unsurprisingly, the estimation procedure becomes more instable as $\tau < 0$ since it implies lighter tails for $\lambda$ and thus more masked observations which reduces the available data sample. On the other hand, a larger tail-index $\gamma$ yields heavier tails for the response and produces extremely large values more frequently, leading to higher variance in tail estimates.
## 6. ILLUSTRATION ON REAL DATA

We exemplify our extreme dimension reduction technique on real-world environmental data from the National Oceanic and Atmospheric Administration (NOAA), the primary U.S. agency for climate observation. Our analysis focuses on the state of Texas, a region frequently affected by extreme weather events, using daily data from NOAA's Global Historical Climatology Network (GHCN) over the five-year period from January 1st, 2020, to December 31st, 2024 (1827 days in total). The data are described in more details in Section 6.1. The competitors as well as the performance assessment method are presented in Section 6.2 and the associated results are provided in Section 6.3.

6.1. **Data selection protocol.** We begin by parsing the complete station inventory provided in ghcnd-stations.txt, which contains metadata (identifier, location, name, etc.) for all GHCN stations worldwide. We retain only stations located within Texas, identified via latitude/longitude filtering or the “TX” state code. From this subset, we further restrict attention to stations with names containing keywords suggestive of institutional or professional meteorological infrastructure (e.g., “airport”, “afb”, “intl”, “weather service”, “observatory”, etc.). This heuristic aims to prioritize stations likely to have more reliable and continuous records, often associated with airports, military bases, or official weather bureaus.

We extract two variables: daily total precipitation (PRCP, in mm) and daily maximum temperature (TMAX, in degrees Celsius). The dataset distinguishes three cases for any variable on a given date: (i) no entry at all, (ii) a NaN value, or (iii) a valid numerical value. We interpret the absence of any entry (i) as indicating that the station was not operational on that day, and exclude such days from the analysis. This simplification removes any masking mechanism related to such absences, under the virtual assumption that no data implies no potential signal. The NaN values, which correspond to known but missing measurements, are imputed to zero. Admittedly, this may lead to an underestimation of the true missing-ness, as certain days with no recorded value could still reflect relevant information. Addressing this issue would require attributing causes to each missing or empty value, potentially through spatial borrowing from nearby stations, which we acknowledge as a meaningful but computationally intensive extension beyond the scope of this study. We retain only those stations where precipitation data is available (i.e., not NaN) for at least 99% of the reduced time range. Any remaining NaN entries are dropped, along with their corresponding dates. This leads to $Y$ having in total $n = 1755$ daily observations. For the maximum temperature, we only consider stations with between 5% and 20% missing values on the same date range as $Y$. From the cleaned dataset, we form all possible triplets of aligned time series $(Y, X^{(1)}, X^{(2)})$, where $Y$ denotes the daily precipitation from one station, and $(X^{(1)}, X^{(2)})$ the maximum temperature from two other stations (which may coincide with $Y$'s location). It yields 27910 possible triplets of stations.

To ensure statistical soundness in our extreme value analysis, we impose a series of filtering criteria on the triplets $(Y, X^{(1)}, X^{(2)})$, concerning both temporal stationarity and tail behavior beyond data completeness. First, we retain only those triplets where the response variable $Y$ (precipitation) satisfies standard weak stationarity conditions. This is verified via the Augmented Dickey-Fuller (ADF) test (testing for the absence of unit roots) and the Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test (testing for stationarity around a constant mean). Both tests must support the stationarity hypothesis at the 5% level. We require joint stationarity for the pair $(X^{(1)}, X^{(2)})$, and
each individual series must again pass both ADF and KPSS tests. Additionally, we enforce cross-correlation stability between $X^{(1)}$ and $X^{(2)}$, as follows: we compute the rolling correlation between the two series (using a fixed window size, e.g., 100 days), and subject the resulting time series of correlations to an ADF test. This tests the temporal constancy of their dependence structure. Only pairs whose rolling correlations are stationary are considered jointly stationary. In addition to the value processes, we retain only the datasets for which the binary masks $M_1$ and $M_2$ flagging missing-ness are jointly stationary. This is assessed via a chi-squared test for homogeneity applied to joint binary patterns over multiple time windows, with Haldane's correction (adding 0.5 to all frequencies to avoid zero counts). The null hypothesis is that the distribution of missing-ness patterns remains stable over time (i.e. the missing-ness process is stationary). We reject it if the p-value is below 5%; thus, only datasets with $p \ge 0.05$ are considered to exhibit stationary missing-ness and are kept for analysis.

Finally, we check whether the marginal distribution of the response variable exhibits heavy-tailed behavior. A standard tool in extreme value theory is the Hill estimator for estimating the tail-index $\gamma$ of a distribution, based on the $k$ largest order statistics of the sample. A Hill plot displays the estimated tail-index as a function of $k$. A visual plateau in such a plot, that is a region where the estimates remain approximately constant across a range of $k$, is typically interpreted as evidence for a regularly varying tail, with the plateau level serving as an estimate of $\gamma$. To enhance reproducibility, we rather implement an automatic procedure to detect heavy-tailness. This method scans the Hill plot for stable plateaus: contiguous regions of order statistics where the estimated tail-index is approximately flat (low slope), low-variance, and exceeds a minimum threshold ($\gamma > 0.2$). A plateau is declared if these criteria are met within a sufficiently long sliding window, ensuring that transient fluctuations are not misclassified as stable signals. If such a region is detected, we conclude that the data are heavy-tailed and retain the mean Hill estimate over the plateau as our tail-index estimate.

Overall, the procedure leads to 5794 admissible triplets $(Y, X^{(1)}, X^{(2)})$ of stations.

## 6.2. Tail covariance evaluation and competitors.
To assess the effectiveness of our method, we propose some custom evaluation pipeline based on the conditional tail covariance and statistics over the admissible triplets of stations.

For a given admissible triplet and a direction $\beta \in \mathbb{R}^2$, we compute the scalar projection $\beta^\top X$ and evaluate its empirical covariance with $Y$, conditionally on $Y$ exceeding a threshold $y$.

A variety of baseline projection competitor strategies is considered: random projections, Extreme PCA (EPCA) computed only on the extreme samples, Extreme Sliced Inverse Regression (ESIR) with a slicing restricted to the top quantile of $Y$, and an extreme variant of Random Forest (ERF) trained to classify extreme events versus non-extreme ones, from which a projection direction is extracted based on feature importance scores. More precisely, fix a threshold $y > 0$ and define the exceedances index set as

$$ \mathcal{I}_y := \{ i \in \{1, \dots, n\} : Y_i > y \}, \quad N_y := |\mathcal{I}_y|. $$

Let us denote by $X_y \in \mathbb{R}^{N_y \times 2}$ the submatrix of $(X_i^{(j)})_{i,j}$ with rows $\{X_i^\top : i \in \mathcal{I}_y\}$ and by $Y_y \in \mathbb{R}^{N_y}$ for the corresponding tail responses. Let $\mathcal{M} := \{\text{EPLS}, \text{EPCA}, \text{LDA}, \text{SIR}, \text{ESIR}, \text{ERF}\}$, then each method $m \in \mathcal{M}$ produces a projection vector score

$$ Z_{m,y} := X_y \beta_m(y) \in \mathbb{R}^{N_y}, $$
where the method-specific direction $\beta_m(y) \in \mathbb{R}^p$ is defined below. For each $m \in \mathcal{M}$, we compute the unbiased empirical covariance between the method's score and the tail responses,

$$ \frac{1}{N_y - 1} \sum_{i \in \mathcal{I}_y} ((Z_{m,y})_i - \bar{Z}_{m,y}) (Y_i - \bar{Y}_y), \qquad (5) $$

with $\bar{Z}_{m,y} = N_y^{-1} \sum_{i \in \mathcal{I}_y} (Z_{m,y})_i$ and $\bar{Y}_y = N_y^{-1} \sum_{i \in \mathcal{I}_y} Y_i$. The four competitors are then computed as follows:

(i) Extreme PCA. The PCA is performed on the tail covariates $\{X_i : i \in \mathcal{I}_y\}$, i.e. $\beta_{EPCA}(y)$ is the top eigenvector of $\hat{\Sigma}_y$ where $\hat{\Sigma}_y = N_y^{-1} \sum_{i \in \mathcal{I}_y} (X_i - \bar{X}_y)(X_i - \bar{X}_y)^\top$. and $\bar{X}_y = N_y^{-1} \sum_{i \in \mathcal{I}_y} X_i$. (ii) Extreme LDA. Let the severity within the tail be $Y_i - y$ for $i \in \mathcal{I}_y$. Partition $\{Y_i - y : i \in \mathcal{I}_y\}$

into five equal-frequency slices and denote the slice index by $s(i) \in \{1, \dots, 5\}$. For each slice $\ell$, define the slice mean as

$$ \mu_{\ell} = \frac{1}{|\{i \in \mathcal{I}_y : s(i) = \ell\}|} \sum_{\substack{i \in \mathcal{I}_y \\ s(i) = \ell}} X_i \in \mathbb{R}^2, $$

and compute the tail covariance $\hat{\Sigma}_y$ as above (on $X_y$). The method forms a contrast vector $\Delta = \mu_K - \frac{1}{K-1} \sum_{k=1}^{K-1} \mu_k$, which captures the shift in $X$ from the less extreme to the most extreme slice. We define $\beta_{LDA}(y) = \hat{\Sigma}_y^+ \Delta$ where $\hat{\Sigma}_y^+$ is the Moore-Penrose pseudo-inverse.

(iii) Extreme Random Forest. A random forest, implemented in scikit-learn, is trained to classify the tail indicator $B_{i,y} := 1_{\{Y_i \ge y\}}$ using all $n$ observations. It yields impurity-based feature importances

$$ \iota = (\iota_1, \iota_2)^\top \in [0, \infty)^2, \quad \iota_1 + \iota_2 = 1, $$

which quantify the relative contribution of each predictor to identifying the extreme region $Y > y$. Using it as a linear proxy of the nonlinear forest, we define the direction as $\beta_{ERF}(y) = \iota/\|\iota\|$.

(iv) Random projections. This naive method consists in sampling 500 directions uniformly on the unit sphere and defining $\beta_R$ as their empirical mean.

(v) Sliced Inverse Regression (SIR). The range of the response $Y$ is divided into $H = 10$ disjoint slices $S_1, \dots, S_H$ (of approximately equal size), based on quantiles of $Y$. Then, for each $h \in \{1, \dots, H\}$, one computes the slice mean $m_h = \frac{1}{N_h} \sum_{i:Y_i \in S_h} X_i$ where $N_h = \#\{i : Y_i \in S_h\}$ is the number of response elements in $S_h$, and also the between-slice covariance matrix

$$ \hat{\Sigma}_B = \sum_{h=1}^{H} \frac{N_h}{n} (m_h - \bar{X})(m_h - \bar{X})^\top, $$

where $\bar{X} = n^{-1} \sum_{i=1}^n X_i$ is the mean of the standardized covariates. The SIR direction $\beta_{SIR}$ is obtained as the leading generalized eigenvector of $\hat{\Sigma}_{XX}^{-1} \hat{\Sigma}_B$, where the sample covariance of $X_i$ is denoted by $\hat{\Sigma}_{XX} = \frac{1}{n} \sum_{i=1}^n (X_i - \bar{X})(X_i - \bar{X})^\top$.

(vi) Extreme Sliced Inverse Regression (Extreme SIR). To obtain a direction $\beta_{ESIR}$, the above procedure is applied only to the exceedances/tail sample $\{(Y_i, X_i) : Y_i > y\}$.
**6.3. Results.** A large-scale analysis is performed over the admissible triplets of stations and across multiple thresholds. The threshold is chosen as $y = q_Y(\alpha)$ where $\alpha \in \{0.90, 0.91, \dots, 0.99\}$. The tail covariances (5) are computed for all methods and ranked from one to six. The evaluation metric then consists in the mean rank (in $\{1, \dots, 6\}$, where 1 is best) achieved by each method over 5794 station triplets. The results are summarized in Figure 12 which represents the evolution of the mean ranks for each method across the different values of the percentile $\alpha$. The mean curves are fairly stable as the percentile increases. Extreme PLS shows consistently superior performance across all thresholds, maintaining the first rank. Next, Extreme PCA and Extreme RF perform similarly, with both methods following the same trend across percentiles. In contrast, Standard SIR, Extreme SIR and Extreme LDA under-perform over most of the threshold range. In the most extreme region, all five competitors tend to converge toward similar performance.

As an illustration, we propose to focus on a single triplet $(Y, X^{(1)}, X^{(2)})$. The threshold considered is the order statistic $Y_{n-k+1,n}$. The tail covariance (5) is plotted in Figure 13 as a function of the threshold index $1 \le k \le n$, for all methods. We also include the random baseline by representing the full range for 500 projections, which displays the range of variation due to randomness. The performance of EPLS is consistent with Figure 12: The EPLS curve closely follows the 'highest random path' in the sense of highest tail covariance among all possibilities. Again, Extreme PCA performs similarly well as Extreme PLS, while Extreme RF does worse and both versions of SIR compare to the average of the random strategy.

The estimates of an extreme conditional quantile (of level $\alpha = 0.99$) associated with the three best performing dimension reduction methods (namely, EPLS, Extreme PCA and Extreme RF) are compared in Figure 14. In all three cases, the explanatory variable is first projected on the estimated dimension reduction axis $\hat{\beta}$, and second, a kernel based estimator of the conditional quantile is computed on the obtained data set $\{(\hat{\beta}^{\top} X_i, Y_i), i = 1, \dots, n\}$. The dimension reduction step has two benefits: it allows us to mitigate the curse of dimensionality in the nonparametric estimator and it moreover yields an easy visualization of the estimated in a two-dimensional space. It appears that EPLS, Extreme PCA and Extreme RF yield similar projected data and thus estimated conditional curves sharing similar shapes.

## 7. APPENDIX: PROOFS

In the sequel we adopt, for any real sequences $A_n$ and $B_n$ such that $B_n \neq 0$ for $n$ large enough, the notation $A_n \lesssim B_n$ if and only if there exists $c \in (0, +\infty)$ independent of $n$ such that $A_n \le cB_n$ for $n$ large enough. Section 7.1 provides three lemmas that will reveal useful to prove the main results in Section 7.2.

**7.1. Preliminary results.** The first lemma is an adaptation of (Bingham et al., 1987, Proposition 1.5.10) to the missing data setting. It provides an asymptotic equivalent of conditional expectations above a large threshold.

**Lemma 1.** Let $\rho \in RV_\mu$ with $\mu \in \mathbb{R}$ and let $Y$ be a random variable satisfying (A$_1$).

(i) If $\gamma\mu < 1$, then, as $y \to \infty$,

$$ \mathbb{E}[\rho(Y) | Y \ge y] = \frac{m_{\rho(Y)}(y)}{\hat{F}(y)} \sim \frac{1}{1 - \gamma\mu} \rho(y). $$
(ii) If (A$_4$) holds, with $\gamma(\mu + \tau) < 1$, then, as $y \to \infty$, for all $(j \neq l) \in \{1, \dots, p\}^2$,
$$ \mathbb{E}[\Lambda^{(j)} \rho(Y) 1_{\{Y \geq y\}}] \sim \frac{c_j}{1 - \gamma(\mu + \tau)} \rho(y) \lambda(y) \bar{F}(y). $$

(iii) If (A$_4$) holds, with $\gamma(\mu + 2\tau) < 1$, then, as $y \to \infty$, for all $(j \neq l) \in \{1, \dots, p\}^2$,
$$ \mathbb{E}[\Lambda^{(j)} \Lambda^{(l)} \rho(Y) 1_{\{Y \geq y\}}] \sim \frac{c_j c_l}{1 - \gamma(\mu + 2\tau)} \rho(y) \lambda^2(y) \bar{F}(y). $$

(iv) If (A$_3$) and (A$_4$) hold with $\gamma(\frac{q}{q-1}\mu + \tau) < 1$, then, as $y \to \infty$, for all $(j \neq l) \in \{1, \dots, p\}^2$,
$$ \mathbb{E}[\Lambda^{(j)} \varepsilon^{(j)} \rho(Y) 1_{\{Y \geq y\}}] = O\left(\rho(y)(\lambda \bar{F})^{1-1/q}(y)\right). $$

(v) If (A$_3$) and (A$_4$) hold with $\gamma(\frac{q}{q-2-t}\mu + \tau) < 1$ for $t \ge 0$, then, as $y \to \infty$, for all $(j \neq l) \in \{1, \dots, p\}^2$,
$$ \mathbb{E}[\Lambda^{(j)}(\varepsilon^{(j)})^{2+t} \rho(Y) 1_{\{Y \geq y\}}] = O\left(\rho(y)(\lambda \bar{F})^{1-(2+t)/q}(y)\right). $$

**Proof.** (i) is established in (Bousebata et al., 2023, Lemma 2).

(ii) Let us consider

$$ \begin{aligned} \mathbb{E}[\Lambda^{(j)} \rho(Y) 1_{\{Y \geq y\}}] &= \mathbb{E}\left\{\mathbb{E}[\Lambda^{(j)} \rho(Y) 1_{\{Y \geq y\}} | Y]\right\} = \mathbb{E}\left\{\rho(Y) 1_{\{Y \geq y\}} \mathbb{E}[\Lambda^{(j)} | Y]\right\} \\ &= \mathbb{E}\left\{\rho(Y) \lambda_j(Y) 1_{\{Y \geq y\}}\right\} = c_j \mathbb{E}\left\{\rho(Y) \lambda(Y) | Y \geq y\right\} \bar{F}(y). \end{aligned} $$

Since $\rho(\cdot)\lambda(\cdot) \in RV_{\mu+\tau}$ and $\gamma(\mu + \tau) < 1$, part (i) of the Lemma proves the result.

(iii) The proof is similar after noting that $\mathbb{E}[\Lambda^{(j)}\Lambda^{(l)} | Y] = \lambda_j(Y)\lambda_l(Y)$ for $j \neq l$.

(iv) By Hölder's, for any $q' \ge 1$ such that $1/q + 1/q' = 1$, one has

$$ \mathbb{E}[\Lambda^{(j)} \varepsilon^{(j)} \rho(Y) 1_{\{Y \geq y\}}] \le \mathbb{E}(|\varepsilon^{(j)}|^q)^{1/q} \mathbb{E}(\rho(Y)^{q'} \Lambda^{(j)} 1_{\{Y \geq y\}})^{1/q'}. $$

In view of (ii) together with $\gamma(\frac{q}{q-1}\mu + \tau) < 1$ and (A$_3$), the result follows.

(v) Again, Hölder inequality implies that for any $q' \ge 1$ such that $(2+t)/q + 1/q' = 1$,

$$ \mathbb{E}[\Lambda^{(j)}(\varepsilon^{(j)})^2 \rho(Y) 1_{\{Y \geq y_n\}}] \le \mathbb{E}(|\varepsilon^{(j)}|^q)^{2/q} \mathbb{E}(|\rho(Y)|^{q'} \Lambda^{(j)} 1_{\{Y \geq y_n\}})^{1/q'}. $$

In view of (ii) together with $\gamma(\frac{q}{q-2-t}\mu + \tau) < 1$ and (A$_3$), the result follows. $\square$

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

The next lemma makes use of Lemma 1 to derive asymptotic equivalents of conditional covariances above a large threshold in the missing data framework.

**Lemma 2.** Suppose $\rho_1 \in RV_{\mu_1}$, $\rho_2 \in RV_{\mu_2}$ with $\mu_1, \mu_2 \in \mathbb{R}$ and let $Y$ be a random variable satisfying (A$_1$). Assume that (A$_4$) holds.

(i) Assume $\gamma(\mu_1 + \mu_2 + 2\tau) < 1$, $\gamma(\mu_s + \tau) < 1$ for $s \in \{1, 2\}$ holds. Then, as $y \to \infty$, for all $(j \neq l) \in \{1, \dots, p\}^2$,
$$ \operatorname{cov}[\Lambda^{(j)} \rho_1(Y) 1_{\{Y \geq y\}}, \Lambda^{(l)} \rho_2(Y) 1_{\{Y \geq y\}}] \sim \frac{c_j c_l}{1 - \gamma(\mu_1 + \mu_2 + 2\tau)} \rho_1(y) \rho_2(y) \lambda^2(y) \bar{F}(y). $$

(ii) Assume $\gamma(\mu_1 + \mu_2 + \tau) < 1$, $\gamma(\mu_s + \tau) < 1$ for $s \in \{1, 2\}$ holds. Then, as $y \to \infty$, for all $(j \neq l) \in \{1, \dots, p\}^2$,
$$ \operatorname{cov}[\Lambda^{(j)} \rho_1(Y) 1_{\{Y \geq y\}}, \Lambda^{(l)} \rho_2(Y) 1_{\{Y \geq y\}}] \sim \frac{c_j}{1 - \gamma(\mu_1 + \mu_2 + \tau)} \rho_1(y) \rho_2(y) \lambda(y) \bar{F}(y). $$

(iii) If, moreover, $\gamma(\frac{q}{q-2}(\mu_1 + \mu_2) + 2\tau) < 1$ and $\gamma(\frac{q}{q-1}\mu_s + \tau) < 1$ for $s \in \{1, 2\}$ and (A$_3$) hold, then, as $y \to \infty$, for all $(j \neq l) \in \{1, \dots, p\}^2$,
$$ \operatorname{cov}[\Lambda^{(j)} \varepsilon^{(j)} \rho_1(Y) 1_{\{Y \geq y\}}, \Lambda^{(l)} \varepsilon^{(l)} \rho_2(Y) 1_{\{Y \geq y\}}] = O\left(\rho_1(y)\rho_2(y)\lambda^{2(1-2/q)}(y)\bar{F}^{1-2/q}(y)\right). $$
(iv) If $\gamma(\frac{q}{q-1}(\mu_1 + \mu_2) + 2\tau) < 1$ and $(A_3)$ hold, then, as $y \to \infty$, for all $(j \neq l) \in \{1, \dots, p\}^2$, 
$$ \operatorname{cov} \left[ \Lambda^{(j)} \varepsilon^{(j)} \rho_1(Y) \mathbf{1}_{\{Y \ge y\}}, \Lambda^{(l)} \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right] = O \left( \rho_1(y) \rho_2(y) \lambda^{2(1-1/q)}(y) \bar{F}^{1-1/q}(y) \right). $$

(v) If $\gamma(\frac{q}{q-1}(\mu_1 + \mu_2) + \mu_2 + \tau) < 1$, $\gamma(\mu_2 + \tau) < 1$ and $(A_3)$ hold, then, as $y \to \infty$, for all $j \in \{1, \dots, p\}$,

$$ \operatorname{cov} \left[ \Lambda^{(j)} \varepsilon^{(j)} \rho_1(Y) \mathbf{1}_{\{Y \ge y\}}, \Lambda^{(j)} \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right] = O \left( \rho_1(y) \rho_2(y) (\lambda \bar{F})^{1-2/q}(y) \right). $$

**Proof.** (i) Taking account of

$$ \operatorname{cov} \left[ \Lambda^{(j)} \rho_1(Y) \mathbf{1}_{\{Y \ge y\}}, \Lambda^{(l)} \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right] = \mathbb{E} \left\{ \Lambda^{(j)} \Lambda^{(l)} \rho_1(Y) \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right\} \\ - \mathbb{E} \left\{ \Lambda^{(j)} \rho_1(Y) \mathbf{1}_{\{Y \ge y\}} \right\} \mathbb{E} \left\{ \Lambda^{(l)} \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right\} $$

and applying Lemma 1(ii,iii) show that, for $j \neq l$:

$$ \operatorname{cov} \left[ \Lambda^{(j)} \rho_1(Y) \mathbf{1}_{\{Y \ge y\}}, \Lambda^{(l)} \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right] \sim \mathbb{E} \left\{ \Lambda^{(j)} \Lambda^{(l)} \rho_1(Y) \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right\} $$

as $y \to \infty$ and the result follows.

(ii) Second,

$$ \operatorname{cov} \left[ \Lambda^{(j)} \rho_1(Y) \mathbf{1}_{\{Y \ge y\}}, \Lambda^{(j)} \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right] = \mathbb{E} \left\{ \Lambda^{(j)} \rho_1(Y) \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right\} \\ - \mathbb{E} \left\{ \Lambda^{(j)} \rho_1(Y) \mathbf{1}_{\{Y \ge y\}} \right\} \mathbb{E} \left\{ \Lambda^{(j)} \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right\} $$

and Lemma 1(ii) implies

$$ \operatorname{cov} \left[ \Lambda^{(j)} \rho_1(Y) \mathbf{1}_{\{Y \ge y\}}, \Lambda^{(j)} \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right] \sim \mathbb{E} \left\{ \Lambda^{(j)} \rho_1(Y) \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right\} $$

as $y \to \infty$ so that the result is proved.

(iii) The covariance

$$ C_{\rho_2, \varepsilon, l}^{\rho_1, \varepsilon, j} := \operatorname{cov} \left[ \Lambda^{(j)} \varepsilon^{(j)} \rho_1(Y) \mathbf{1}_{\{Y \ge y\}}, \Lambda^{(l)} \varepsilon^{(l)} \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right] $$

may be expanded as:

$$ \mathbb{E} \left\{ \Lambda^{(j)} \Lambda^{(l)} \varepsilon^{(j)} \varepsilon^{(l)} \rho_1(Y) \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right\} - \mathbb{E} \left\{ \Lambda^{(j)} \varepsilon^{(j)} \rho_1(Y) \mathbf{1}_{\{Y \ge y\}} \right\} \mathbb{E} \left\{ \Lambda^{(l)} \varepsilon^{(l)} \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right\}. $$

The first term is bounded using Hölder inequality, applied for all $q' \ge 1$ such that $2/q + 1/q' = 1$:

$$ \mathbb{E}(|\varepsilon^{(j)}|^q)^{1/q} \mathbb{E}(|\varepsilon^{(l)}|^q)^{1/q} \mathbb{E}((\rho_1(Y)\rho_2(Y))^{q'} \Lambda^{(j)} \Lambda^{(l)} \mathbf{1}_{\{Y \ge y\}})^{1/q'}. $$

Similarly, for all $q'' \ge 1$ such that $1/q + 1/q'' = 1$, the product of the second and third terms is bounded above by

$$ \mathbb{E}(|\varepsilon^{(j)}|^q)^{1/q} \mathbb{E}(|\varepsilon^{(l)}|^q)^{1/q} \mathbb{E}(\rho_1(Y)^{q''} \Lambda^{(j)} \mathbf{1}_{\{Y \ge y\}})^{1/q''} \mathbb{E}(\rho_2(Y)^{q''} \Lambda^{(l)} \mathbf{1}_{\{Y \ge y\}})^{1/q''}. $$

In view of Lemma 1(ii,iii) together with $\gamma(\frac{q}{q-2}(\mu_1+\mu_2)+2\tau) < 1$ and $\gamma(\frac{q}{q-1}\mu_s+\tau) < 1$ for $s \in \{1, 2\}$ and $(A_3)$, it follows that for any $q' > 2, q'' > 1$ such that $2/q + 1/q' = 1$ and $1/q + 1/q'' = 1$,

$$ \begin{aligned} C_{\rho_2, \varepsilon, l}^{\rho_1, \varepsilon, j} &= O((\rho_1 \cdot \rho_2 \cdot \lambda^{2/q'} \cdot \bar{F}^{1/q'}(y)) \lor (\rho_1 \cdot \rho_2 \cdot \lambda^{2/q''} \cdot \bar{F}(y)^{2/q''})) \\ &= O((\rho_1 \cdot \rho_2 \cdot \lambda^{2/q'} \cdot \bar{F}^{1/q'}(y)). \end{aligned} $$

The result follows by using $2/q + 1/q' = 1$.

(iv) It readily follows from the same Hölder type argument with $1/p + 1/q = 1$ and an application of Lemma 1(iii).
(v) The covariance

$$C_{\rho_2, j}^{\rho_1, \varepsilon, j} := \operatorname{cov} \left[ \Lambda^{(j)} \varepsilon^{(j)} \rho_1(Y) \mathbf{1}_{\{Y \ge y\}}, \Lambda^{(j)} \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right]$$

may be expanded as:

$$\mathbb{E} \left\{ \Lambda^{(j)} \varepsilon^{(j)} \rho_1(Y) \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right\} - \mathbb{E} \left\{ \Lambda^{(j)} \varepsilon^{(j)} \rho_1(Y) \mathbf{1}_{\{Y \ge y\}} \right\} \mathbb{E} \left\{ \Lambda^{(j)} \rho_2(Y) \mathbf{1}_{\{Y \ge y\}} \right\}.$$

The first term is bounded using Hölder inequality, applied for all $q' \ge 1$ such that $1/q + 1/q' = 1$:

$$|\mathbb{E}\{\Lambda^{(j)} \varepsilon^{(j)} \rho_1(Y) \rho_2(Y) \mathbf{1}_{\{Y \ge y\}}\}| \le \mathbb{E}(|\varepsilon^{(j)}|^q)^{1/q} \mathbb{E}((\rho_1(Y)\rho_2(Y))^{q'} \Lambda^{(j)} \mathbf{1}_{\{Y \ge y\}})^{1/q'}.$$

Similarly, the second term is bounded above:

$$|\mathbb{E}\{\Lambda^{(j)} \varepsilon^{(j)} \rho_1(Y) \mathbf{1}_{\{Y \ge y\}}\}| \le \mathbb{E}(|\varepsilon^{(j)}|^q)^{1/q} \mathbb{E}(\rho_1(Y)^{q'} \Lambda^{(j)} \mathbf{1}_{\{Y \ge y\}})^{1/q'}.$$

In view of Lemma 1(ii) together with $\gamma(\frac{q}{q-1}(\mu_1 + \mu_2) \lor \mu_2 + \tau) < 1$, $\gamma(\mu_2 + \tau) < 1$ and (A3), it follows that for any $q' > 1$ such that $1/q + 1/q' = 1$,

$$C_{\rho_2, j}^{\rho_1, \varepsilon, j} = O \left( \rho_1(y) \rho_2(y) (\lambda \bar{F})^{1-2/q}(y) \right).$$

The result is thus proved.

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

Finally, Lemma 1 also allows us to derive asymptotic equivalents of conditional tail moments in the missing data framework.

**Lemma 3.** Let $Y$ be a random variable satisfying (A1). Let $\iota \in \{0, 1\}$.

(i) Suppose (A4) holds. Then, as $y \to \infty$,

$$m_{\lambda^{(j)}}(y) \sim \frac{c_j}{1 - \gamma\tau} \lambda(y) \bar{F}(y),$$

for all $j \in \{1, \dots, p\}$.

(ii) Assume (M1), (A2) and (A3) hold with $\gamma(\kappa + \iota) < 1$. Then, as $y \to \infty$,

$$m_{Y^\iota X^{(j)}}(y) \sim \frac{\beta_j}{1 - \gamma(\kappa + \iota)} y^\iota g(y) \bar{F}(y),$$

for all $j \in \{1, \dots, d\}$.

(iii) Assume (M1), (A2), (A3) and (A4) hold with $\gamma(\kappa + \tau + \iota) < 1$, and $1 < \gamma(q\kappa + \tau)$. Then, as $y \to \infty$,

$$m_{Y^\iota \Lambda^{(j)} X^{(j)}}(y) \sim \frac{\beta_j c_j}{1 - \gamma(\kappa + \tau + \iota)} y^\iota g(y) \lambda(y) \bar{F}(y),$$

for all $j \in \{1, \dots, d\}$.

**Proof.** (i) One may apply Lemma 1(i) since $\gamma\tau < 1$ in view of $\tau < 0$, to get for $j \in \{1, \dots, p\}$ as $y \to \infty$:

$$m_{\lambda^{(j)}}(y) \sim \frac{c_j}{1 - \gamma\tau} \lambda(y) \bar{F}(y).$$

(ii) (Bousebata et al., 2023, Lemma 7(i)) yields for $j \in \{1, \dots, d\}$, as long as $\gamma(\kappa + \iota) < 1$,

$$m_{Y^\iota X^{(j)}}(y) \sim \frac{\beta_j}{1 - \gamma(\kappa + s)} y^\iota g(y) \bar{F}(y),$$

as $y \to \infty$.

(iii) According to model (M1), the following decomposition holds for $j \in \{1, \dots, d\}$:

$$m_{Y\Lambda^{(j)}X^{(j)}}(y) = \beta_j m_{Y\Lambda^{(j)}g(Y)}(y) + m_{Y\Lambda^{(j)}\varepsilon^{(j)}}(y).$$
The asymptotics of Lemma 1(ii,iv) together with assumptions (A3) and (A4), readily yield, as $y \to \infty$, under $\gamma(\kappa + \tau + \iota) < 1$ and $\gamma(\frac{q}{q-1} + \tau) < 1$,

$$m_{Y^\iota \Lambda^{(j)} g(Y)}(y) \sim \frac{\beta_j c_j}{1 - \gamma(\kappa + \tau + \iota)} y^\iota g(y) \lambda(y) \bar{F}(y),$$

$$m_{Y^\iota \Lambda^{(j)} \varepsilon^{(j)}(y)} = O\left(y^\iota \lambda(y) (y^{1-1/q} \bar{F}(y))^{1-1/q}\right) = o(m_{Y^\iota \Lambda^{(j)} g(Y)}(y)),$$

where the last negligibility is induced by the condition $1 < \gamma(q\kappa + \tau)$. Thus, it follows that, for any $j \in \{1, \dots, d\}$, as $y \to \infty$,

$$m_{Y^\iota \Lambda^{(j)} X^{(j)}(y)} = \frac{\beta_j c_j}{1 - \gamma(\kappa + \tau + \iota)} y^\iota g(y) \lambda(y) \bar{F}(y) (1 + o(1))$$

and the result is proved.

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

## 7.2. Proof of main results.

**Proof of Proposition 2.** Consider $\delta > 0$ given in condition (3) of strong mixing, $\theta_0 := 1 + \frac{\delta}{2+\delta}$ and $\theta \ge \theta_0$. In virtue of the Cramér–Wold device and to establish the joint asymptotic normality of the $(3d+1)$- random vector $\Xi_n$, we shall prove that any non-zero linear combination of its components is asymptotically Gaussian. Consider then $(a_{1,1}, \dots, a_{1,d}, a_{2,1}, \dots, a_{2,d}, a_{3,1}, \dots, a_{3,d}, a_{4}, a_{5})^\top \in \mathbb{R}^{3d+2}$ deterministic and let us investigate the asymptotic distribution of $\tilde{\chi}_n$ defined as follows:

$$\chi_n = \sqrt{n(\lambda \bar{F})^\theta(y_n)} \left\{ \sum_{k=1}^{3} \sum_{j=1}^{d} a_{k,j} \chi_{j,n}^{(k)} + a_4 \chi_n^{(4)} + a_5 \chi_n^{(5)} \right\} := \sum_{i=1}^{n} \chi_{i,n},$$

where $\chi_{i,n} := \sqrt{\frac{(\lambda \bar{F})^\theta(y_n)}{n}} (\sum_{k=1}^3 \chi_{i,n}^{(k)} + \chi_{i,n}^{(4)} + \chi_{n}^{(5)})$ with, for any $1 \le j \le d$ and for $\iota \in \{0,1\}$,

$$\chi_{i,n}^{(1+\iota)} := \sum_{j=1}^{d} a_{1,j} \left( \frac{Y_i^\iota \Lambda_i^{(j)} X_i^{(j)}}{m_{Y^\iota \Lambda^{(j)} X^{(j)}(y_n)}} \mathbb{1}\{Y_i \ge y_n\} - 1 \right), \quad \chi_{i,n}^{(2)} := \sum_{j=1}^{d} a_{2,j} \left( \frac{\Lambda_i^{(j)}}{m_{\lambda^{(j)}(y_n)}} \mathbb{1}\{Y_i \ge y_n\} - 1 \right),$$

$$\chi_{i,n}^{(4+\iota)} := a_{4+\iota} \left( Y_i^\iota \frac{\mathbb{1}\{Y_i \ge y_n\}}{m_{Y^\iota}(y_n)} - 1 \right).$$

By measurable mapping, the $\alpha$-mixing hypothesis implies that $(\chi_{i,n})_{1 \le i \le n}$ is a $\alpha$-mixing tri- angular, row-strictly stationary array of real random variables. Define for $1 \le i < j \le n$, $\mathcal{F}_{n,i,j} = \sigma(\chi_{n,i}, \chi_{n,i+1}, \dots, \chi_{n,j})$ and the associated strong mixing coefficient by

$$\forall n \ge 2, \quad \forall \ell < n, \quad \alpha_n(\ell) := \sup\{|\mathbb{P}(A \cap B) - \mathbb{P}(A)\mathbb{P}(B)|, A \in \mathcal{F}_{n,1,k}, B \in \mathcal{F}_{n,k+\ell,n}, 1 \le k \le n - \ell\}.$$

Then, it readily follows that $\alpha_n \le \alpha$ and one is in the right setting to apply (Davison et al., 2023, Lemma C.7) to the sequence $(\chi_{i,n})_{1 \le i \le n}$. The four requirements of this result are

$$\frac{n}{r_n} \operatorname{var} \left( \sum_{i=1}^{\ell_n} \chi_{i,n} \right) = \frac{n}{r_n} \ell_n \operatorname{var}(\chi_{1,n}) + 2 \frac{n}{r_n} \sum_{1 \le i_1 < i_2 \le \ell_n} \operatorname{cov}(\chi_{i_1,n}, \chi_{i_2,n}) \ll 1, \quad (6)$$

$$\frac{n}{r_n} \operatorname{var} \left( \sum_{i=1}^{r_n} \chi_{i,n} \right) = n \operatorname{var}(\chi_{1,n}) + 2 \frac{n}{r_n} \sum_{1 \le i_1 < i_2 \le r_n} \operatorname{cov}(\chi_{i_1,n}, \chi_{i_2,n}) \to \sigma^2 < \infty, \quad (7)$$

$$\operatorname{var} \left( \sum_{i=1}^{n-r_n \lfloor n/r_n \rfloor} \chi_{i,n} \right) = (n-r_n \lfloor n/r_n \rfloor) \operatorname{var}(\chi_{1,n}) + 2 \sum_{1 \le i_1 < i_2 \le n-r_n \lfloor n/r_n \rfloor} \operatorname{cov}(\chi_{i_1,n}, \chi_{i_2,n}) \\ \le r_n \operatorname{var}(\chi_{1,n}) + 2 \sum_{1 \le i_1 < i_2 \le n-r_n \lfloor n/r_n \rfloor} \operatorname{cov}(\chi_{i_1,n}, \chi_{i_2,n}) \ll 1, \quad (8)$$
and

$$ \frac{n}{r_n} \mathbb{E} \left( \left| \sum_{i=1}^{r_n} \chi_{i,n} \right|^{2+\delta} \right) \ll 1. \qquad (9) $$

The first step is to study the variance $\text{var}(\chi_{1,n})$ involved in (6)–(8). An expansion yields:

$$ \text{var}(\chi_{1,n}) = \frac{(\lambda \bar{F})^{\theta}(y_n)}{n} \sum_{k=1}^{5} \text{var}(\chi_{1,n}^{(k)}) + 2 \frac{(\lambda \bar{F})^{\theta}(y_n)}{n} \sum_{k_1 < k_2} \text{cov}(\chi_{1,n}^{(k_1)}, \chi_{1,n}^{(k_2)}), \quad (10) $$

where each term of variance and covariance will be addressed separately. We begin by collecting the variances, starting with the first case $k = 1$. Using the standard variance/covariance properties, as well as the identity $\text{cov}(X - \mathbb{E}X, Y - \mathbb{E}Y) = \text{cov}(X, Y)$, one gets

$$ \begin{aligned} \text{var}(\chi_{i,n}^{(1)}) &= \text{var} \left( \sum_{j=1}^{d} a_{1,j} \left( \frac{\Lambda_i^{(j)} X_i^{(j)}}{m_{\Lambda^{(j)}X^{(j)}}(y_n)} \mathbb{1}_{\{Y_i \ge y_n\}} - 1 \right) \right) \\ &= \sum_{j=1}^{d} a_{1,j}^2 \frac{\text{var}(\Lambda_i^{(j)} X_i^{(j)} \mathbb{1}_{\{Y_i \ge y_n\}})}{m_{\Lambda^{(j)}X^{(j)}}^2(y_n)} \\ &\quad + 2 \sum_{1 \le j_1 < j_2 \le d} a_{1,j_1} a_{1,j_2} \frac{\text{cov}(\Lambda_i^{(j_1)} X_i^{(j_1)} \mathbb{1}_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} X_i^{(j_2)} \mathbb{1}_{\{Y_i \ge y_n\}})}{m_{\Lambda^{(j_1)}X^{(j_1)}}(y_n) m_{\Lambda^{(j_2)}X^{(j_2)}}(y_n)}. \end{aligned} $$

Under (M1), the first term $\text{var}(\Lambda_i^{(j)} X_i^{(j)} \mathbb{1}_{\{Y_i \ge y_n\}})/m_{\Lambda^{(j)}X^{(j)}}^2(y_n)$ may be factorized as

$$ \beta_j^2 \frac{m_{g^2(Y)\Lambda^{(j)}(y_n)}}{m_{\Lambda^{(j)}X^{(j)}}^2(y_n)} \left\{ 1 + \frac{m_{\Lambda^{(j)}\varepsilon^{(j)}(y_n)}}{\beta_j^2 m_{g^2(Y)\Lambda^{(j)}(y_n)}} + 2 \frac{m_{g(Y)\Lambda^{(j)}\varepsilon^{(j)}(y_n)}}{\beta_j m_{g^2(Y)\Lambda^{(j)}(y_n)}} - \frac{m_{\Lambda^{(j)}X^{(j)}}^2(y_n)}{\beta_j^2 m_{g^2(Y)\Lambda^{(j)}(y_n)}} \right\}. $$

Using Lemma 3(iii) under $\gamma(\kappa + \tau) < 1$ and $1 < \gamma(q\kappa + \tau)$, Lemma 1(ii) under $\gamma(2\kappa + \tau) < 1$ and Lemma 1(iv) under $\gamma(\frac{q}{q-1}\kappa + \tau) < 1$, it follows:

$$ \sum_{j=1}^{d} a_{1,j}^{2} \frac{\text{var}(\Lambda_{i}^{(j)} X_{i}^{(j)} \mathbb{1}_{\{Y_{i} \ge y_{n}\}})}{m_{\Lambda^{(j)}X^{(j)}}^{2}(y_{n})} = \frac{(1+o(1))}{\lambda(y_{n})\bar{F}(y_{n})} \frac{(1-\gamma(\kappa+\tau))^{2}}{1-\gamma(2\kappa+\tau)} \sum_{j=1}^{d} \frac{a_{1,j}^{2}}{c_{j}}. $$

Besides, under (M1), the covariance $\text{cov}(\Lambda_i^{(j_1)} X_i^{(j_1)} \mathbb{1}_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} X_i^{(j_2)} \mathbb{1}_{\{Y_i \ge y_n\}})$ is

$$ \begin{aligned} & \beta_{j_1} \beta_{j_2} \text{cov}(g(Y_i)\Lambda_i^{(j_1)} \mathbb{1}_{\{Y_i \ge y_n\}}, g(Y_i)\Lambda_i^{(j_2)} \mathbb{1}_{\{Y_i \ge y_n\}}) \\ & + \beta_{j_1} \text{cov}(g(Y_i)\Lambda_i^{(j_1)} \mathbb{1}_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} \varepsilon_i^{(j_2)} \mathbb{1}_{\{Y_i \ge y_n\}}) \\ & + \beta_{j_2} \text{cov}(\Lambda_i^{(j_1)} \varepsilon_i^{(j_1)} \mathbb{1}_{\{Y_i \ge y_n\}}, g(Y_i)\Lambda_i^{(j_2)} \mathbb{1}_{\{Y_i \ge y_n\}}) \\ & + \text{cov}(\Lambda_i^{(j_1)} \varepsilon_i^{(j_1)} \mathbb{1}_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} \varepsilon_i^{(j_2)} \mathbb{1}_{\{Y_i \ge y_n\}}). \end{aligned} $$

The first term is $\frac{c_{j_1} c_{j_2} \beta_{j_1} \beta_{j_2}}{1-2\gamma(\kappa+\tau)} g^2(y_n) \lambda^2(y_n) \bar{F}(y_n)(1+o(1))$ by Lemma 2(i) under $2\gamma(\kappa+\tau) < 1$. The second and third terms are $O(g(y_n)\lambda^{2(1-1/q)}(y_n)\bar{F}^{1-1/q}(y_n))$ by Lemma 2(iv) under $\gamma(\frac{q}{q-1}\kappa+2\tau) < 1$. The fourth term is $O(\lambda^{2(1-2/q)}(y_n)\bar{F}^{1-2/q}(y_n))$ using Lemma 2(iii) and is negligible compared to the second one. The latter is in turn negligible compared to the first one under $1 < \gamma(2\kappa+2q+\tau)$. Gathering everything, it follows that

$$ \text{cov}(\Lambda_i^{(j_1)} X_i^{(j_1)} \mathbb{1}_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} X_i^{(j_2)} \mathbb{1}_{\{Y_i \ge y_n\}}) = \frac{\beta_{j_1} \beta_{j_2} c_{j_1} c_{j_2}}{1 - 2\gamma(\kappa + \tau)} g^2(y_n) (\lambda \bar{F})(y_n) (1 + o(1)), $$
as $n \to \infty$. Then, by Lemma 3(iii) under $\gamma(\kappa + \tau) < 1$ and $1 < \gamma(q\kappa + \tau)$ for the denominator, it yields

$$ \mathrm{var}(\chi_{i,n}^{(1)}) = \frac{(1+o(1))}{\lambda(y_n)\bar{F}(y_n)} \left\{ \frac{(1-\gamma(\kappa+\tau))^2}{1-\gamma(2\kappa+\tau)} \sum_{j=1}^{d} \frac{a_{1,j}^2}{c_j} + \frac{(1-\gamma(\kappa+\tau))^2}{1-2\gamma(\kappa+\tau)} \sum_{1\le j_1 < j_2 \le d} 2a_{1,j_1}a_{1,j_2} \right\}. \quad (11) $$

We proceed in the same way for $\chi_{i,n}^{(2)}$. One has to begin:

$$ \begin{aligned} \mathrm{var}(\chi_{i,n}^{(2)}) &= \mathrm{var} \left( \sum_{j=1}^{d} a_{2,j} \left( \frac{Y_i \Lambda_i^{(j)} X_i^{(j)}}{m_{Y\Lambda^{(j)}X^{(j)}}(y_n)} \mathbb{1}_{\{Y_i \ge y_n\}} - 1 \right) \right) \\ &= \sum_{j=1}^{d} a_{2,j}^2 \frac{\mathrm{var}(Y_i \Lambda_i^{(j)} X_i^{(j)} \mathbb{1}_{\{Y_i \ge y_n\}})}{m_{Y\Lambda^{(j)}X^{(j)}}^2(y_n)} \\ &\quad + 2 \sum_{1 \le j_1 < j_2 \le d} a_{2,j_1} a_{2,j_2} \frac{\mathrm{cov}(Y_i \Lambda_i^{(j_1)} X_i^{(j_1)} \mathbb{1}_{\{Y_i \ge y_n\}}, Y_i \Lambda_i^{(j_2)} X_i^{(j_2)} \mathbb{1}_{\{Y_i \ge y_n\}})}{m_{Y\Lambda^{(j_1)}X^{(j_1)}}(y_n) m_{Y\Lambda^{(j_2)}X^{(j_2)}}(y_n)}. \end{aligned} $$

Under (M₁), the term $\mathrm{var}(Y_i\Lambda_i^{(j)}X_i^{(j)}\mathbb{1}_{\{Y_i \ge y_n\}})/m_{Y\Lambda^{(j)}X^{(j)}}^2(y_n)$ may be factorized as

$$ \beta_j^2 \frac{m_{Y^2 g^2(Y)\Lambda^{(j)}}(y_n)}{m_{Y\Lambda^{(j)}X^{(j)}}^2(y_n)} \left\{ 1 + \frac{m_{Y\Lambda^{(j)}\varepsilon^{(j)}}(y_n)}{\beta_j^2 m_{Y^2 g^2(Y)\Lambda^{(j)}}(y_n)} + 2 \frac{m_{Yg(Y)\Lambda^{(j)}\varepsilon^{(j)}}(y_n)}{\beta_j m_{Y^2 g^2(Y)\Lambda^{(j)}}(y_n)} - \frac{m_{Y\Lambda^{(j)}X^{(j)}}^2(y_n)}{\beta_j^2 m_{Y^2 g^2(Y)\Lambda^{(j)}}(y_n)} \right\}. $$

Using Lemma 3(iii) under $\gamma(\kappa+\tau+1) < 1$ and $1 < \gamma(q(\kappa+1)+\tau)$, Lemma 1(ii) under $\gamma(2\kappa+\tau+2) < 1$ and Lemma 1(iv) under $\gamma(\frac{q}{q-1}\kappa + \tau) < 1$, it follows that

$$ \sum_{j=1}^{d} a_{2,j}^2 \frac{\mathrm{var}(Y_i \Lambda_i^{(j)} X_i^{(j)} \mathbb{1}_{\{Y_i \ge y_n\}})}{m_{Y\Lambda^{(j)}X^{(j)}}^2(y_n)} = \frac{(1+o(1))}{\lambda(y_n)\bar{F}(y_n)} \frac{(1-\gamma(\kappa+\tau+1))^2}{1-\gamma(2\kappa+\tau+2)} \sum_{j=1}^{d} \frac{a_{2,j}^2}{c_j}. $$

Besides, under (M₁), the covariance $\mathrm{cov}(Y_i\Lambda_i^{(j_1)}X_i^{(j_1)}\mathbb{1}_{\{Y_i \ge y_n\}}, Y_i\Lambda_i^{(j_2)}X_i^{(j_2)}\mathbb{1}_{\{Y_i \ge y_n\}})$ is

$$ \begin{aligned} & \beta_{j_1} \beta_{j_2} \mathrm{cov}(Y_i g(Y_i) \Lambda_i^{(j_1)} \mathbb{1}_{\{Y_i \ge y_n\}}, Y_i g(Y_i) \Lambda_i^{(j_2)} \mathbb{1}_{\{Y_i \ge y_n\}}) \\ & + \beta_{j_1} \mathrm{cov}(Y_i g(Y_i) \Lambda_i^{(j_1)} \mathbb{1}_{\{Y_i \ge y_n\}}, Y_i \Lambda_i^{(j_2)} \varepsilon_i^{(j_2)} \mathbb{1}_{\{Y_i \ge y_n\}}) \\ & + \beta_{j_2} \mathrm{cov}(Y_i \Lambda_i^{(j_1)} \varepsilon_i^{(j_1)} \mathbb{1}_{\{Y_i \ge y_n\}}, Y_i g(Y_i) \Lambda_i^{(j_2)} \mathbb{1}_{\{Y_i \ge y_n\}}) \\ & + \mathrm{cov}(Y_i \Lambda_i^{(j_1)} \varepsilon_i^{(j_1)} \mathbb{1}_{\{Y_i \ge y_n\}}, Y_i \Lambda_i^{(j_2)} \varepsilon_i^{(j_2)} \mathbb{1}_{\{Y_i \ge y_n\}}). \end{aligned} $$

Factorizing by the first and highest term, using Lemma 2(i) under $2\gamma(\kappa + \tau + 1) < 1$, Lemma 2(iv) under $\gamma(\frac{q}{q-1}\kappa + 2\tau) < 1$ and Lemma 2(iii), one obtains

$$ \mathrm{cov}(Y_i \Lambda_i^{(j_1)} X_i^{(j_1)} \mathbb{1}_{\{Y_i \ge y_n\}}, Y_i \Lambda_i^{(j_2)} X_i^{(j_2)} \mathbb{1}_{\{Y_i \ge y_n\}}) = \frac{\beta_{j_1} \beta_{j_2} c_{j_1} c_{j_2}}{1 - 2\gamma(\kappa + \tau + 1)} y_n^2 g^2(y_n) (\lambda \bar{F})(y_n) (1 + o(1)), $$

as $n \to \infty$. Then, by Lemma 3(iii) under $\gamma(\kappa + \tau + 1) < 1$ and $1 < \gamma(q\kappa + \tau)$ for the denominator, it yields

$$ \mathrm{var}(\chi_{i,n}^{(2)}) = \frac{(1+o(1))}{\lambda(y_n)\bar{F}(y_n)} \left\{ \frac{(1-\gamma(\kappa+\tau+1))^2}{1-\gamma(2\kappa+\tau+2)} \sum_{j=1}^{d} \frac{a_{2,j}^2}{c_j} + \frac{(1-\gamma(\kappa+\tau+1))^2}{1-2\gamma(\kappa+\tau+1)} \sum_{1\le j_1 < j_2 \le d} 2a_{2,j_1}a_{2,j_2} \right\}. \quad (12) $$
Similarly, using Lemma 3(i) for the denominator and Lemma 2(i) for the covariance entails

$$
\begin{align}
\mathrm{var}(\chi_{i,n}^{(3)}) &= \sum_{j=1}^{d} a_{3,j}^2 \frac{\mathrm{var}(\Lambda_i^{(j)} 1_{\{Y_i \ge y_n\}})}{m_{\lambda^{(j)}(y_n)}^2} + 2 \sum_{1 \le j_1 < j_2 \le d} a_{3,j_1} a_{3,j_2} \frac{\mathrm{cov}(\Lambda_i^{(j_1)} 1_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} 1_{\{Y_i \ge y_n\}})}{m_{\lambda^{(j_1)}(y_n)} m_{\lambda^{(j_2)}(y_n)}} \nonumber \\
&= \frac{(1 + o(1))}{\lambda(y_n) \bar{F}(y_n)} (1 - \gamma \tau) \sum_{j=1}^{d} \frac{a_{3,j}^2}{c_j} + O\left(\frac{1}{\bar{F}(y_n)}\right) = \frac{(1 + o(1))}{\lambda(y_n) \bar{F}(y_n)} (1 - \gamma \tau) \sum_{j=1}^{d} \frac{a_{3,j}^2}{c_j}, \tag{13}
\end{align}
$$

since $\tau < 0$. Finally, the last sub-variance is easily treated with the variance of a Bernoulli distribution:

$$
\mathrm{var}(\chi_{i,n}^{(4)}) = \mathrm{var}\left(a_4 \left(\frac{1}{\bar{F}(y_n)} 1_{\{Y_i \ge y_n\}} - 1\right)\right) = \frac{a_4^2 (1 - \bar{F}(y_n))}{\bar{F}(y_n)} \sim \frac{a_4^2}{\bar{F}(y_n)}, \quad (14)
$$

while, using Lemma 1(i) under $\gamma < 1/2$,

$$
\mathrm{var}(\chi_{i,n}^{(5)}) = a_5^2 \frac{\mathrm{var}(Y_i 1_{\{Y_i \ge y_n\}})}{m_Y^2(y_n)} = a_5^2 \left( \frac{m_{Y^2}(y_n)}{m_Y^2(y_n)} - 1 \right) \sim \frac{(1-\gamma)^2}{1-2\gamma} \frac{1}{\bar{F}(y_n)}. \quad (15)
$$

Now, we focus on the sub-covariances in $\mathrm{var}(\chi_{1,n})$. The usual properties of the covariance operator entail

$$
\begin{align*}
\mathrm{cov}(\chi_{i,n}^{(1)}, \chi_{i,n}^{(3)}) &= \mathrm{cov} \left( \sum_{j=1}^{d} a_{1,j} \left( \frac{\Lambda_i^{(j)} X_i^{(j)}}{m_{\Lambda^{(j)}X^{(j)}(y_n)} 1_{\{Y_i \ge y_n\}} - 1 \right), \sum_{j=1}^{d} a_{3,j} \left( \frac{\Lambda_i^{(j)}}{m_{\Lambda^{(j)}(y_n)} 1_{\{Y_i \ge y_n\}} - 1 \right) \right) \\
&= \sum_{j_1,j_2=1}^{d} a_{1,j_1} a_{3,j_2} \frac{\mathrm{cov}(\Lambda_i^{(j_1)} X_i^{(j_1)} 1_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} 1_{\{Y_i \ge y_n\}})}{m_{\Lambda^{(j_1)}X^{(j_1)}(y_n)} m_{\Lambda^{(j_2)}(y_n)}} \\
&= \sum_{j_1 \neq j_2} a_{1,j_1} a_{3,j_2} \frac{\mathrm{cov}(\Lambda_i^{(j_1)} X_i^{(j_1)} 1_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} 1_{\{Y_i \ge y_n\}})}{m_{\Lambda^{(j_1)}X^{(j_1)}(y_n)} m_{\Lambda^{(j_2)}(y_n)}} \\
&\quad + \sum_{j=1}^{d} a_{1,j} a_{3,j} \frac{\mathrm{cov}(\Lambda_i^{(j)} X_i^{(j)} 1_{\{Y_i \ge y_n\}}, \Lambda_i^{(j)} 1_{\{Y_i \ge y_n\}})}{m_{\Lambda^{(j)}X^{(j)}(y_n)} m_{\Lambda^{(j)}(y_n)}}.
\end{align*}
$$

Let us start with the double sum in $\mathrm{cov}(\chi_{i,n}^{(1)}, \chi_{i,n}^{(3)})$ by writing, under (M$_1$), that

$$
\begin{align*}
\mathrm{cov}(\Lambda_i^{(j_1)} X_i^{(j_1)} 1_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} 1_{\{Y_i \ge y_n\}}) &= \beta_{j_1} \mathrm{cov}(g(Y_i) \Lambda_i^{(j_1)} 1_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} 1_{\{Y_i \ge y_n\}}) \\
&\quad + \mathrm{cov}(\Lambda_i^{(j_1)} \varepsilon_i^{(j_1)} 1_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} 1_{\{Y_i \ge y_n\}}).
\end{align*}
$$

By Lemma 2(i) under $\gamma(\kappa + 2\tau) < 1$, $\gamma\kappa < 1$, the first term is $O((g\lambda^2\bar{F})(y_n))$; while Lemma 2(iv) under $\gamma(\frac{q}{q-1} + 2\tau) < 1$ shows that the second term is $O((\lambda^{2(1-1/q)}\bar{F}^{1-1/q})(y_n))$. Moreover, $(\lambda^{2(1-1/q)}\bar{F}^{1-1/q})(y_n) \ll (g\lambda^2\bar{F})(y_n)$ under $1 < \gamma(q\kappa + 2\tau)$, so that using Lemma 3(i), Lemma 3(iii) and Lemma 2(iv) under $\gamma(\frac{q}{q-1} + 2\tau) < 1$ and $1 < \gamma(q\kappa + \tau)$, one has, since $d < +\infty$,

$$
\sum_{j_1 \neq j_2} a_{1,j_1} a_{3,j_2} \frac{\mathrm{cov}(\Lambda_i^{(j_1)} X_i^{(j_1)} 1_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} 1_{\{Y_i \ge y_n\}})}{m_{\Lambda^{(j_1)}X^{(j_1)}(y_n)} m_{\Lambda^{(j_2)}(y_n)}} = O\left(\frac{1}{\bar{F}(y_n)}\right).
$$

The second double sum in $\mathrm{cov}(\chi_{i,n}^{(1)}, \chi_{i,n}^{(3)})$ is treated in the similar fashion, by first decomposing under (M$_1$) and then using this time Lemma 2(v). All-in-one, it yields:

$$
\mathrm{cov}(\chi_{i,n}^{(1)}, \chi_{i,n}^{(3)}) = \frac{1}{\lambda(y_n) \bar{F}(y_n)} \sum_{j=1}^{d} \frac{a_{1,j} a_{3,j}}{c_j}. \qquad (16)
$$
Concerning cov($\chi_{i,n}^{(2)}$, $\chi_{i,n}^{(3)}$), we begin by writing:

$$
\begin{align*}
\operatorname{cov}(\chi_{i,n}^{(2)}, \chi_{i,n}^{(3)}) &= \operatorname{cov} \left( \sum_{j=1}^{d} a_{2,j} \left( \frac{Y_i \Lambda_i^{(j)} X_i^{(j)}}{m_{Y\Lambda^{(j)}X^{(j)}}(y_n)} \mathbf{1}_{\{Y_i \ge y_n\}} - 1 \right), \sum_{j=1}^{d} a_{3,j} \left( \frac{\Lambda_i^{(j)}}{m_{\Lambda^{(j)}}(y_n)} \mathbf{1}_{\{Y_i \ge y_n\}} - 1 \right) \right) \\
&= \sum_{j_1, j_2=1}^{d} a_{2,j_1} a_{3,j_2} \frac{\operatorname{cov}(Y_i \Lambda_i^{(j_1)} X_i^{(j_1)} \mathbf{1}_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} \mathbf{1}_{\{Y_i \ge y_n\}})}{m_{Y\Lambda^{(j_1)}X^{(j_1)}}(y_n) m_{\Lambda^{(j_2)}}(y_n)} \\
&= \sum_{j_1 \neq j_2} a_{2,j_1} a_{3,j_2} \frac{\operatorname{cov}(Y_i \Lambda_i^{(j_1)} X_i^{(j_1)} \mathbf{1}_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} \mathbf{1}_{\{Y_i \ge y_n\}})}{m_{Y\Lambda^{(j_1)}X^{(j_1)}}(y_n) m_{\Lambda^{(j_2)}}(y_n)} \\
&\quad + \sum_{j=1}^{d} a_{2,j} a_{3,j} \frac{\operatorname{cov}(Y_i \Lambda_i^{(j)} X_i^{(j)} \mathbf{1}_{\{Y_i \ge y_n\}}, \Lambda_i^{(j)} \mathbf{1}_{\{Y_i \ge y_n\}})}{m_{Y\Lambda^{(j)}X^{(j)}}(y_n) m_{\Lambda^{(j)}}(y_n)}.
\end{align*}
$$

Let us start with the double sum in cov($\chi_{i,n}^{(2)}$, $\chi_{i,n}^{(3)}$) by remarking that, under (M₁),

$$
\begin{align*}
\operatorname{cov}(Y_i \Lambda_i^{(j_1)} X_i^{(j_1)} \mathbf{1}_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} \mathbf{1}_{\{Y_i \ge y_n\}}) ={}& \beta_{j_1} \operatorname{cov}(Y_i g(Y_i) \Lambda_i^{(j_1)} \mathbf{1}_{\{Y_i \ge y_n\}}, \Lambda^{(j_2)} \mathbf{1}_{\{Y_i \ge y_n\}}) \\
& + \operatorname{cov}(Y_i \Lambda_i^{(j_1)} \varepsilon_i^{(j_1)} \mathbf{1}_{\{Y_i \ge y_n\}}, \Lambda^{(j_2)} \mathbf{1}_{\{Y_i \ge y_n\}}).
\end{align*}
$$

By Lemma 2(i) under $\gamma(\kappa + 2\tau + 1) < 1$, $\gamma(\kappa + 1) < 1$, the first term is $O(y_n(g\lambda^2\bar{F})(y_n))$; while Lemma 2(iv) under $\gamma(\frac{q}{q-1} + 2\tau) < 1$ shows that the second term is $O(y_n(\lambda^{2(1-1/q)}\bar{F}^{1-1/q})(y_n))$. Moreover, $(\lambda^{2(1-1/q)}\bar{F}^{1-1/q})(y_n) \ll (g\lambda^2\bar{F})(y_n)$ under $1 < \gamma(q\kappa + 2\tau)$, so that using Lemma 3(i), Lemma 3(iii) and Lemma 2(iv) under $\gamma(\frac{q}{q-1} + 2\tau) < 1$ and $1 < \gamma(q\kappa + \tau)$, the leading term is the first one and by comparison, since $d < +\infty$,

$$
\sum_{j_1 \neq j_2} a_{2,j_1} a_{3,j_2} \frac{\operatorname{cov}(Y_i \Lambda_i^{(j_1)} X_i^{(j_1)} \mathbf{1}_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} \mathbf{1}_{\{Y_i \ge y_n\}})}{m_{\Lambda^{(j_1)}X^{(j_1)}}(y_n) m_{\Lambda^{(j_2)}}(y_n)} = O\left(\frac{1}{\bar{F}(y_n)}\right).
$$

The second double sum in cov($\chi_{i,n}^{(2)}$, $\chi_{i,n}^{(3)}$) is treated in the similar fashion, by first decomposing under (M₁) and then using this time Lemma 2(v) under $\gamma(\frac{q}{q-1} + \tau) < 1$. All-in-one, it yields:

$$
\operatorname{cov}(\chi_{i,n}^{(2)}, \chi_{i,n}^{(3)}) = \frac{1}{\lambda(y_n)\bar{F}(y_n)} \sum_{j=1}^{d} \frac{a_{2,j} a_{3,j}}{c_j}. \qquad (17)
$$

Regarding cov($\chi_{i,n}^{(2)}$, $\chi_{i,n}^{(1)}$), we begin by writing:

$$
\begin{align*}
\operatorname{cov}(\chi_{i,n}^{(2)}, \chi_{i,n}^{(1)}) &= \operatorname{cov} \left( \sum_{j=1}^{d} a_{2,j} \left( \frac{Y_i \Lambda_i^{(j)} X_i^{(j)}}{m_{Y\Lambda^{(j)}X^{(j)}}(y_n)} \mathbf{1}_{\{Y_i \ge y_n\}} - 1 \right), \sum_{j=1}^{d} a_{1,j} \left( \frac{\Lambda_i^{(j)} X_i^{(j)}}{m_{\Lambda^{(j)}X^{(j)}}(y_n)} \mathbf{1}_{\{Y_i \ge y_n\}} - 1 \right) \right) \\
&= \sum_{j_1, j_2=1}^{d} a_{2,j_1} a_{1,j_2} \frac{\operatorname{cov}(Y_i \Lambda_i^{(j_1)} X_i^{(j_1)} \mathbf{1}_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} X_i^{(j_2)} \mathbf{1}_{\{Y_i \ge y_n\}})}{m_{Y\Lambda^{(j_1)}X^{(j_1)}}(y_n) m_{\Lambda^{(j_2)}X^{(j_2)}}(y_n)} \\
&= \sum_{j_1 \neq j_2} a_{2,j_1} a_{1,j_2} \frac{\operatorname{cov}(Y_i \Lambda_i^{(j_1)} X_i^{(j_1)} \mathbf{1}_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} X_i^{(j_2)} \mathbf{1}_{\{Y_i \ge y_n\}})}{m_{Y\Lambda^{(j_1)}X^{(j_1)}}(y_n) m_{\Lambda^{(j_2)}X^{(j_2)}}(y_n)} \\
&\quad + \sum_{j=1}^{d} a_{1,j} a_{2,j} \frac{\operatorname{cov}(Y_i \Lambda_i^{(j)} X_i^{(j)} \mathbf{1}_{\{Y_i \ge y_n\}}, \Lambda_i^{(j)} X_i^{(j)} \mathbf{1}_{\{Y_i \ge y_n\}})}{m_{Y\Lambda^{(j)}X^{(j)}}(y_n) m_{\Lambda^{(j)}X^{(j)}}(y_n)}.
\end{align*}
$$
We start with the double sum in $cov(\chi_{i,n}^{(2)}, \chi_{i,n}^{(1)})$ by writing, under (M₁), that

$$
\begin{align*}
\operatorname{cov}(Y_i \Lambda_i^{(j_1)} X_i^{(j_1)} 1_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} X_i^{(j_2)} 1_{\{Y_i \ge y_n\}}) &= \beta_{j_1} \beta_{j_2} \operatorname{cov}(Y_i g(Y_i) \Lambda_i^{(j_1)} 1_{\{Y_i \ge y_n\}}, g(Y_i) \Lambda^{(j_2)} 1_{\{Y_i \ge y_n\}}) \\
&\quad + \beta_{j_1} \operatorname{cov}(Y_i g(Y_i) \Lambda_i^{(j_1)} 1_{\{Y_i \ge y_n\}}, \Lambda^{(j_2)} \varepsilon_i^{(j_2)} 1_{\{Y_i \ge y_n\}}) \\
&\quad + \beta_{j_2} \operatorname{cov}(\Lambda_i^{(j_1)} \varepsilon_i^{(j_1)} 1_{\{Y_i \ge y_n\}}, Y_i g(Y_i) \Lambda^{(j_2)} 1_{\{Y_i \ge y_n\}}) \\
&\quad + \operatorname{cov}(Y_i \Lambda_i^{(j_1)} \varepsilon_i^{(j_1)} 1_{\{Y_i \ge y_n\}}, \Lambda^{(j_2)} \varepsilon_i^{(j_2)} 1_{\{Y_i \ge y_n\}}).
\end{align*}
$$

By Lemma 2(i) under $\gamma(2\kappa+2\tau+1) < 1$, the first term is $O(y_n(g^2\lambda^2\bar{F})(y_n))$; while Lemma 2(iv) under $\gamma(\frac{q}{q-1}(\kappa+1)+2\tau) < 1$ shows that the second and third term are $O(y_n g(y_n)\lambda^{2(1-1/q)}(y_n)\bar{F}^{1-1/q}(y_n))$.
The fourth term is $O(y_n\lambda^{2(1-2/q)}(y_n)\bar{F}^{1-2/q}(y_n))$ thanks to Lemma 2(iii) under $\gamma(\frac{q}{q-2} + 2\tau) < 1$ and $\gamma(\frac{q}{q-1} + \tau) < 1$. Doing again a comparison, it holds that $(\lambda^{2(1-1/q)}\bar{F}^{1-1/q})(y_n) \ll (g\lambda^2\bar{F})(y_n)$ under $1 < \gamma(q\kappa+2\tau)$, so that using Lemma 3(i), Lemma 3(iii) and Lemma 2(iv) under $\gamma(\frac{q}{q-1}+2\tau) < 1$ and $1 < \gamma(q\kappa + \tau)$, the leading term is the first one and, since $d < +\infty$,

$$
\sum_{j_1 \neq j_2} a_{2,j_1} a_{1,j_2} \frac{\operatorname{cov}(Y_i \Lambda_i^{(j_1)} X_i^{(j_1)} 1_{\{Y_i \ge y_n\}}, \Lambda_i^{(j_2)} X_i^{(j_2)} 1_{\{Y_i \ge y_n\}})}{m_{Y\Lambda^{(j_1)}X^{(j_1)}}(y_n) m_{\Lambda^{(j_2)}X^{(j_2)}}(y_n)} = O\left(\frac{1}{\bar{F}(y_n)}\right).
$$

The case $cov(\chi_{i,n}^{(k_1)}, \chi_{i,n}^{(k_2)})$ with $k_1 \in \{1,2\}$, $k_2 \in \{4,5\}$ is treated as follows. Let $\iota, s \in \{0,1\}$ and using the definition of the tail-moments, one has, since $d < +\infty$,

$$
\begin{align*}
\operatorname{cov}(\chi_{i,n}^{(1+\iota)}, \chi_{i,n}^{(4+s)}) &= \operatorname{cov} \left( \sum_{j=1}^{d} a_{1+\iota,j} \left( \frac{Y_i^{\iota} \Lambda_i^{(j)} X_i^{(j)}}{m_{Y^{\iota} \Lambda^{(j)} X^{(j)}}(y_n) 1_{\{Y_i \ge y_n\}} - 1} \right), a_{4+s} \left( \frac{Y_i^{s} 1_{\{Y_i \ge y_n\}}}{m_{Y^s}(y_n)} - 1 \right) \right) \\
&= a_{4+s} \sum_{j=1}^{d} a_{1+\iota,j} \frac{\operatorname{cov}(Y_i^{\iota} \Lambda_i^{(j)} X_i^{(j)} 1_{\{Y_i \ge y_n\}}, Y_i^s 1_{\{Y_i \ge y_n\}})}{m_{Y^{\iota} \Lambda^{(j)} X^{(j)}}(y_n) m_{Y^s}(y_n)} \\
&= a_{4+s} \sum_{j=1}^{d} a_{1+\iota,j} \frac{\mathbb{E}(Y_i^{\iota+s} \Lambda_i^{(j)} X_i^{(j)} 1_{\{Y \ge y\}}) - \mathbb{E}(Y_i^{\iota} \Lambda_i^{(j)} X_i^{(j)} 1_{\{Y \ge y\}}) \mathbb{E}(Y_i^s 1_{\{Y \ge y\}})}{m_{Y^{\iota} \Lambda^{(j)} X^{(j)}}(y_n) m_{Y^s}(y_n)} \\
&= a_{4+s} \sum_{j=1}^{d} a_{1+\iota,j} \frac{\mathbb{E}(Y_i^{\iota+s} \Lambda_i^{(j)} X_i^{(j)} 1_{\{Y \ge y\}})}{m_{Y^{\iota} \Lambda^{(j)} X^{(j)}}(y_n) m_{Y^s}(y_n)} (1 + o(1)) = O\left(\frac{1}{\bar{F}(y_n)}\right).
\end{align*}
$$

Besides, for the case $k_1 = 3$ and $k_2 \in \{4,5\}$, since $d < +\infty$,

$$
\begin{align*}
\operatorname{cov}(\chi_{i,n}^{(3)}, \chi_{i,n}^{(4+s)}) &= \operatorname{cov} \left( \sum_{j=1}^{d} a_{3,j} \left( \frac{\Lambda_i^{(j)}}{m_{\lambda^{(j)}}(y_n)} 1_{\{Y_i \ge y_n\}} - 1 \right), a_{4+s} \left( \frac{Y_i^s 1_{\{Y_i \ge y_n\}}}{m_{Y^s}(y_n)} - 1 \right) \right) \\
&= a_{4+s} \sum_{j=1}^{d} a_{3,j} \frac{\operatorname{cov}(\Lambda_i^{(j)} 1_{\{Y_i \ge y_n\}}, Y_i^s 1_{\{Y_i \ge y_n\}})}{m_{\lambda^{(j)}}(y_n) m_{Y^s}(y_n)} \\
&= a_{4+s} \sum_{j=1}^{d} a_{3,j} \frac{m_{Y^s \Lambda^{(j)}}(y_n) - m_{\lambda^{(j)}}(y_n) m_{Y^s}(y_n)}{m_{\lambda^{(j)}}(y_n) m_{Y^s}(y_n)} = O\left(\frac{1}{\bar{F}(y_n)}\right).
\end{align*}
$$

Finally,

$$
\begin{align*}
\operatorname{cov}(\chi_{i,n}^{(5)}, \chi_{i,n}^{(4)}) &= \operatorname{cov} \left( a_5 \left( \frac{Y_i}{m_Y(y_n)} 1_{\{Y_i \ge y_n\}} - 1 \right), a_4 \left( \frac{1}{\bar{F}(y_n)} - 1 \right) \right) \\
&= a_4 a_5 \frac{\operatorname{cov}(Y_i 1_{\{Y_i \ge y_n\}}, 1_{\{Y_i \ge y_n\}})}{m_Y(y_n) \bar{F}(y_n)} = a_4 a_5 \frac{\mathbb{E}(Y_i 1_{\{Y \ge y\}})(1 - \bar{F}(y_n))}{m_Y(y_n) \bar{F}(y_n)} = O\left(\frac{1}{\bar{F}(y_n)}\right).
\end{align*}
$$
Combining (11)–(17) into (10), it follows that, for some $\varphi_1$ depending on $d$, $\gamma$, $\kappa$, $\tau$ and $a_{k,j}$ for $1 \le k \le 5$ and $1 \le j \le d$, since $d < +\infty$,

$$ \mathrm{var}(\tilde{\chi}_{1,n}) = \frac{(\lambda \bar{F})^{\theta-1}(y_n)}{n} \varphi_1(1 + o(1)). \quad (18) $$

The second step is the control of the covariances in (6)–(8) and, to this aim, let us consider the expansion, for $i_1 \neq i_2$,

$$ \mathrm{cov}(\chi_{i_1,n}, \chi_{i_2,n}) = \frac{(\lambda \bar{F})^{\theta}(y_n)}{n} \sum_{k=1}^{4} \mathrm{cov}(\chi_{i_1,n}^{(k)}, \chi_{i_2,n}^{(k)}) + 2 \frac{(\lambda \bar{F})^{\theta}(y_n)}{n} \sum_{k_1<k_2} \mathrm{cov}(\chi_{i_1,n}^{(k_1)}, \chi_{i_2,n}^{(k_2)}). $$

Our goal now is to prove that, for any $k_1, k_2 \in \{1, 2, 3, 4, 5\}$ and $i_1, i_2 \le n$, it holds that

$$ \mathrm{cov}(\chi_{i_1,n}^{(k_1)}, \chi_{i_2,n}^{(k_2)}) \lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) (\lambda \bar{F})^{\frac{2}{2+\delta}-2} (y_n). $$

Invoking (Doukhan, 1994, Theorem 3) with herein $p = q = 2 + \delta$, one may write, for $\iota \in \{0, 1\}$,

$$ \begin{aligned} \mathrm{cov}(\chi_{i_1,n}^{(1+\iota)}, \chi_{i_2,n}^{(1+\iota)}) &= \sum_{j_1,j_2=1}^{d} a_{1+\iota,j_1} a_{1+\iota,j_2} \frac{\mathrm{cov}(Y_{i_1}^{\iota} \Lambda_{i_1}^{(j_1)} X_{i_1}^{(j_1)} 1_{\{Y_{i_1} \ge y_n\}}, Y_{i_2}^{\iota} \Lambda_{i_2}^{(j_2)} X_{i_2}^{(j_2)} 1_{\{Y_{i_2} \ge y_n\}})}{m_{Y^{\iota} \Lambda^{(j_1)} X^{(j_1)}(y_n)} m_{Y^{\iota} \Lambda^{(j_2)} X^{(j_2)}(y_n)}} \\ &\le 8 \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \sum_{j_1,j_2=1}^{d} a_{1+\iota,j_1} a_{1+\iota,j_2} \frac{\prod_{s=1}^{2} \frac{\mathrm{E}(Y^{\iota} \Lambda^{(j_s)} (X^{(j_s)})^{2+\delta} 1_{\{Y \ge y_n\}})}{m_{Y^{\iota} \Lambda^{(j_1)} X^{(j_1)}(y_n)} m_{Y^{\iota} \Lambda^{(j_2)} X^{(j_2)}(y_n)}}^{\frac{1}{2+\delta}}. \end{aligned} $$

Then, under model (M1), and by a Jensen-type argument, it follows that

$$ \begin{aligned} \mathrm{cov}(\chi_{i_1,n}^{(1+\iota)}, \chi_{i_2,n}^{(1+\iota)}) &\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \sum_{j_1,j_2=1}^{d} a_{1+\iota,j_1} a_{1+\iota,j_2} \frac{\max_{j \le d} \frac{\mathrm{E}(Y^{\iota} g^{2+\delta} (Y) \Lambda^{(j)} 1_{\{Y \ge y_n\}})^{\frac{2}{2+\delta}}}{m_{Y^{\iota} \Lambda^{(j_1)} X^{(j_1)}(y_n)} m_{Y^{\iota} \Lambda^{(j_2)} X^{(j_2)}(y_n)}} \\ &\quad + \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \sum_{j_1,j_2=1}^{d} a_{1+\iota,j_1} a_{1+\iota,j_2} \frac{\max_{j \le d} \frac{\mathrm{E}(Y^{\iota} \Lambda^{(j)} (\varepsilon^{(j)})^{2+\delta} 1_{\{Y \ge y_n\}})^{\frac{2}{2+\delta}}}{m_{Y^{\iota} \Lambda^{(j_1)} X^{(j_1)}(y_n)} m_{Y^{\iota} \Lambda^{(j_2)} X^{(j_2)}(y_n)}}. \end{aligned} $$

Next, using Lemma 3(iii) under $1 < \gamma(q\kappa + \tau)$ for the denominator and Lemma 1(ii) under $\gamma((2+\delta)\kappa + \tau + 1) < 1$ together with Lemma 1(v) under $\gamma(\frac{q}{q-2-\delta} + \tau) < 1$ for respectively the first and second numerator terms, one may write under (A3) for the negligibility in the last line, for $n$ large enough, for any $i_2 \neq i_2$, since $d < +\infty$ and $1 < \gamma(q\kappa + \tau)$,

$$ \begin{aligned} \mathrm{cov}(\chi_{i_1,n}^{(1+\iota)}, \chi_{i_2,n}^{(1+\iota)}) &\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \left\{ \frac{y_n^{\iota \frac{2}{2+\delta}} g^2(y_n) (\lambda \bar{F})^{\frac{2}{2+\delta}} (y_n)}{y_n^2 g^2(y_n) (\lambda \bar{F})^2(y_n)} + \frac{y_n^{\iota \frac{2}{2+\delta}} (\lambda \bar{F})^{\frac{2}{2+\delta} - \frac{2}{q}} (y_n)}{y_n^2 g^2(y_n) (\lambda \bar{F})^2(y_n)} \right\} \\ &= \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) y_n^{\iota \frac{2}{2+\delta} - 2} (\lambda \bar{F})^{\frac{2}{2+\delta} - 2} (y_n) \left\{ 1 + g^{-2}(y_n) (\lambda \bar{F})^{-\frac{2}{q}} (y_n) \right\} \\ &\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) (\lambda \bar{F})^{\frac{2}{2+\delta} - 2} (y_n). \end{aligned} \quad (19) $$
Lemma 3(i) and (Doukhan, 1994, Theorem 3) give that, for $n$ large enough, for any $i_2 \neq i_2$, since $d < +\infty$,

$$
\begin{align*}
\operatorname{cov}(\chi_{i_1,n}^{(3)}, \chi_{i_2,n}^{(3)}) &= \operatorname{cov} \left( \sum_{j=1}^{d} a_{3,j} \left( \frac{\Lambda_{i_1}^{(j)}}{m_{\lambda^{(j)}}(y_n)} \mathbf{1}_{\{Y_{i_1} \ge y_n\}} - 1 \right), \sum_{j=1}^{d} a_{3,j} \left( \frac{\Lambda_{i_2}^{(j)}}{m_{\lambda^{(j)}}(y_n)} \mathbf{1}_{\{Y_{i_2} \ge y_n\}} - 1 \right) \right) \\
&= \sum_{j_1, j_2=1}^{d} a_{3,j_1} a_{3,j_2} \frac{\operatorname{cov}(\Lambda_{i_1}^{(j_1)} \mathbf{1}_{\{Y_{i_1} \ge y_n\}}, \Lambda_{i_2}^{(j_2)} \mathbf{1}_{\{Y_{i_2} \ge y_n\}})}{m_{\lambda^{(j_1)}}(y_n) m_{\lambda^{(j_2)}}(y_n)} \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \sum_{j_1, j_2=1}^{d} a_{3,j_1} a_{3,j_2} \frac{\max_{j \le d} \{m_{\lambda^{(j)}}(y_n)\}^{\frac{2}{2+\delta}}}{m_{\lambda^{(j_1)}}(y_n) m_{\lambda^{(j_2)}}(y_n)} \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) (\lambda \bar{F})^{\frac{2}{2+\delta} - 2}(y_n). \tag{20}
\end{align*}
$$

Once more, (Doukhan, 1994, Theorem 3) herein with $p = q = 2 + \delta$ provides, for $\iota \in \{0, 1\}$,

$$
\operatorname{cov}(\chi_{i_1,n}^{(4+\iota)}, \chi_{i_2,n}^{(4+\iota)}) = a_{4+\iota}^2 \frac{\operatorname{cov}(Y_{i_1}^{\iota} \mathbf{1}_{\{Y_{i_1} \ge y_n\}}, Y_{i_2}^{\iota} \mathbf{1}_{\{Y_{i_2} \ge y_n\}})}{m_{Y^{\iota}}^2(y_n)} \le 8a_{4+\iota}^2 \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \bar{F}^{\frac{2}{2+\delta}-2}(y_n). \quad (21)
$$

Concerning the cross-covariances, one may write in virtue of (Doukhan, 1994, Theorem 3),

$$
\begin{align*}
\operatorname{cov}(\chi_{i_1,n}^{(1)}, \chi_{i_2,n}^{(3)}) &= \operatorname{cov} \left( \sum_{j=1}^{d} a_{1,j} \left( \frac{\Lambda_{i_1}^{(j)} X_{i_1}^{(j)}}{m_{\Lambda^{(j)} X^{(j)}(y_n)} \mathbf{1}_{\{Y_{i_1} \ge y_n\}} - 1 \right), \sum_{j=1}^{d} a_{3,j} \left( \frac{\Lambda_{i_2}^{(j)}}{m_{\Lambda^{(j)}(y_n)} \mathbf{1}_{\{Y_{i_2} \ge y_n\}} - 1 \right) \right) \\
&= \sum_{j_1, j_2=1}^{d} \alpha_{1,j_1} \alpha_{2,j_2} \frac{\operatorname{cov}(\Lambda_{i_1}^{(j_1)} X_{i_1}^{(j_1)} \mathbf{1}_{\{Y_{i_1} \ge y_n\}}, \Lambda_{i_2}^{(j_2)} \mathbf{1}_{\{Y_{i_2} \ge y_n\}})}{m_{\Lambda^{(j_1)} X^{(j_1)}(y_n)} m_{\Lambda^{(j_2)}(y_n)} \\
&\le 8 \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \sum_{j_1, j_2=1}^{d} a_{1,j_1} a_{1,j_2} \frac{\mathbb{E}(\Lambda^{(j_1)}(X^{(j_1)})^{2+\delta} \mathbf{1}_{\{Y \ge y_n\}})^{\frac{1}{2+\delta}} \mathbb{E}(\Lambda^{(j_2)} \mathbf{1}_{\{Y \ge y_n\}})^{\frac{1}{2+\delta}}}{m_{\Lambda^{(j_1)} X^{(j_1)}(y_n)} m_{\Lambda^{(j_2)}(y_n)}}.
\end{align*}
$$

Now, in view of (M1), one can decompose the first expectation to get

$$
\begin{align*}
\operatorname{cov}(\chi_{i_1,n}^{(1)}, \chi_{i_2,n}^{(3)}) &\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \sum_{j_1, j_2=1}^{d} a_{1,j_1} a_{3,j_2} \frac{\max_{j \le d} \{\mathbb{E}(g^{2+\delta}(Y) \Lambda^{(j)} \mathbf{1}_{\{Y \ge y_n\}}) \cdot m_{\lambda^{(j)}}(y_n)\}^{\frac{1}{2+\delta}}}{m_{\Lambda^{(j_1)} X^{(j_1)}(y_n)} m_{\lambda^{(j_2)}}(y_n)} \\
&\quad + \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \sum_{j_1, j_2=1}^{d} a_{1,j_1} a_{3,j_2} \frac{\max_{j \le d} \{\mathbb{E}(\Lambda^{(j)}(\varepsilon^{(j)})^{2+\delta} \mathbf{1}_{\{Y \ge y_n\}}) \cdot m_{\lambda^{(j)}}(y_n)\}^{\frac{1}{2+\delta}}}{m_{\Lambda^{(j_1)} X^{(j_1)}(y_n)} m_{\lambda^{(j_2)}}(y_n)}.
\end{align*}
$$

By Lemma 3(i,iii) since $\tau < 0$ and $1 < \gamma(q\kappa + \tau)$ for the denominator, by Lemma 1(ii) under
$\gamma((2 + \delta)\kappa + \tau) < 1$ for the first numerator and Lemma 1(v) for the second numerator, one may
write, since $d < +\infty$,

$$
\begin{align*}
\operatorname{cov}(\chi_{i_1,n}^{(1)}, \chi_{i_2,n}^{(3)}) &\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \left\{ \frac{g(y_n)(\lambda \bar{F})^{\frac{2}{2+\delta}}(y_n)}{g(y_n)(\lambda \bar{F})^2(y_n)} + \frac{(\lambda \bar{F})^{\frac{2}{2+\delta} - \frac{2}{q}}(y_n)}{g(y_n)(\lambda \bar{F})^2(y_n)} \right\} \\
&= \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \{(\lambda \bar{F})^{\frac{2}{2+\delta} - 2}(y_n) + g^{-1}(y_n)(\lambda \bar{F})^{\frac{2}{2+\delta} - \frac{2}{q} - 2}(y_n)\} \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) (\lambda \bar{F})^{\frac{2}{2+\delta} - 2}(y_n), \tag{22}
\end{align*}
$$
where we used $1 < \gamma(q\kappa/2 + \tau)$ for negligibility in the last line. By (Doukhan, 1994, Theorem 3),

$$
\begin{align*}
\operatorname{cov}(\chi_{i_1,n}^{(1)}, \chi_{i_2,n}^{(2)}) &= \operatorname{cov} \left( \sum_{j=1}^{d} a_{1,j} \left( \frac{\Lambda_{i_1}^{(j)} X_{i_1}^{(j)}}{m_{\Lambda^{(j)}X^{(j)}(y_n)}} \mathbf{1}_{\{Y_{i_1} \ge y_n\}} - 1 \right), \sum_{j=1}^{d} a_{2,j} \left( \frac{Y_{i_2} \Lambda_{i_2}^{(j)} X_{i_2}^{(j)}}{m_{Y\Lambda^{(j)}X^{(j)}(y_n)}} \mathbf{1}_{\{Y_{i_2} \ge y_n\}} - 1 \right) \right) \\
&= \sum_{j_1,j_2=1}^{d} \alpha_{1,j_1} \alpha_{2,j_2} \frac{\operatorname{cov}(\Lambda_{i_1}^{(j_1)} X_{i_1}^{(j_1)} \mathbf{1}_{\{Y_{i_1} \ge y_n\}}, Y_{i_2} \Lambda_{i_2}^{(j_2)} X_{i_2}^{(j_2)} \mathbf{1}_{\{Y_{i_2} \ge y_n\}})}{m_{\Lambda^{(j_1)}X^{(j_1)}(y_n)} m_{Y\Lambda^{(j_2)}X^{(j_2)}(y_n)}} \\
&\le 8\alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \sum_{j_1,j_2=1}^{d} a_{1,j_1} a_{3,j_2} \prod_{s \in \{0,1\}} \frac{\mathbb{E}(Y^s \Lambda^{(j_{s+1})} (X^{(j_{s+1})})^{2+\delta} \mathbf{1}_{\{Y \ge y_n\}})^{\frac{1}{2+\delta}}}{m_{Y^s \Lambda^{(j_{s+1})} X^{(j_{s+1})}(y_n)}} \\
&=: 8\alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \sum_{j_1,j_2=1}^{d} a_{1,j_1} a_{3,j_2} \prod_{s \in \{0,1\}} E_{j_s, \delta, n}.
\end{align*}
$$

Now, in view of (M₁) and a Jensen-type argument, one can decompose:

$$
\prod_{s \in \{0,1\}} E_{j_s, \delta, n} \lesssim \frac{\{\mathbb{E}(g^{2+\delta}(Y)\Lambda^{(j_1)}\mathbf{1}_{\{Y \ge y_n\}}) \cdot \mathbb{E}(Yg^{2+\delta}(Y)\Lambda^{(j_2)}\mathbf{1}_{\{Y \ge y_n\}})\}^{\frac{1}{2+\delta}}}{m_{\Lambda^{(j_1)}X^{(j_1)}(y_n)} m_{Y\Lambda^{(j_2)}X^{(j_2)}(y_n)}} \\
+ \frac{\{\mathbb{E}(\Lambda^{(j_1)}(\varepsilon^{(j_1)})^{2+\delta}\mathbf{1}_{\{Y \ge y_n\}}) \cdot \mathbb{E}(Y\Lambda^{(j_2)}(\varepsilon^{(j_2)})^{2+\delta}\mathbf{1}_{\{Y \ge y_n\}})\}^{\frac{1}{2+\delta}}}{m_{\Lambda^{(j_1)}X^{(j_1)}(y_n)} m_{Y\Lambda^{(j_2)}X^{(j_2)}(y_n)}} \\
+ \frac{\{\mathbb{E}(\Lambda^{(j_1)}(\varepsilon^{(j_1)})^{2+\delta}\mathbf{1}_{\{Y \ge y_n\}}) \cdot \mathbb{E}(Yg^{2+\delta}(Y)\Lambda^{(j_2)}\mathbf{1}_{\{Y \ge y_n\}})\}^{\frac{1}{2+\delta}}}{m_{\Lambda^{(j_1)}X^{(j_1)}(y_n)} m_{Y\Lambda^{(j_2)}X^{(j_2)}(y_n)}} \\
+ \frac{\{\mathbb{E}(g^{2+\delta}(Y)\Lambda^{(j_1)}\mathbf{1}_{\{Y \ge y_n\}}) \cdot \mathbb{E}(Y\Lambda^{(j_2)}(\varepsilon^{(j_2)})^{2+\delta}\mathbf{1}_{\{Y \ge y_n\}})\}^{\frac{1}{2+\delta}}}{m_{\Lambda^{(j_1)}X^{(j_1)}(y_n)} m_{Y\Lambda^{(j_2)}X^{(j_2)}(y_n)}}.
$$

By Lemma 3(iii) under $\gamma(\kappa + \tau + 1) < 1$, the common denominator is, up to a multiplicative constant, asymptotically equal to $y_n g^2(y_n)(\lambda \bar{F})^2(y_n)$. The first term is then $O\left(y_n^{\frac{1}{2+\delta}-1}(\lambda \bar{F})^{\frac{2}{2+\delta}-2}(y_n)\right)$
by Lemma 1(ii) under $\gamma((2+\delta)\kappa+\tau+1) < 1$. The second term is $O\left(y_n^{\frac{1}{2+\delta}-1}g^{-2}(y_n)(\lambda \bar{F})^{\frac{2}{2+\delta}-\frac{2}{q}-2}(y_n)\right)$
thanks to Lemma 1(v) under $\gamma(\frac{q}{q-2-\delta} + \tau) < 1$. Using Lemma 1(ii,v), the third and fourth terms
are both $O\left(y_n^{\frac{1}{2+\delta}-1}g^{-1}(y_n)(\lambda \bar{F})^{\frac{2}{2+\delta}-\frac{1}{q}-2}(y_n)\right)$. Notice that the leading term is the first one under
$1 < \gamma(q\kappa + \tau)$ while being bounded by $(\lambda \bar{F})^{\frac{2}{2+\delta}-2}(y_n)$. Therefore, it follows that, since $d < +\infty$,

$$
\operatorname{cov}(\chi_{i_1,n}^{(1)}, \chi_{i_2,n}^{(2)}) \lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) (\lambda \bar{F})^{\frac{2}{2+\delta}-2} (y_n). \quad (23)
$$
Once more, by (Doukhan, 1994, Theorem 3), using (M₁) and a Jensen-type argument, it yields:

$$
\begin{align*}
\operatorname{cov}(\chi_{i_1,n}^{(3)}, \chi_{i_2,n}^{(2)}) &= \operatorname{cov} \left( \sum_{j=1}^{d} a_{3,j} \left( \frac{\Lambda_{i_1}^{(j)}}{m_{\Lambda^{(j)}(y_n)}} \mathbb{1}_{\{Y_{i_1} \ge y_n\}} - 1 \right), \sum_{j=1}^{d} a_{2,j} \left( \frac{Y_{i_2} \Lambda_{i_2}^{(j)} X_{i_2}^{(j)}}{m_{Y \Lambda^{(j)} X^{(j)}(y_n)}} \mathbb{1}_{\{Y_{i_2} \ge y_n\}} - 1 \right) \right) \\
&= \sum_{j_1,j_2=1}^{d} a_{3,j_1} a_{2,j_2} \frac{\operatorname{cov}(\Lambda_{i_1}^{(j_1)} \mathbb{1}_{\{Y_{i_1} \ge y_n\}}, Y_{i_2} \Lambda_{i_2}^{(j_2)} X_{i_2}^{(j_2)} \mathbb{1}_{\{Y_{i_2} \ge y_n\}})}{m_{\Lambda^{(j_1)}(y_n)} m_{Y \Lambda^{(j_2)} X^{(j_2)}(y_n)} \\
&\le 8 \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \sum_{j_1,j_2=1}^{d} a_{3,j_1} a_{2,j_2} \frac{m_{\Lambda^{(j_1)}(y_n)}^{\frac{1}{2+\delta}-1}}{m_{Y \Lambda^{(j_2)} X^{(j_2)}(y_n)}^{\frac{1}{2+\delta}} m_{Y \Lambda^{(j_2)}(X^{(j_2)})^{2+\delta}(y_n)} \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \sum_{j_1,j_2=1}^{d} a_{3,j_1} a_{2,j_2} \frac{m_{\Lambda^{(j_1)}(y_n)}^{\frac{1}{2+\delta}-1}}{m_{Y \Lambda^{(j_2)} X^{(j_2)}(y_n)}^{\frac{1}{2+\delta}} m_{Y \Lambda^{(j_2)} g^{2+\delta}(Y)(y_n)} + m_{Y \Lambda^{(j_2)}(\varepsilon^{(j_2)})^{2+\delta}(y_n)}
\end{align*}
$$

By Lemma 3(i) and Lemma 1(ii,v), since $d < +\infty$,

$$
\operatorname{cov}(\chi_{i_1,n}^{(3)}, \chi_{i_2,n}^{(2)}) \lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) y_n^{\frac{1}{2+\delta}-1} (\lambda \bar{F})^{\frac{2}{2+\delta}-2} (y_n) \{1 + g^{-1}(y_n) (\lambda \bar{F})^{-\frac{1}{q}}(y_n)\}.
$$

Since $1 < \gamma(\kappa q + \tau)$, we conclude that:

$$
\operatorname{cov}(\chi_{i_1,n}^{(3)}, \chi_{i_2,n}^{(2)}) \lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) (\lambda \bar{F})^{\frac{2}{2+\delta}-2} (y_n). \quad (24)
$$

It only remains to investigate the cases with $\chi_{i,n}^{(k)}$, $k \in \{4, 5\}$. This time, we employ (Hall et al., 2014, Corollary A.1) with herein $p = 2 + \delta$ to get, using (M₁) and Jensen's inequality,

$$
\begin{align*}
\operatorname{cov}(\chi_{i_1,n}^{(1)}, \chi_{i_2,n}^{(4)}) &= \operatorname{cov} \left( \sum_{j=1}^{d} a_{1,j} \left( \frac{\Lambda_{i_1}^{(j)} X_{i_1}^{(j)}}{m_{\Lambda^{(j)} X^{(j)}(y_n)}} \mathbb{1}_{\{Y_{i_1} \ge y_n\}} - 1 \right), a_4 \left( \frac{\mathbb{1}_{\{Y_{i_2} \ge y_n\}}}{\bar{F}(y_n)} - 1 \right) \right) \\
&= a_4 \sum_{j=1}^{d} a_{1,j} \frac{\operatorname{cov}(\Lambda_{i_1}^{(j_1)} X_{i_1}^{(j_1)} \mathbb{1}_{\{Y_{i_1} \ge y_n\}}, \mathbb{1}_{\{Y_{i_2} \ge y_n\}})}{m_{\Lambda^{(j_1)} X^{(j_1)}(y_n)} \bar{F}(y_n)} \\
&\le 6a_4 \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \sum_{j=1}^{d} a_{1,j} \frac{\mathbb{E}(\Lambda^{(j)}(X^{(j)})^{2+\delta} \mathbb{1}_{\{Y \ge y_n\}})^{\frac{1}{2+\delta}}}{m_{\Lambda^{(j_1)} X^{(j_1)}(y_n)}} \bar{F}^{-1}(y_n) \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \bar{F}^{-1}(y_n) \frac{m_{\Lambda^{(j)} g^{2+\delta}(Y)}^{\frac{1}{2+\delta}} (y_n) + m_{\Lambda^{(j)}(\varepsilon^{(j)})^{2+\delta}}^{\frac{1}{2+\delta}} (y_n)}{m_{\Lambda^{(j_1)} X^{(j_1)}(y_n)}}.
\end{align*}
$$

By Lemma 1(ii,v), Lemma 3(iii) and the fact that $1 < \gamma(\kappa q + \tau)$, it follows that,

$$
\begin{align*}
\operatorname{cov}(\chi_{i_1,n}^{(1)}, \chi_{i_2,n}^{(4)}) &\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \bar{F}^{-1}(y_n) \frac{g(y_n)(\lambda \bar{F})^{\frac{1}{2+\delta}}(y_n) + (\lambda \bar{F})^{\frac{1}{2+\delta}-\frac{1}{q}}(y_n)}{g(y_n)(\lambda \bar{F})(y_n)} \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \bar{F}^{-1}(y_n) (\lambda \bar{F})^{\frac{1}{2+\delta}-1}(y_n) \{1 + g^{-1}(\lambda \bar{F})^{-\frac{1}{q}}(y_n)\} \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \bar{F}^{-1}(y_n) (\lambda \bar{F})^{\frac{1}{2+\delta}-1}(y_n).
\end{align*}
$$
Again, by Lemma 1(i,ii,v) under $(2+\delta)\gamma < 1$, $\gamma((2+\delta)\kappa + \tau + 1) < 1$ and $\gamma(\frac{q}{q-2-\delta} + \tau) < 1$, for $\iota \in \{0,1\}$, we have that, using (M1) and Jensen's inequality:

$$
\begin{align*}
\operatorname{cov}(\chi_{i_1,n}^{(1+\iota)}, \chi_{i_2,n}^{(5)}) &= \operatorname{cov} \left( \sum_{j=1}^{d} a_{1+\iota,j} \left( \frac{Y_{i_1}^{\iota} \Lambda_{i_1}^{(j)} X_{i_1}^{(j)}}{m_{Y^{\iota} \Lambda^{(j)} X^{(j)}(y_n)} (Y_{i_1} \ge y_n)} - 1 \right), a_5 \left( \frac{Y_{i_2} 1_{\{Y_{i_2} \ge y_n\}}}{m_Y(y_n)} - 1 \right) \right) \\
&= a_5 \sum_{j=1}^{d} a_{1+\iota,j} \frac{\operatorname{cov}(Y_{i_1}^{\iota} \Lambda_{i_1}^{(j_1)} X_{i_1}^{(j_1)} 1_{\{Y_{i_1} \ge y_n\}}, Y_{i_2} 1_{\{Y_{i_2} \ge y_n\}})}{m_{Y^{\iota} \Lambda^{(j_1)} X^{(j_1)}(y_n)} m_Y(y_n)} \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \max_{1 \le j \le d} \frac{\operatorname{E}(Y^{\iota} \Lambda^{(j)} (X^{(j)})^{2+\delta} 1_{\{Y \ge y_n\}})^{\frac{1}{2+\delta}} m_{Y^{2+\delta}}^{\frac{1}{2+\delta}}(y_n)}{m_{Y^{\iota} \Lambda^{(j_1)} X^{(j_1)}(y_n)} m_Y(y_n)} \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \frac{m_{Y^{2+\delta}}^{\frac{1}{2+\delta}}(y_n)}{m_Y(y_n)} \frac{m_{Y^{\iota} \Lambda^{(j)} g^{2+\delta}(Y)}^{\frac{1}{2+\delta}}(y_n) + m_{Y^{\iota} \Lambda^{(j)} (\varepsilon^{(j)})^{2+\delta}}^{\frac{1}{2+\delta}}(y_n)}{m_{Y^{\iota} \Lambda^{(j_1)} X^{(j_1)}(y_n)}} \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \bar{F}^{\frac{1}{2+\delta}-1}(y_n) y_n^{\iota(\frac{1}{2+\delta}-1)} \frac{g(y_n)(\lambda \bar{F})^{\frac{1}{2+\delta}}(y_n) + (\lambda \bar{F})^{\frac{1}{2+\delta}-\frac{1}{q}}(y_n)}{g(y_n)(\lambda \bar{F})(y_n)} \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) y_n^{\iota(\frac{1}{2+\delta}-1)} (\lambda \bar{F}^2)^{\frac{1}{2+\delta}-1}(y_n).
\end{align*}
$$

Using in addition Lemma 3(i), we get:

$$
\begin{align*}
\operatorname{cov}(\chi_{i_1,n}^{(4)}, \chi_{i_2,n}^{(5)}) &= \operatorname{cov} \left( a_4 \left( \frac{\Lambda_{i_1}^{(j)}}{m_{\lambda^{(j)}}(y_n)} 1_{\{Y_{i_1} \ge y_n\}} - 1 \right), a_5 \left( \frac{Y_{i_2} 1_{\{Y_{i_2} \ge y_n\}}}{m_Y(y_n)} - 1 \right) \right) \\
&= a_4 a_5 \frac{\operatorname{cov}(\Lambda_{i_1}^{(j_1)} 1_{\{Y_{i_1} \ge y_n\}}, Y_{i_2} 1_{\{Y_{i_2} \ge y_n\}})}{m_{\lambda^{(j_1)}}(y_n) m_Y(y_n)} \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \frac{m_{Y^{2+\delta}}^{\frac{1}{2+\delta}}(y_n)}{m_Y(y_n)} m_{\lambda^{(j)}}^{\frac{1}{2+\delta}-1}(y_n) \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \bar{F}^{\frac{1}{2+\delta}-1}(y_n) (\lambda \bar{F})^{\frac{1}{2+\delta}-1}(y_n) \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) (\lambda \bar{F}^2)^{\frac{1}{2+\delta}-1}(y_n).
\end{align*}
$$

We may conclude that, since $\tau < 0$,

$$
\operatorname{cov}(\chi_{i_1,n}^{(1)}, \chi_{i_2,n}^{(4)}) \lor \operatorname{cov}(\chi_{i_1,n}^{(1)}, \chi_{i_2,n}^{(5)}) \lor \operatorname{cov}(\chi_{i_1,n}^{(4)}, \chi_{i_2,n}^{(5)}) \lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) (\lambda \bar{F})^{\frac{2}{2+\delta}-2}(y_n). \quad (25)
$$
Samewise, using (M₁) and Jensen's inequality,

$$
\begin{align*}
\operatorname{cov}(\chi_{i_1,n}^{(2)}, \chi_{i_2,n}^{(4)}) &= \operatorname{cov} \left( \sum_{j=1}^{d} a_{2,j} \left( \frac{Y_{i_1} \Lambda_{i_1}^{(j)} X_{i_1}^{(j)}}{m_{Y\Lambda^{(j)}X^{(j)}}(y_n)} \mathbb{1}_{\{Y_{i_1} \ge y_n\}} - 1 \right), a_4 \left( \frac{\mathbb{1}_{\{Y_{i_2} \ge y_n\}}}{\bar{F}(y_n)} - 1 \right) \right) \\
&= a_4 \sum_{j=1}^{d} a_{2,j} \frac{\operatorname{cov}(Y_{i_1} \Lambda_{i_1}^{(j)} X_{i_1}^{(j)} \mathbb{1}_{\{Y_{i_1} \ge y_n\}}, \mathbb{1}_{\{Y_{i_2} \ge y_n\}})}{m_{Y\Lambda^{(j)}X^{(j)}}(y_n) \bar{F}(y_n)} \\
&\le 6a_4 \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \sum_{j=1}^{d} a_{2,j} \frac{\mathbb{E}(Y\Lambda^{(j)}(X^{(j)})^{2+\delta} \mathbb{1}_{\{Y \ge y_n\}})^{\frac{1}{2+\delta}}}{m_{Y\Lambda^{(j)}X^{(j)}}(y_n)} \bar{F}^{-1}(y_n) \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \bar{F}^{-1}(y_n) \frac{m_{Y\Lambda^{(j)}X^{(j)}}^{\frac{1}{2+\delta}}(y_n) + m_{Y\Lambda^{(j)}(\varepsilon^{(j)})^{2+\delta}}^{\frac{1}{2+\delta}}(y_n)}{m_{Y\Lambda^{(j)}X^{(j)}}(y_n)} \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \bar{F}^{-1}(y_n) y_n^{\frac{1}{2+\delta}-1} \frac{g(y_n)(\lambda \bar{F})^{\frac{1}{2+\delta}}(y_n) + (\lambda \bar{F})^{\frac{1}{2+\delta}-\frac{1}{q}}(y_n)}{g(y_n)(\lambda \bar{F})(y_n)} \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \bar{F}^{-1}(y_n) y_n^{\frac{1}{2+\delta}-1} (\lambda \bar{F})^{\frac{1}{2+\delta}-1}(y_n) \{1 + g^{-1}(\lambda \bar{F})^{-\frac{1}{q}}(y_n)\} \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \bar{F}^{-1}(y_n) y_n^{\frac{1}{2+\delta}-1} (\lambda \bar{F})^{\frac{1}{2+\delta}-1}(y_n).
\end{align*}
$$

On the other hand, we have according to (Doukhan, 1994, Theorem 3),

$$
\begin{align*}
\operatorname{cov}(\chi_{i_1,n}^{(3)}, \chi_{i_2,n}^{(5)}) &= \operatorname{cov} \left( \sum_{j=1}^{d} a_{3,j} \left( \frac{\Lambda_{i_1}^{(j)}}{m_{\lambda^{(j)}}(y_n)} \mathbb{1}_{\{Y_{i_1} \ge y_n\}} - 1 \right), a_5 \left( \frac{Y_{i_2} \mathbb{1}_{\{Y_{i_2} \ge y_n\}}}{m_Y(y_n)} - 1 \right) \right) \\
&= a_5 \sum_{j=1}^{d} a_{3,j} \frac{\operatorname{cov}(\Lambda_{i_1}^{(j)} \mathbb{1}_{\{Y_{i_1} \ge y_n\}}, Y_{i_2} \mathbb{1}_{\{Y_{i_2} \ge y_n\}})}{m_{\lambda^{(j)}}(y_n) m_Y(y_n)} \\
&\lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \frac{m_{Y^{2+\delta}}^{\frac{1}{2+\delta}}(y_n)}{m_Y(y_n)} m_{\lambda^{(j)}}^{\frac{1}{2+\delta}-1}(y_n) \lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) (\lambda \bar{F}^2)^{\frac{1}{2+\delta}-1}(y_n).
\end{align*}
$$

Combining the previous results, it follows that, since $\tau < 0$,

$$
\operatorname{cov}(\chi_{i_1,n}^{(2)}, \chi_{i_2,n}^{(4)}) \lor \operatorname{cov}(\chi_{i_1,n}^{(3)}, \chi_{i_2,n}^{(5)}) \lor \operatorname{cov}(\chi_{i_1,n}^{(2)}, \chi_{i_2,n}^{(5)}) \lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) (\lambda \bar{F})^{\frac{2}{2+\delta}-2}(y_n). \quad (26)
$$

In the same fashion as before, (Hall et al., 2014, Corollary A.1) with herein $p = 2 + \delta$ yields:

$$
\begin{align*}
\operatorname{cov}(\chi_{i_1,n}^{(3)}, \chi_{i_2,n}^{(4)}) &= \operatorname{cov} \left( \sum_{j=1}^{d} a_{3,j} \left( \frac{\Lambda_{i_1}^{(j)}}{m_{\lambda^{(j)}}(y_n)} \mathbb{1}_{\{Y_{i_1} \ge y_n\}} - 1 \right), a_4 \left( \frac{\mathbb{1}_{\{Y_{i_2} \ge y_n\}}}{\bar{F}(y_n)} - 1 \right) \right) \\
&= a_4 \sum_{j=1}^{d} a_{3,j} \frac{\operatorname{cov}(\Lambda_{i_1}^{(j)} \mathbb{1}_{\{Y_{i_1} \ge y_n\}}, \mathbb{1}_{\{Y_{i_2} \ge y_n\}})}{m_{\lambda^{(j)}}(y_n) \bar{F}(y_n)} \\
&\le 6a_4 \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \bar{F}^{-1}(y_n) \sum_{j=1}^{d} a_{3,j} m_{\lambda^{(j)}}^{\frac{1}{2+\delta}-1}(y_n) \\
&\lesssim a_4 \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \bar{F}^{-1}(y_n) (\lambda \bar{F})^{\frac{1}{2+\delta}-1}(y_n).
\end{align*}
$$

Therefore, it yields since $\tau < 0$,

$$
\operatorname{cov}(\chi_{i_1,n}^{(3)}, \chi_{i_2,n}^{(4)}) \lesssim \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) (\lambda \bar{F})^{\frac{2}{2+\delta}-2}(y_n). \quad (27)
$$
Combining the bounds (19)–(27), it follows that, when $n$ is large enough and for $u_n \le n$,

$$ \sum_{1 \le i_1 < i_2 \le u_n} \operatorname{cov}(\chi_{i_1,n}, \chi_{i_2,n}) \lesssim \frac{1}{n} (\lambda \bar{F})^{\frac{2}{2+\delta}-2+\theta}(y_n) \sum_{1 \le i_1 < i_2 \le u_n} \alpha^{\frac{\delta}{2+\delta}}(|i_1 - i_2|). \quad (28) $$

Furthermore, the following inequality holds under (3),

$$ \sum_{1 \le i_1 < i_2 \le u_n} \alpha^{\frac{\delta}{2+\delta}}(i_2 - i_1) \lesssim u_n. \quad (29) $$

Indeed, a derivation writes as,

$$ \begin{aligned} \sum_{1 \le i_1 < i_2 \le u_n} \alpha^{\frac{\delta}{2+\delta}}(i_2 - i_1) &= \sum_{i=1}^{u_n} \sum_{j=1}^{u_n-i} \alpha^{\frac{\delta}{2+\delta}}(j) = \sum_{i=1}^{u_n} \sum_{j=1}^{u_n-i} \alpha^{\frac{\delta}{2+\delta}}(j) = \sum_{i,j} \alpha^{\frac{\delta}{2+\delta}}(j) \mathbf{1}_{\{1 \le i \le u_n, 1 \le j \le u_n - i\}} \\ &= \sum_{j=1}^{u_n-1} \sum_{i=1}^{u_n-j} \alpha^{\frac{\delta}{2+\delta}}(j) = \sum_{i,j} \alpha^{\frac{\delta}{2+\delta}}(j) \mathbf{1}_{\{1 \le i \le u_n - j, 1 \le j \le u_n - 1\}} \\ &= \sum_{j=1}^{u_n-1} (u_n - j) \alpha^{\frac{\delta}{2+\delta}}(j) \le u_n \sum_{j=1}^{u_n-1} \alpha^{\frac{\delta}{2+\delta}}(j). \end{aligned} $$

Once that the control of the variance (18) and of the covariance (28) are established, we may go back to the conditions (6)–(9) required for applying (Davison et al., 2023, Lemma C.7). Firstly, in the case $\theta = \theta_0$, the asymptotic (18) implies that $n \operatorname{var}(\chi_{1,n}) \ll 1$. Moreover, the case $\delta = 0$, namely $\theta_0 = 1$, corresponds to $n \operatorname{var}(\chi_{1,n})$ being convergent.

Secondly, (3), (18) and (28) readily yield the first statement (6) and the third one (8) since $n - r_n|n/r_n| \le r_n$. The second statement (7) also follows from (3), (18), (28) and (29) with $u_n = r_n$ since we have shown that:

$$ \frac{n}{r_n} \sum_{1 \le i_1 < i_2 \le r_n} \operatorname{cov}(\chi_{i_1,n}, \chi_{i_2,n}) \lesssim (\lambda \bar{F})^{\frac{2}{2+\delta}-2+\theta}(y_n) = (\lambda \bar{F})^{\theta-\theta_0}(y_n) < \infty. $$

When $\theta > \theta_0 := 1 + \frac{\delta}{2+\delta}$, the latest bound in the previous display is converging to zero, since $\tau < 0$, which in turn implies that the limiting variance is zero. Furthermore, it is not mandatory to verify the Lyapunov condition (9) according to the statement of (Davison et al., 2023, Lemma C.7). As such, the proof in the case $\theta > \theta_0$ is finished.

It only remains to fulfill the Lyapunov condition (9) in the case $\theta = \theta_0$. Let us start with the triangular inequality for the $L^{2+\delta}(\Omega)$-norm and Jensen's inequality to write

$$ \mathbb{E} \left( \left| \sum_{i=1}^{r_n} \chi_{i,n} \right|^{2+\delta} \right) = \left\| \sum_{i=1}^{r_n} \chi_{i,n} \right\|_{L^{2+\delta}(\Omega)}^{2+\delta} \le \sum_{i=1}^{r_n} \left\| \chi_{i,n} \right\|_{L^{2+\delta}(\Omega)}^{2+\delta} = \frac{(\lambda \bar{F})^{\theta}(y_n)}{n} \sum_{i=1}^{r_n} \sum_{k=1}^{5} \left\| \chi_{i,n}^{(k)} \right\|_{L^{2+\delta}(\Omega)}^{2+\delta}. $$
Denote in the sequel $\rho_n = (\lambda \bar{F})^\theta(y_n)/n$. Therefore, the homogeneity of the norm and the triangular inequality entail, for $n$ large enough,

$$
\begin{align*}
& \mathbb{E}\left(\left|\sum_{i=1}^{r_n} \chi_{i,n}\right|^{2+\delta}\right) \\
& \le \rho_n^{1+\delta/2} \sum_{i=1}^{r_n} \sum_{j=1}^{d} \left( a_{1,j}^{2+\delta} \left\| \frac{\Lambda_i^{(j)} X_i^{(j)}}{m_{\Lambda^{(j)}X^{(j)}}(y_n)} \mathbb{1}_{\{Y_i \ge y_n\}} \right\|_{L^{2+\delta}(\Omega)}^{2+\delta} + a_{2,j}^{2+\delta} \left\| \frac{Y_i \Lambda_i^{(j)} X_i^{(j)}}{m_{Y\Lambda^{(j)}X^{(j)}}(y_n)} \mathbb{1}_{\{Y_i \ge y_n\}} \right\|_{L^{2+\delta}(\Omega)}^{2+\delta} \right) \\
& \quad + \rho_n^{1+\delta/2} \sum_{i=1}^{r_n} a_5^{2+\delta} \left\| \frac{Y_i \mathbb{1}_{\{Y_i \ge y_n\}}}{m_Y(y_n)} \right\|_{L^{2+\delta}(\Omega)}^{2+\delta} + \rho_n^{1+\delta/2} \sum_{i=1}^{r_n} \left( 2a_4^{2+\delta} + 2a_5^{2+\delta} + \sum_{j=1}^{d} (a_{1,j}^{2+\delta} + 2a_{3,j}^{2+\delta} + a_{2,j}^{2+\delta}) \right) \\
& \lesssim \rho_n^{1+\delta/2} \sum_{i=1}^{r_n} \sum_{j=1}^{d} \left( a_{1,j}^{2+\delta} \frac{m_{\Lambda^{(j)}(X^{(j)})2+\delta}(y_n)}{m_{\Lambda^{(j)}X^{(j)}}^{2+\delta}(y_n)} + a_{2,j}^{2+\delta} \frac{m_{\Lambda^{(j)}(YX^{(j)})2+\delta}(y_n)}{m_{Y\Lambda^{(j)}X^{(j)}}^{2+\delta}(y_n)} \right) + \rho_n^{1+\delta/2} \sum_{i=1}^{r_n} a_5^{2+\delta} \frac{m_{Y^{2+\delta}}(y_n)}{m_{Y^{2+\delta}}^{2+\delta}(y_n)} + r_n \rho_n^{1+\delta/2} \\
& \lesssim \rho_n^{1+\delta/2} \sum_{i=1}^{r_n} \sum_{j=1}^{d} a_{1,j}^{2+\delta} \frac{(m_{g^{2+\delta}(Y_i)\Lambda^{(j)}} + m_{\Lambda^{(j)}(\varepsilon^{(j)})2+\delta})(y_n)}{m_{\Lambda^{(j)}X^{(j)}}^{2+\delta}(y_n)} \\
& \quad + \rho_n^{1+\delta/2} \sum_{i=1}^{r_n} \sum_{j=1}^{d} a_{2,j}^{2+\delta} \frac{(m_{Y^{2+\delta}g^{2+\delta}(Y)\Lambda^{(j)}} + m_{\Lambda^{(j)}(Y\varepsilon^{(j)})2+\delta})(y_n)}{m_{Y\Lambda^{(j)}X^{(j)}}^{2+\delta}(y_n)} + r_n \rho_n^{1+\delta/2} \bar{F}^{-(1+\delta)}(y_n) + r_n \rho_n^{1+\delta/2}.
\end{align*}
$$

By Lemma 1(ii) under $\gamma((2+\delta)(\kappa+1)+\tau) < 1$ and Lemma 1(v) under $\gamma(\frac{q}{q-2-\delta}(2+\delta)+\tau) < 1$ for the numerators, and by Lemma 3(i,iii) for the denominators, it follows that

$$
\begin{align*}
\mathbb{E}\left(\left|\sum_{i=1}^{r_n} \chi_{i,n}\right|^{2+\delta}\right) & \lesssim r_n \rho_n^{1+\delta/2} \left\{ (\lambda \bar{F})^{-(1+\delta)}(y_n) + g^{-(2+\delta)}(y_n) (\lambda \bar{F})^{-(1+\delta)-(2+\delta)/q}(y_n) \right\} \\
& \lesssim r_n \rho_n^{1+\delta/2} (\lambda \bar{F})^{-(1+\delta)}(y_n) = n^{-\delta/2} \ll 1,
\end{align*}
$$

where we used in the last line the fact that $1 < \gamma(\kappa q + \tau)$ and $\theta_0(1+\delta/2) = 1+\delta$. This concludes (9) and thus the proof.

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

**Proof of Theorem 1.** We start by considering the expansion $\hat{\beta}_{\Lambda}(y_n) - \beta = (\hat{\beta}_{\Lambda}(y_n) - w(y_n)) + (w(y_n) - \beta)$ where $w$ is given in (2). Then, by triangle inequality and in view of $g(y_n) \bar{F}^{1/q}(y_n) \|w(y_n) - \beta\| = O(1)$ provided by (Bousebata et al., 2023, Proposition 2) under $\gamma(\kappa+1) < 1$, it is sufficient to show that, for any $\theta > 1 + \frac{\delta}{2+\delta}$, for some Gaussian random vector $G$,

$$ n^{1/2}(\lambda \bar{F})^{\theta/2}(y_n) \left( \hat{\beta}_{\Lambda}(y_n) - w(y_n) \right) \xrightarrow[n \to +\infty]{\mathbb{P}} \beta \odot G, \quad \text{in } \mathbb{R}^p. \quad (30) $$

Let us introduce

$$ \psi_1 := \frac{1 - \gamma\tau}{1 - \gamma(\kappa + \tau + 1)} - \frac{1 - \gamma\tau}{(1 - \gamma)(1 - \gamma(\kappa + \tau))} \quad \text{and} \quad \psi_2 := \frac{\gamma^2\kappa}{(1 - \gamma(\kappa + 1))(1 - \gamma\kappa)(1 - \gamma)}. $$

Let $j \in \{1, \dots, d\}$. Applying the Delta-method together with Proposition 2 and the identity given in (4), there exists a sequence of real random variables $(\vartheta_{j,n})_j$ such that $n^{1/2}(\lambda \bar{F})^{\theta/2}(y_n)(\vartheta_{j,n})_j$ is asymptotically a centered Gaussian vector $G \in \mathbb{R}^d$ when $\theta = 1 + \frac{\delta}{2+\delta}$; in the case $\theta > 1 + \frac{\delta}{2+\delta}$, the variance of $G$ is zero, i.e., $n^{1/2}(\lambda \bar{F})^{\theta/2}(y_n)\vartheta_{j,n} \ll 1$ in probability. Additionally,

$$ \hat{v}_{\Lambda}^{(j)}(y_n) = \frac{\bar{F}(y_n)}{m_{\lambda^{(j)}(y_n)}} \left( \bar{F}(y_n) m_{Y\Lambda^{(j)}X^{(j)}}(y_n) - m_Y(y_n) m_{\Lambda^{(j)}X^{(j)}}(y_n) \right) (1 + \vartheta_{j,n}). $$
Next, using Lemma 1(i), Lemma 3(i,iii) and Lemma 3(ii) under $\gamma(\kappa + 1) < 1$, one obtains, since $\beta_j \neq 0$ when $j \in \{1, \dots, d\}$:

$$ \hat{v}_{\Lambda}^{(j)}(y_n) = \psi_1 \beta_j y_n g(y_n) \bar{F}(y_n) (1 + \vartheta_{j,n}). $$

Furthermore, going back to (4) and expanding $\|\hat{v}_{\Lambda}(y_n)\|^2 = \hat{v}_{\Lambda}(y_n)^{\top} \hat{v}_{\Lambda}(y_n)$ thanks to Proposition 2, Lemma 1(i) and Lemma 3(i,iii) yield

$$ \|\hat{v}_{\Lambda}(y_n)\| = \psi_1 y_n g(y_n) \bar{F}(y_n) (1 + o_{\mathbb{P}}(1)), \quad (31) $$

since $\beta$ is of unit norm, Besides, by Lemma 1(i) and Lemma 3(ii) under $\gamma(\kappa + 1) < 1$,

$$ v^{(j)}(y_n) = \psi_2 \beta_j y_n g(y_n) \bar{F}^2(y_n) (1 + o(1)) \text{ and } \|v(y_n)\| = \psi_2 y_n g(y_n) \bar{F}^2(y_n) (1 + o(1)). $$

Combining the last three displays, it yields:

$$ n^{1/2}(\lambda \bar{F})^{\theta/2}(y_n) \left( \frac{\hat{v}_{\Lambda}^{(j)}(y_n)}{\|\hat{v}_{\Lambda}(y_n)\|} - \frac{v^{(j)}(y_n)}{\|v(y_n)\|} \right)_{1 \le j \le d} = n^{1/2}(\lambda \bar{F})^{\theta/2}(y_n) (\beta_j \vartheta_{j,n})_{1 \le j \le d} (1 + o_{\mathbb{P}}(1)). $$

Now, let $d+1 \le j \le p$, meaning that $\beta_j = 0$. It remains to show that for any such $j$ and any $\theta \ge 1 + \frac{\delta}{2+\delta}$,

$$ n^{1/2}(\lambda \bar{F})^{\theta/2}(y_n) \frac{\hat{v}_{\Lambda}^{(j)}(y_n)}{\|\hat{v}_{\Lambda}(y_n)\|} \xrightarrow[n \to +\infty]{\mathbb{P}} 0. \quad (32) $$

When $\beta_j = 0$, using Proposition 2 with (4), one may write:

$$ \hat{v}_{\Lambda}^{(j)}(y_n) = \frac{1}{m_{\Lambda^{(j)}}(y_n)} \left( \bar{F}(y_n) \hat{m}_{Y_{\Lambda^{(j)}} \varepsilon^{(j)}}(y_n) - \hat{m}_Y(y_n) \hat{m}_{\Lambda^{(j)} \varepsilon^{(j)}}(y_n) \right) (1 + O_{\mathbb{P}}\{n(\lambda \bar{F})^{\theta}(y_n)\}^{1/2}). $$

The tool to control the random terms in the latter display will be Tchebychev's inequality. As such, we begin by studying the variance of each quantities. First, $\mathrm{var}(\hat{m}_{Y^t \Lambda^{(j)} \varepsilon^{(j)}}(y_n)) = \mathrm{var}(\frac{1}{n} \sum_{i=1}^n Y_i^t \Lambda^{(j)} \varepsilon^{(j)} 1_{\{Y_i \ge y_n\}})$ may be decomposed as

$$ \frac{1}{n} \mathrm{var}(Y^t \Lambda^{(j)} \varepsilon^{(j)} 1_{\{Y \ge y_n\}}) + \frac{2}{n^2} \sum_{i_1 < i_2} \mathrm{cov}(Y_{i_1}^t \Lambda_{i_1}^{(j)} \varepsilon_{i_1}^{(j)} 1_{\{Y_{i_1} \ge y_n\}}, Y_{i_2}^t \Lambda_{i_2}^{(j)} \varepsilon_{i_2}^{(j)} 1_{\{Y_{i_2} \ge y_n\}}). $$

By (Doukhan, 1994, Theorem 3) and Lemma 1(v) under $\gamma(\frac{q}{q-2-\delta}(2+\delta) + \tau) < 1$, we may bound the double sum of covariances and the variance as follows, using (29),

$$ \sum_{i_1 < i_2} \alpha^{\frac{\delta}{2+\delta}} (|i_1 - i_2|) \mathbb{E}(\Lambda^{(j)} (Y^t \varepsilon^{(j)})^{2+\delta} 1_{\{Y \ge y_n\}})^{\frac{2}{2+\delta}} \lesssim n y_n^{2t} (\lambda \bar{F})^{\frac{2}{2+\delta} - \frac{2}{q}} (y_n), \\ \mathrm{var}(Y^t \Lambda^{(j)} \varepsilon^{(j)} 1_{\{Y \ge y\}}) \le m_{\Lambda^{(j)}(Y^t \varepsilon^{(j)})} (y_n) \lesssim y_n^{2t} (\lambda \bar{F})^{1-2/q} (y_n). $$

Thus, by comparison, one gets:

$$ \mathrm{var}(\hat{m}_{Y^t \Lambda^{(j)} \varepsilon^{(j)}}(y_n)) \lesssim \frac{1}{n} y_n^{2t} (\lambda \bar{F})^{\frac{2}{2+\delta} - \frac{2}{q}} (y_n). $$

Finally, using Tchebychev's inequality together with $\mathbb{E}(\hat{m}_{Y^t \Lambda^{(j)} \varepsilon^{(j)}}(y_n)) = \mathbb{E}(Y^t \Lambda^{(j)} \varepsilon^{(j)} 1_{\{Y \ge y_n\}})$ upon which we apply Lemma 1(iv), and $n(\lambda \bar{F})^{\theta}(y_n) \gg 1$ for the comparison, it follows

$$ \hat{m}_{Y^t \Lambda^{(j)} \varepsilon^{(j)}}(y_n) = \mathbb{E}(Y^t \Lambda^{(j)} \varepsilon^{(j)} 1_{\{Y \ge y_n\}}) + O\left(\frac{1}{\sqrt{n}} y_n^t (\lambda \bar{F})^{\frac{1}{2+\delta} - \frac{1}{q}} (y_n)\right) \\ = O\left(y_n^t (\lambda \bar{F})^{1-1/q} (y_n)\right). $$
Besides, another application of Tchebychev's inequality, together with a similar argument as above leading to $\text{var}(\hat{m}_Y(y_n)) \lesssim y_n^{2+\delta} \bar{F}^{\frac{2}{2+\delta}}(y_n)$, gives after comparison,

$$ \hat{m}_Y(y_n) = O(y_n \bar{F}(y_n)) + O\left(\frac{1}{\sqrt{n}} y_n^{1+\frac{\delta}{2}} \bar{F}^{\frac{1}{2+\delta}}(y_n)\right) = O(y_n \bar{F}(y_n)). $$

Therefore, under (M₁), plugging the above asymptotics in (4) entails, using also (31) and Lemma 3(iii),

$$ n^{1/2}(\lambda \bar{F})^{\theta/2}(y_n) \frac{\hat{v}_{\Lambda}^{(j)}(y_n)}{\|\hat{v}_{\Lambda}(y_n)\|} = n^{1/2}(\lambda \bar{F})^{\theta/2}(y_n) O_{\mathbb{P}}\left(g^{-1}(y_n)(\lambda \bar{F})^{-1/q}(y_n)\right). $$

Note that the quantity in the brackets converges to zero under $1 < \gamma(\kappa q + \tau)$. The assumption that $ng^{-2}(y_n)(\lambda \bar{F})(y_n)^{\theta-2/q} \ll 1$ concludes the proof in the case $j \in \{d+1, \dots, p\}$, regardless of $\theta = 1 + \frac{\delta}{2+\delta}$ or not. □

## REFERENCES

Aghbalou, A., F. Portier, A. Sabourin, and C. Zhou (2024). Tail inverse regression: Dimension reduction for prediction of extremes. *Bernoulli* 30(1), 503–533.

Arbel, J., S. Girard, and H. Lorenzo (2024). Shrinkage for Extreme Partial Least Squares. *Statistics and Computing* 34, 181.

Basrak, B., R. A. Davis, and T. Mikosch (2002). Regular variation of GARCH processes. *Stochastic Processes and their Applications* 99(1), 95–115.

Beirlant, J., M. Bladt, G. Maribe, and A. Verster (2023). Estimation of tail parameters with missing largest observations. *Electronic Journal of Statistics* 17(2), 3728–3761.

Beirlant, J., Y. Goegebeur, J. Segers, and J. L. Teugels (2004). *Statistics of extremes: theory and applications*, Volume 558. New-York: John Wiley & Sons.

Bernard-Michel, C., L. Gardes, and S. Girard (2009). Gaussian regularized sliced inverse regression. *Statistics and Computing* 19(1), 85–98.

Bingham, N., C. Goldie, and J. Teugels (1987). *Regular Variation*, Volume 27 of *Encyclopedia of Mathematics and its application*. Cambridge University Press.

Bousebata, M., G. Enjolras, and S. Girard (2023). Extreme Partial Least-Squares. *Journal of Multivariate Analysis* 194, 105101.

Bradley, R. C. (2005). Basic properties of strong mixing conditions. a survey and some open questions. *Probability Surveys* 2, 107–144.

Cai, J.-J., J. H. Einmahl, L. de Haan, and C. Zhou (2015). Estimation of the marginal expected shortfall: the mean when a related variable is extreme. *Journal of the Royal Statistical Society, Series B* 77(2), 417–442.

Carrasco, M. and X. Chen (2002). Mixing and moment properties of various GARCH and stochastic volatility models. *Econometric Theory* 18(1), 17–39.

Cook, R. D. (2007). Fisher lecture: Dimension reduction in regression. *Statistical Science* 22(1), 1–26.

Davison, A. C., S. A. Padoan, and G. Stupfler (2023). Tail risk inference via expectiles in heavy-tailed time series. *Journal of Business & Economic Statistics* 41(3), 876–889.

de Haan, L. and A. Ferreira (2007). *Extreme value theory: an introduction*. New-York: Springer Science & Business Media.

Doukhan, P. (1994). *Mixing*. New York: Springer.

Drees, H. and H. Rootzén (2010). Limit theorems for empirical processes of cluster functionals. *The Annals of Statistics* 38(4), 2145–2186.
Engle, R. F. and T. Bollerslev (1986). Modelling the persistence of conditional variances. *Econometric Reviews* 5(1), 1–50.

Gardes, L. (2018). Tail dimension reduction for extreme quantile estimation. *Extremes* 21(1), 57–95.

Girard, S., G. Stupfler, and A. Usseglio-Carleve (2021). Extreme conditional expectile estimation in heavy-tailed heteroscedastic regression models. *The Annals of Statistics* 49(6), 3358–3382.

Hall, P., C. Heyde, Z. Birnbaum, and E. Lukacs (2014). *Martingale Limit Theory and Its Application*. Communication and Behavior. Academic Press.

Ibragimov, I. A. (1962). Some limit theorems for stationary processes. *Theory of Probability & Its Applications* 7(4), 349–382.

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

Mayer, I., A. Sportisse, N. Tierney, N. Vialaneix, and J. Josse (2022). R-miss-tastic: a unified platform for missing values methods and workflows. *The R Journal* 14(2), 244–266.

Mokkadem, A. (1988). Mixing properties of ARMA processes. *Stochastic Processes and their Applications* 29(2), 309–315.

Nelson, D. B. (1990). Stationarity and persistence in the GARCH(1,1) model. *Econometric Theory* 6(3), 318–334.

Rubin, D. B. (1976). Inference and missing data. *Biometrika* 63(3), 581–592.

Saracco, J. (1997). An asymptotic theory for sliced inverse regression. *Communications in Statistics – Theory Methods* 26(9), 2141–2171.

Scarrott, C. and A. MacDonald (2012). A review of extreme value threshold estimation and uncertainty quantification. *REVSTAT-Statistical Journal* 10(1), 33–60.

van Dijk, D., T. Teräsvirta, and P. H. Franses (2002). Smooth transition autoregressive models - a survey of recent developments. *Econometric Reviews* 21(1), 1–47.

Wold, H. (1975). Soft modelling by latent variables: The non-linear iterative partial least squares (nipals) approach. *Journal of Applied Probability* 12(S1), 117–142.

Xu, H., R. Davis, and G. Samorodnitsky (2022). Handling missing extremes in tail estimation. *Extremes* 25(2), 199–227.

Xu, W., H. J. Wang, and D. Li (2022). Extreme quantile estimation based on the tail single-index model. *Statistica Sinica* 32(2), 893–914.

**Acknowledgments.** This work is supported by the French National Research Agency (ANR) in the framework of the Investissements d'Avenir Program (ANR-15-IDEX-02). S. Girard also acknowledges the support of the Chair Stress Test, Risk Management and Financial Steering, led by the French Ecole Polytechnique and its Foundation and sponsored by BNP Paribas.
<figure><img src="image_7.png" /><figcaption>(A) $\gamma = 0.1, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(B) $\gamma = 0.1, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(C) $\gamma = 0.1, \tau = -1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(D) $\gamma = 0.5, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>(E) $\gamma = 0.5, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(F) $\gamma = 0.5, \tau = -1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(G) $\gamma = 0.9, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(H) $\gamma = 0.9, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(I) $\gamma = 0.9, \tau = -1, \kappa = 0.5$.</figcaption></figure>

<caption>FIGURE 1. Simulation results on the inverse model without serial dependence (i.i.d. case), i.e., $(\phi_{\text{resp}}, \theta_{\text{resp}}) = (0,0)$ and $(\omega_{\text{noise}}, \alpha_{\text{noise}}, \beta_{\text{noise}}) = (1,0,0)$.</caption>
<figure><img src="image_16.png" /><figcaption>(A) $\gamma = 0.1, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(B) $\gamma = 0.1, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(C) $\gamma = 0.1, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(D) $\gamma = 0.4, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(E) $\gamma = 0.4, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(F) $\gamma = 0.4, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>(G) $\gamma = 0.7, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>(H) $\gamma = 0.7, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>(I) $\gamma = 0.7, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<caption>FIGURE 2. Simulation results with serial dependence 'classic ARMA-resp. + classic GARCH-noise', i.e., ($\phi_{\text{resp}}, \theta_{\text{resp}}$) = (0.8, -0.3) and ($\omega_{\text{noise}}, \alpha_{\text{noise}}, \beta_{\text{noise}}$) = (0.05, 0.1, 0.85).</caption>
<figure><img src="image_25.png" /><figcaption>(A) $\gamma = 0.1, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(B) $\gamma = 0.1, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(C) $\gamma = 0.1, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(D) $\gamma = 0.4, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>(E) $\gamma = 0.4, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(F) $\gamma = 0.4, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>(G) $\gamma = 0.7, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>(H) $\gamma = 0.7, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>(I) $\gamma = 0.7, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

FIGURE 3. Simulation results on the inverse model with serial dependence of type 'classic ARMA-resp. + IGARCH-like noise', i.e., ($\phi_{\text{resp}}$, $\theta_{\text{resp}}$) = (0.8, -0.3) and ($\omega_{\text{noise}}$, $\alpha_{\text{noise}}$, $\beta_{\text{noise}}$) = (0.05, 0.05, 0.94).
<figure><img src="image_34.png" /><figcaption>(A) $\gamma = 0.1, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_35.png" /><figcaption>(B) $\gamma = 0.1, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_36.png" /><figcaption>(C) $\gamma = 0.1, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>(D) $\gamma = 0.4, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_38.png" /><figcaption>(E) $\gamma = 0.4, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>(F) $\gamma = 0.4, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_40.png" /><figcaption>(G) $\gamma = 0.7, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>(H) $\gamma = 0.7, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>(I) $\gamma = 0.7, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<caption>FIGURE 4. Simulation results on the inverse model with serial dependence of type 'pathological ARMA-resp. + classic GARCH noise', i.e., ($\phi_{\text{resp}}$, $\theta_{\text{resp}}$) = (0.99, -0.98) and ($\omega_{\text{noise}}$, $\alpha_{\text{noise}}$, $\beta_{\text{noise}}$) = (0.05, 0.1, 0.85).</caption>
<figure><img src="image_43.png" /><figcaption>(A) $\gamma = 0.1, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_44.png" /><figcaption>(B) $\gamma = 0.1, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>(C) $\gamma = 0.1, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>(D) $\gamma = 0.4, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_47.png" /><figcaption>(E) $\gamma = 0.4, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_48.png" /><figcaption>(F) $\gamma = 0.4, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_49.png" /><figcaption>(G) $\gamma = 0.7, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>(H) $\gamma = 0.7, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>(I) $\gamma = 0.7, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

FIGURE 5. Simulation results on the inverse model with serial dependence of type 'pathological ARMA-resp. + IGARCH-like noise', i.e., ($\phi_{\text{resp}}$, $\theta_{\text{resp}}$) = (0.99, -0.98) and ($\omega_{\text{noise}}$, $\alpha_{\text{noise}}$, $\beta_{\text{noise}}$) = (0.05, 0.05, 0.94).
<figure><img src="image_52.png" /><figcaption>FIGURE 6. Simulation results on the inverse model with serial dependence of type standard GARCH-resp and standard ARMA-noise, i.e., ($\omega_{\text{resp}}$, $\alpha_{\text{resp}}$, $\beta_{\text{resp}}$) = (0.05, 0.1, 0.85) and ($\phi_{\text{noise}}$, $\theta_{\text{noise}}$) = (0.8, -0.3).</figcaption></figure>
<figure><img src="image_53.png" /><figcaption>(A) $\gamma = 0.1, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>(B) $\gamma = 0.1, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>(C) $\gamma = 0.1, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_56.png" /><figcaption>(D) $\gamma = 0.4, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>(E) $\gamma = 0.4, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_58.png" /><figcaption>(F) $\gamma = 0.4, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_59.png" /><figcaption>(G) $\gamma = 0.5, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_60.png" /><figcaption>(H) $\gamma = 0.5, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>(I) $\gamma = 0.5, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

FIGURE 7. Simulation results on the inverse model with serial dependence of type IGARCH-like response and standard ARMA-noise, i.e., ($\omega_{\text{resp}}$, $\alpha_{\text{resp}}$, $\beta_{\text{resp}}$) = (0.05, 0.05, 0.94) and ($\phi_{\text{noise}}$, $\theta_{\text{noise}}$) = (0.8, -0.3).
<figure><img src="image_62.png" /><figcaption>(A) $\gamma = 0.1, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_63.png" /><figcaption>(B) $\gamma = 0.1, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_64.png" /><figcaption>(C) $\gamma = 0.1, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_65.png" /><figcaption>(D) $\gamma = 0.4, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_66.png" /><figcaption>(E) $\gamma = 0.4, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_67.png" /><figcaption>(F) $\gamma = 0.4, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_68.png" /><figcaption>(G) $\gamma = 0.5, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_69.png" /><figcaption>(H) $\gamma = 0.5, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_70.png" /><figcaption>(I) $\gamma = 0.5, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<caption>FIGURE 8. Simulation results on the inverse model with serial dependence of type standard GARCH-response and pathological ARMA-noise, i.e., $(\omega_{\text{resp}}, \alpha_{\text{resp}}, \beta_{\text{resp}}) = (0.05, 0.1, 0.85)$ and $(\phi_{\text{noise}}, \theta_{\text{noise}}) = (0.99, -0.98)$.</caption>
<figure><img src="image_71.png" /><figcaption>(A) $\gamma = 0.1, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_72.png" /><figcaption>(B) $\gamma = 0.1, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_73.png" /><figcaption>(C) $\gamma = 0.1, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_74.png" /><figcaption>(D) $\gamma = 0.4, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_75.png" /><figcaption>(E) $\gamma = 0.4, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_76.png" /><figcaption>(F) $\gamma = 0.4, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_77.png" /><figcaption>(G) $\gamma = 0.5, \tau = -0.1, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_78.png" /><figcaption>(H) $\gamma = 0.5, \tau = -0.5, \kappa = 0.5$.</figcaption></figure>

<figure><img src="image_79.png" /><figcaption>(I) $\gamma = 0.5, \tau = -0.9, \kappa = 0.5$.</figcaption></figure>

<caption>FIGURE 9. Simulation results on the inverse model with serial dependence of type IGARCH-like response and pathological ARMA-noise, i.e., $(\omega_{\text{resp}}, \alpha_{\text{resp}}, \beta_{\text{resp}}) = (0.05, 0.05, 0.94)$ and $(\phi_{\text{noise}}, \theta_{\text{noise}}) = (0.99, -0.98)$.</caption>
<figure><img src="image_80.png" /><figcaption>FIGURE 10. Simulation results with serial dependence of type ESTAR for the response and standard GARCH for the noise.</figcaption></figure>

<figure><img src="image_81.png" /><figcaption>FIGURE 11. Simulation results with serial dependence of type ESTAR for the response and both cases pathological GARCH for the noise.</figcaption></figure>
<figure><img src="image_82.png" /><figcaption>FIGURE 12. Evolution of mean rank performance across threshold percentiles for the different methods on real data. Results are based on evaluation across 5,794 admissible weather station triplets. Lower values indicate better performance (closer to rank one). The random projection results are not displayed here, the mean rank being constant equal to six.</figcaption></figure>

<figure><img src="image_83.png" /><figcaption>FIGURE 13. Tail covariance achieved by EPLS (blue) on the triplet (USC00416892, US1TXMRR003, USC00410058) of weather stations compared to the competitors and a random projection baseline. The *x*-axis is the index *k* associated to the order statistic $Y_{n-k+1,n}$ used as threshold. The grey shaded region indicates the full variability range of the random baseline.</figcaption></figure>
<figure><img src="image_84.png" /><figcaption>FIGURE 14. Comparison on the triplet (USC00415766, US1TXTN0058, US1TXWA0042) of conditional extreme quantile ($\alpha = 0.99$) estimation across three projection methods. Each panel displays the projected point cloud ($\hat{\beta}^\top X_i, Y_i$) with the corresponding estimated conditional quantile curve.</figcaption></figure>