# Bias-reduced and variance-corrected asymptotic Gaussian inference about extreme expectiles

Abdelaati Daouia, Gilles Stupfler, Antoine Usseglio-Carleve

► To cite this version:

Abdelaati Daouia, Gilles Stupfler, Antoine Usseglio-Carleve. Bias-reduced and variance-corrected asymptotic Gaussian inference about extreme expectiles. 2023. [https://hal.science/hal-04097302v2](https://hal.science/hal-04097302v2)

HAL Id: hal-04097302

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

Preprint submitted on 11 Dec 2023

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_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Bias-reduced and variance-corrected asymptotic Gaussian inference about extreme expectiles

Abdelaati Daouia<sup>a</sup>, Gilles Stupfler<sup>b</sup> & Antoine Usseglio-Carleve<sup>c</sup>

<sup>a</sup> Toulouse School of Economics, University of Toulouse Capitole, France

ORCID: 0000-0003-2621-8860

<sup>b</sup> Univ Angers, CNRS, LAREMA, SFR MATHSTIC, F-49000 Angers, France

ORCID: 0000-0003-2497-9412

<sup>c</sup> Avignon Université, LMA UPR 2151, 84000 Avignon, France

ORCID: 0000-0002-8148-3758

## Abstract

The expectile is a prime candidate for being a standard risk measure in actuarial and financial contexts, for its ability to recover information about probabilities and typical behavior of extreme values, as well as its excellent axiomatic properties. A series of recent papers has focused on expectile estimation at extreme levels, with a view on gathering essential information about low-probability, high-impact events that are of most interest to risk managers. The obtention of accurate confidence intervals for extreme expectiles is paramount in any decision process in which they are involved, but actual inference on these tail risk measures is still a difficult question due to their least squares nature and sensitivity to tail heaviness. This article focuses on asymptotic Gaussian inference about tail expectiles in the challenging context of heavy-tailed observations. We use an in-depth analysis of the proofs of asymptotic normality results for two classes of extreme expectile estimators to derive bias-reduced and variance-corrected Gaussian confidence intervals. These, unlike previous attempts in the literature, are well-rooted in statistical theory and can accommodate underlying distributions that display a wide range of tail behaviors. A large-scale simulation study and three real data analyses confirm the versatility of the proposed technique.

**MSC 2010 subject classifications:** 62E20, 62G15, 62G20, 62G30, 62G32, 62P05

**Keywords:** Asymptotic normality, Bias correction, Expectiles, Extreme values, Heavy tails, Inference, Variance correction

## 1 Introduction

The problem of correctly estimating and inferring extreme risk, carried by low-probability high-impact events such as systemic financial crises and natural disasters, arises in a large range of applications such as insurance and finance, where it is crucial to correctly evaluate and manage the risk carried by a portfolio of claimants or stocks. Several risk measures have been fruitfully used in this context. Among these, expectiles have recently grown increasingly popular for a number of reasons, including the fact that they induce the only law-invariant, coherent (Artzner et al., 1999) and elicitable (Gneiting, 2011) risk measure, see Bellini et al. (2014) and Ziegel (2016). Consequently, a straightforward backtesting methodology can be developed for expectiles, which allows one to rank expectile forecasts by their accuracy (see Theorem 10 in Gneiting, 2011). This crucially hinges on their
formulation as minimizers of an asymmetric squared loss function (Newey and Powell, 1987), as follows:

$$ \xi_{\tau} = \arg \min_{\theta \in \mathbb{R}} \mathbb{E}(\eta_{\tau}(Y - \theta) - \eta_{\tau}(Y)), \quad (1.1) $$

where Y denotes the loss variable, assumed to have a finite first moment, $\tau \in (0, 1)$ is the asymmetry level, and $\eta_{\tau}(u) = |\tau - 1\{u \le 0\}|u^2$ is the so-called expectile check function (throughout, $1\{\cdot\}$ denotes the indicator function).

Expectile estimation and inference was first developed in Newey and Powell (1987) in the context of testing for homoskedasticity and conditional symmetry in linear regression problems. Recent contributions include deep asymptotic results for the estimation of central, non-tail expectiles of fixed order $\tau$, see for example Holzmann and Klar (2016) and Krätschmer and Zähle (2017). The problem of estimating extreme expectiles, whose order $\tau \uparrow 1$, is more difficult and has been studied only fairly recently, even though it constitutes the right framework for the assessment of extreme risk. The starting point for extreme expectile estimation appears to be a series of papers focusing on the challenging case when the underlying distribution has a heavy right tail (Daouia et al., 2018, 2019, 2020). In this series of articles, two classes of estimators are developed: the first class extrapolates to the far tail a Least Asymmetrically Weighted Squares (LAWS) estimator, obtained through the minimization of the empirical counterpart of problem (1.1), and the second extrapolates a quantile-based (or indirect) estimator whose construction rests upon a remarkable asymptotic proportionality relationship linking extreme expectiles to their quantile analogs.

It was noted later by Padoan and Stupfler (2022) that asymptotic Gaussian inference of extreme expectiles using these two classes of estimators was a difficult question, due to the fact that the asymptotic variances of the Gaussian limiting distributions of the estimators tend to provide a poor representation of the actual uncertainty in finite samples. A solution put forward in Padoan and Stupfler (2022) is to construct corrected confidence intervals that better approximate this uncertainty. However, because their built-in bias correction operates under the restriction that finite-sample bias due to the second order approximation does not dominate, there is no guarantee that they perform well in the wider heavy tail framework. It is now known that comprehensive bias correction is necessary if a reasonable degree of finite-sample accuracy of the extreme expectile estimators is to be ensured (Girard et al., 2022b). The variability of the bias-corrected versions introduced in Girard et al. (2022b) is, however, similarly hard to handle using straightforward plug-in estimators of the asymptotic variances arising in the Gaussian limiting distributions, and no satisfactory solution for inference is provided therein. More generally, since the precise form of the corrections of Padoan and Stupfler (2022) is motivated only by intensive simulations in certain models, a reverse engineering of their construction is difficult, thus making a potential extension to other types of distributional tails (such as light tails) impossible.

The contribution of the present paper is to develop corrected Gaussian confidence intervals for extreme expectiles of heavy-tailed distributions, with a rigorous theoretical foundation and whose coverage is close to the nominal confidence level even in moderately large samples. Our method essentially consists in carefully identifying, and then correcting, the approximation errors made in the proofs of the asymptotic normality of the extrapolated LAWS and quantile-based estimators. These approximation errors are typically due to either (i) using the asymptotic connection between extreme quantiles and expectiles while ignoring higher-order error terms, (ii) incorrectly neglecting correlations between two estimators when the asymptotic behavior of one of them dominates, or (iii) employing the delta-method for linearization purposes and, in doing so, incurring variance
distortions that are not accounted for. We provide successive corrections for each of these
types of errors, resulting in refined variance estimators that converge to the asymptotic
variances. In order to do so, we start by carefully approximating the finite-sample variabil-
ity of the standard (i.e. without bias reduction) versions of the LAWS and quantile-based
estimators. Then, we remark that the bias correction procedures of Girard et al. (2022b)
introduce further variability, and we design a simple but nonetheless very effective extra
correction of the asymptotic variance that is able to push the coverage of the resulting
asymptotic Gaussian confidence intervals close to the nominal rate irrespective of the
degree of tail heaviness in the underlying distribution.

Even though our corrections are technically involved and rely on a very careful investigation of the probabilistic behavior of the expectile estimators, they are conceptually very simple and result in confidence intervals that are computationally cheap. In particular, we avoid resorting to bootstrap, which is computationally expensive and known to be difficult to calibrate when the underlying distribution is heavy-tailed: Athreya (1987), Knight (1989) and Hall (1990) show that the traditional bootstrap is not consistent for the distribution of a sample mean, and Angus (1993) proves that it is not consistent either for the distribution of a sample minimum or maximum (i.e. an extreme value). Being extensions of the mean and extreme values, tail expectiles are similarly difficult to infer using resampling schemes. A potential consistent alternative is subsample bootstrap, but this has been shown to perform poorly both for inference about the mean (Hall and Jing, 1998) and about the tail index of heavy-tailed distributions (Guillou, 2000). Parametric bootstrap in the spirit of Cornea-Madeira and Davidson (2015) is another possibility, but it has not, to the best of our knowledge, been implemented for heavy-tailed distributions. The approach we propose bypasses all these difficulties in using bootstrap by staying within the familiar realm of asymptotic Gaussian inference.

The outline of this paper is the following. Section 2 spells out in detail our statistical
framework as well as the two classes of expectile estimators that we will focus on. We then
work out bias-reduced and variance-corrected asymptotic Gaussian confidence intervals
built on the LAWS and quantile-based estimators in Sections 3 and 4, respectively, and
show that they have asymptotically correct coverage. Section 5 examines their finite-
sample performance on simulated data and on three samples of real data from insurance
and finance. Appendix A describes the estimators of the second-order parameters used
in our implementation. Appendix B gives further details as to how our corrections are
calculated. Appendix C contains all necessary proofs and Appendix D provides extra
finite-sample results about our simulation study and real data analyses. Our methods are
implemented in the freely available R package *Expectrem*, which can be downloaded at
https://github.com/AntoineUC/Expectrem.

## 2 Statistical framework and inferential problem

Let $F: y \mapsto \mathbb{P}(Y \le y)$ and $\bar{F} = 1-F$ denote, respectively, the distribution and survival functions of the loss variable $Y$, whose large values represent extreme losses. We focus on heavy-tailed distributions that are commonplace in insurance and finance (Embrechts et al., 1997), which amounts to assuming that $\bar{F}$ is regularly varying with negative index in a neighborhood of $+\infty$, that is, there is $\gamma > 0$ such that $\bar{F}(ty)/\bar{F}(t) \to y^{-1/\gamma}$ as $t \to \infty$ for any $y > 0$. Then $\mathbb{E}|Y| < \infty$ provided $\mathbb{E}|\min(Y, 0)| < \infty$ and $\gamma < 1$, so that the expectile $\xi_\tau$ is well-defined by (1.1) for any $\tau \in (0, 1)$.

Let the data points $Y_1, \dots, Y_n$ be independent copies of $Y$, and $\tau_n, \tau'_n \uparrow 1$ denote high asymmetry levels such that $n(1-\tau_n) \to \infty$ and $(1-\tau'_n)/(1-\tau_n) \to 0$ as $n \to \infty$. It is
useful to think of $\tau'_n$ as the target expectile level satisfying $n(1 - \tau'_n) \to c < \infty$, so that the extreme expectile $\xi_{\tau'_n}$ is located in a region where very few or no data points lie, and of $\tau_n$ as an intermediate level, i.e. “extreme, but not too much”, so that $\xi_{\tau_n}$ can be estimated nonparametrically. Then, according to Daouia et al. (2018), an extreme expectile $\xi_{\tau'_n}$ can be estimated by, first, estimating an intermediate expectile $\xi_{\tau_n}$ (by, say, $\bar{\xi}_{\tau_n}$) and the tail index $\gamma$ (by, say, $\bar{\gamma}$), before combining them to obtain an extrapolated estimator of the form

$$ \bar{\xi}_{\tau'_n}^* = \left( \frac{1 - \tau'_n}{1 - \tau_n} \right)^{-\bar{\gamma}} \bar{\xi}_{\tau_n}. $$

A first approach to estimate $\xi_{\tau_n}$ is by using the empirical counterpart of the minimization problem (1.1), giving rise to the LAWS estimator $\hat{\xi}_{\tau_n} = \arg\min_{\theta \in \mathbb{R}} \sum_{i=1}^n \eta_{\tau_n}(Y_i - \theta)$. An alternative option is to use the asymptotic proportionality relationship $\xi_{\tau_n} \sim (\gamma^{-1} - 1)^{-\gamma} q_{\tau_n}$ (see Bellini and Di Bernardino, 2017), and to estimate $q_{\tau_n}$ therein by $\hat{q}_{\tau_n} = Y_{n-\lfloor n(1-\tau_n) \rfloor, n}$, where $Y_{1,n} \le Y_{2,n} \le \dots \le Y_{n,n}$ are the order statistics relative to $Y_1, \dots, Y_n$. Combined with the Hill estimator (Hill, 1975)

$$ \hat{\gamma}_{\tau_n}^H = \frac{1}{\lfloor n(1-\tau_n) \rfloor} \sum_{i=1}^{\lfloor n(1-\tau_n) \rfloor} \log Y_{n-i+1,n} - \log Y_{n-\lfloor n(1-\tau_n) \rfloor, n}, $$

this results in a quantile-based (or indirect) estimator $\tilde{\xi}_{\tau_n} = (1/\hat{\gamma}_{\tau_n}^H - 1)^{-\hat{\gamma}_{\tau_n}^H} \hat{q}_{\tau_n}$. An expectile-based alternative estimator of $\gamma$, pioneered by Girard et al. (2022a,b), is

$$ \hat{\gamma}_{\tau_n}^E = \left( 1 + \frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1 - \tau_n} \right)^{-1} \quad \text{where} \quad \hat{F}_n(x) = \frac{1}{n} \sum_{i=1}^n 1\{Y_i > x\}. $$

All in all, the present article will focus on inference about $\xi_{\tau'_n}$ based on the estimators

$$ \hat{\xi}_{\tau'_n}^* = \left( \frac{1 - \tau'_n}{1 - \tau_n} \right)^{-\hat{\gamma}_{\tau_n}^E} \hat{\xi}_{\tau_n} \qquad (2.1) $$

$$ \text{and } \tilde{\xi}_{\tau'_n}^* = \left( \frac{1 - \tau'_n}{1 - \tau_n} \right)^{-\hat{\gamma}_{\tau_n}^H} \tilde{\xi}_{\tau_n} = \left( \frac{1 - \tau'_n}{1 - \tau_n} \right)^{-\hat{\gamma}_{\tau_n}^H} (1/\hat{\gamma}_{\tau_n}^H - 1)^{-\hat{\gamma}_{\tau_n}^H} \hat{q}_{\tau_n}. \qquad (2.2) $$

Their asymptotic analysis requires quantifying the gap between the extremes of $Y$ and pure Pareto extremes through the following second-order regular variation condition.

$C_2(\gamma, \rho, A)$ There exist $\gamma > 0$, $\rho \le 0$ and a measurable auxiliary function $A$ having constant sign and converging to 0 at infinity such that

$$ \forall y > 0, \lim_{t \to \infty} \frac{1}{A(1/\bar{F}(t))} \left( \frac{\bar{F}(ty)}{\bar{F}(t)} - y^{-1/\gamma} \right) = y^{-1/\gamma} \frac{y^{\rho/\gamma} - 1}{\gamma \rho}. $$

Here and throughout the ratio $(y^a - 1)/a$ should be read as $\log y$ when $a = 0$. In this context it was shown in Daouia et al. (2018) and Girard et al. (2022b) that, if $\rho < 0$ and under the conditions $\lambda_1 = \lim_{n \to \infty} \sqrt{n(1 - \tau_n)}A((1 - \tau_n)^{-1}) \in \mathbb{R}$, $\lambda_2 = \lim_{n \to \infty} \sqrt{n(1 - \tau_n)}/q_{\tau_n} \in \mathbb{R}$ and $\sqrt{n(1 - \tau_n)}/\log((1 - \tau_n)/(1 - \tau'_n)) \to \infty$, one has

$$ \frac{\sqrt{n(1 - \tau_n)}}{\log((1 - \tau_n)/(1 - \tau'_n))} \log \left( \frac{\hat{\xi}_{\tau'_n}^*}{\hat{\xi}_{\tau'_n}} \right) \xrightarrow{d} \mathcal{N} \left( \frac{\gamma(\gamma^{-1} - 1)^{1-\rho}}{1 - \gamma - \rho} \lambda_1 + \gamma^2(\gamma^{-1} - 1)^{\gamma+1}\mathbb{E}(Y)\lambda_2, \frac{\gamma^3(1 - \gamma)}{1 - 2\gamma} \right) \qquad (2.3) $$
when $\mathbb{E}|\min(Y, 0)|^2 < \infty$ and $\gamma < 1/2$, and

$$ \frac{\sqrt{n(1-\tau_n)}}{\log((1-\tau_n)/(1-\tau'_n))} \log\left(\frac{\tilde{\xi}_{\tau'_n}^*}{\xi_{\tau'_n}}\right) \xrightarrow{d} \mathcal{N}\left(\frac{\lambda_1}{1-\rho}, \gamma^2\right). \quad (2.4) $$

The asymptotic distributions in (2.3) and (2.4) feature bias components due to the semi-parametric heavy tail framework. This has motivated the recent work of Girard et al. (2022b) on bias reduction procedures for extreme expectile estimation. They assume that $A(t) = b\gamma t^\rho$, for a certain constant $b \neq 0$ and $\rho < 0$; this condition is satisfied by most classical continuous heavy-tailed distributions, see Table 4. Then, given estimates $\bar{\rho}$ of $\rho$ and $\bar{b}$ of $b$ (we use here the estimators of $\rho$ and $b$ provided by the R package Expectrem$^1$, whose construction we recall in Appendix A), they first introduce a bias-reduced version of the extrapolated LAWS estimator having the form

$$ \hat{\xi}_{\tau'_n}^{\star, \text{BR}} = \left( \frac{1 - \tau'_n}{1 - \tau_n} \right)^{-\hat{\gamma}_{\tau_n}^{\text{E,BR}}} \hat{\xi}_{\tau_n} \times (1 + \hat{B}_{1,n})(1 + \hat{B}_{2,n})(1 + \hat{B}_{3,n}), \quad (2.5) $$

where $\hat{\gamma}_{\tau_n}^{\text{E,BR}}$ is a bias-reduced version of $\hat{\gamma}_{\tau_n}^{\text{E}}$ defined as

$$ \hat{\gamma}_{\tau_n}^{\text{E,BR}} = \left( 1 + \frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1 - \tau_n} \times \left( 1 - \frac{\bar{Y}_n}{\hat{\xi}_{\tau_n}} \right)^{-1} (2\tau_n - 1) \left( 1 + \frac{\bar{b}(\hat{F}_n(\hat{\xi}_{\tau_n}))^{-\bar{\rho}}}{1 - \hat{\gamma}_{\tau_n}^{\text{E,BR}} - \bar{\rho}} \right) \right)^{-1}, $$

with $\bar{Y}_n$ being the sample mean of $Y_1, \dots, Y_n$. [In Girard et al. (2022b), a bias-reduced version of the Hill estimator is used in place of $\hat{\gamma}_{\tau_n}^{\text{E}}$ in the right-hand side above; the R package Expectrem has been updated to reflect this change, motivated by the construction of an estimator which relies on expectiles as the main tool.] The quantities $\hat{B}_{j,n}$, meanwhile, are defined as follows:

$$ \hat{B}_{1,n} = \frac{((1 - \tau'_n)/(1 - \tau_n))^{-\bar{\rho}} - 1}{\bar{\rho}} \bar{b} \hat{\gamma}_{\tau_n}^{\text{E,BR}} (1 - \tau_n)^{-\bar{\rho}}, $$

$$ 1 + \hat{B}_{2,n} = (1 + \hat{r}(\tau_n))^{\hat{\gamma}_{\tau_n}^{\text{E,BR}}} \left( 1 + \frac{\frac{(1/\hat{\gamma}_{\tau_n}^{\text{E,BR}} - 1)^{-\bar{\rho}}}{(1+\hat{r}(\tau_n))^\bar{\rho}} - 1}{\bar{\rho}} \bar{b} \hat{\gamma}_{\tau_n}^{\text{E,BR}} (1 - \tau_n)^{-\bar{\rho}} \right)^{-1} $$

$$ \text{where } 1 + \hat{r}(\tau_n) = \left(1 - \frac{\bar{Y}_n}{\hat{\xi}_{\tau_n}}\right) \frac{1}{2\tau_n - 1} \left(1 + \frac{\bar{b}(\hat{F}_n(\hat{\xi}_{\tau_n}))^{-\bar{\rho}}}{1 - \hat{\gamma}_{\tau_n}^{\text{E,BR}} - \bar{\rho}}\right)^{-1}, $$

$$ \text{and } 1 + \hat{B}_{3,n} = (1 + \hat{r}^*(\tau'_n))^{-\hat{\gamma}_{\tau_n}^{\text{E,BR}}} \left( 1 + \frac{\frac{(1/\hat{\gamma}_{\tau_n}^{\text{E,BR}} - 1)^{-\bar{\rho}}}{(1+\hat{r}^*(\tau'_n))^\bar{\rho}} - 1}{\bar{\rho}} \bar{b} \hat{\gamma}_{\tau_n}^{\text{E,BR}} (1 - \tau'_n)^{-\bar{\rho}} \right) $$

$$ \text{where } 1 + \hat{r}^*(\tau'_n) = \left(1 - \frac{\bar{Y}_n}{\hat{\xi}_{\tau'_n}^*}\right) \frac{1}{2\tau'_n - 1} \left(1 + \frac{\bar{b}(1/\hat{\gamma}_{\tau_n}^{\text{E,BR}} - 1)^{-\bar{\rho}}}{1 - \hat{\gamma}_{\tau_n}^{\text{E,BR}} - \bar{\rho}} (1 - \tau'_n)^{-\bar{\rho}}\right)^{-1}. $$

Besides, a bias-reduced version of the indirect, extrapolated quantile-based estimator is given by

$$ \tilde{\xi}_{\tau'_n}^{\star, \text{BR}} = \left( \frac{1 - \tau'_n}{1 - \tau_n} \right)^{-\hat{\gamma}_{\tau_n}^{\text{H,BR}}} \left\{ (1/\hat{\gamma}_{\tau_n}^{\text{H,BR}} - 1)^{-\hat{\gamma}_{\tau_n}^{\text{H,BR}}} \hat{q}_{\tau_n} \right\} \times (1 + \tilde{B}_{1,n})(1 + \tilde{B}_{3,n}). \quad (2.6) $$

<footnote>$^1$Credit goes to B.G. Manjunath and F. Caeiro for the original R implementation in package evt0, which is unfortunately unavailable from CRAN as of April 2023.</footnote>
Here $\tilde{B}_{1,n}$ is defined as $\hat{B}_{1,n}$, but with

$$ \hat{\gamma}_{\tau_n}^{\mathrm{H,BR}} = \hat{\gamma}_{\tau_n}^{\mathrm{H}} \left( 1 - \frac{\bar{b}}{1 - \bar{\rho}} (1 - \tau_n)^{-\bar{\rho}} \right), $$

which is a bias-reduced version of the Hill estimator proposed by Caeiro et al. (2005), in place of $\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}$. Similarly for $\tilde{B}_{3,n}$, where in addition $\hat{r}^*(\tau_n')$ is replaced by $\tilde{r}^*(\tau_n')$, in which $\hat{\xi}_{\tau_n'}^*$ is replaced by $\tilde{\xi}_{\tau_n'}^*$. Then, $\hat{\xi}_{\tau_n}^{\mathrm{BR}}$ and $\tilde{\xi}_{\tau_n'}^{\mathrm{BR}}$ in (2.5) and (2.6) have the same asymptotic behavior as $\hat{\xi}_{\tau_n'}^*$ and $\tilde{\xi}_{\tau_n'}^*$ in (2.3) and (2.4), respectively, the only difference being that they are asymptotically unbiased, see Theorem 2 in Girard et al. (2022b). A Gaussian 100(1 - $\alpha$)% asymptotic confidence interval for $\xi_{\tau_n'}$ based on $\hat{\xi}_{\tau_n'}^{\mathrm{BR}}$ is then

$$ \begin{gathered} \hat{I}_{\tau_n'}^{(1)}(\alpha) = \left[ \hat{\xi}_{\tau_n'}^{\mathrm{BR}} \exp \left( \pm \frac{\log((1 - \tau_n)/(1 - \tau_n'))}{\sqrt{n(1 - \tau_n)}} \sqrt{\hat{s}_n^{2,\mathrm{BR}}} \times z_{1-\alpha/2} \right) \right] \\ \text{with } \hat{s}_n^{2,\mathrm{BR}} = \frac{(\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}})^3 (1 - \hat{\gamma}_{\tau_n}^{\mathrm{E,BR}})}{1 - 2\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}}, \end{gathered} \quad (2.7) $$

where $z_{1-\alpha/2}$ is the quantile of level $1 - \alpha/2$ for the standard Gaussian distribution. An alternative confidence interval based on $\hat{\xi}_{\tau_n'}^{\mathrm{BR}}$ is

$$ \begin{gathered} \hat{I}_{\tau_n'}^{(1)}(\alpha) = \left[ \hat{\xi}_{\tau_n'}^{\mathrm{BR}} \exp \left( \pm \frac{\log((1 - \tau_n)/(1 - \tau_n'))}{\sqrt{n(1 - \tau_n)}} \sqrt{\hat{\sigma}_n^{2,\mathrm{BR}}} \times z_{1-\alpha/2} \right) \right] \\ \text{with } \hat{\sigma}_n^{2,\mathrm{BR}} = (\hat{\gamma}_{\tau_n}^{\mathrm{H,BR}})^2. \end{gathered} \quad (2.8) $$

The rationale behind our inference methods is that the asymptotic results recalled thus far ignore the correlations between the intermediate expectile estimator and the tail index estimator, as well as the inherent variability of the intermediate expectile estimator which can be very substantial in the heavy right tail setup. The intervals $\hat{I}_{\tau_n'}^{(1)}(\alpha)$ and $\hat{I}_{\tau_n'}^{(1)}(\alpha)$ then tend to have poor coverage even in the ideal scenario when the underlying distribution is very close to the Pareto distribution, as shown in the middle panels of Figures 1 and 2 below. By comparing their left and right panels, it can also be inferred that the bias correction is reasonably effective at least for sufficiently large $\tau_n$, whereas the estimation of the variances of $(\sqrt{n(1 - \tau_n)}/\log((1 - \tau_n)/(1 - \tau_n')))\log(\hat{\xi}_{\tau_n'}^{\mathrm{BR}}/\xi_{\tau_n'})$ and $(\sqrt{n(1 - \tau_n)}/\log((1 - \tau_n)/(1 - \tau_n')))\log(\hat{\xi}_{\tau_n'}^{\mathrm{BR}}/\xi_{\tau_n'})$, by $\hat{s}_n^{2,\mathrm{BR}}$ and $\hat{\sigma}_n^{2,\mathrm{BR}}$ respectively, is a long way off the truth. It is on the accurate estimation of these variances that our contributions below will focus.

### 3 Bias-reduced and variance-corrected LAWS-based inference

Assume throughout this section that $\mathbb{E}|\min(Y, 0)|^2 < \infty$ and $\gamma < 1/2$. The basic idea to show the joint asymptotic normality of $\hat{\gamma}_{\tau_n}^{\mathrm{E}}$ and $\hat{\xi}_{\tau_n}$, which are the key building blocks in the estimator $\hat{\xi}_{\tau_n'}^{\mathrm{BR}}$, is to note that this amounts to showing the joint asymptotic normality of intermediate expectiles and quantiles at suitably chosen intermediate levels. The crucial point then is that the expectile $\xi_{\tau}$ is just the $\tau$th quantile of the distribution function $E = 1 - \bar{E}$ defined as

$$ \bar{E}(y) = \frac{\mathbb{E}(|Y - y|1\{Y > y\})}{\mathbb{E}(|Y - y|)} = \frac{\varphi^{(1)}(y)}{2\varphi^{(1)}(y) + y - \mathbb{E}(Y)}, \quad \text{with } \varphi^{(a)}(y) = \mathbb{E}([Y - y]^a 1\{Y > y\}). $$
<figure><img src="image_2.png" /><figcaption>Figure 1: Left panels: Comparison of the empirical variance of $(\sqrt{n(1-\tau_n)}/\log((1-\tau_n)/(1-\tau'_n))) \log(\hat{\xi}_{\tau'_n}^{\star,\text{BR}}/\xi_{\tau'_n})$ (red curve) with boxplots of the asymptotic variance estimate $\hat{s}_n^{2,\text{BR}}$. Middle panels: Empirical coverage probability of the associated 95% confidence interval $\hat{I}_{\tau'_n}^{(1)}(\alpha)$ with the nominal level $1-\alpha = 0.95$ in red. Right panels: Target value $\xi_{\tau'_n}$ (red line) and medians of the lower and upper bounds of $\hat{I}_{\tau'_n}^{(1)}(\alpha)$. These results were obtained from $N = 5,000$ replicated samples of $n = 1,000$ independent observations from the following distributions: the Burr distribution with parameters $\gamma = 0.2$ and $\rho = -5$ (first row) and the Fréchet distribution with parameter $\gamma = 0.4$ (second row). The target expectile level is $\tau'_n = 1 - 1/n = 0.999$, and the results are represented as functions of the discretized level $k = k_n = n(1-\tau_n)$ throughout.</figcaption></figure>

See Jones (1994). Since $\hat{\xi}_{\tau}$ is the $\tau$th expectile relative to the empirical distribution function $\hat{F}_n = 1 - \hat{F}_n$, it is the $\tau$th quantile of the distribution function $\hat{E}_n = 1 - \hat{E}_n$ defined as

$$ \hat{E}_n(y) = \frac{\hat{\varphi}_n^{(1)}(y)}{2\hat{\varphi}_n^{(1)}(y) + y - \bar{Y}_n}, \quad \text{with} \quad \hat{\varphi}_n^{(a)}(y) = \frac{1}{n} \sum_{i=1}^{n} (Y_i - y)^a \mathbb{1}\{Y_i > y\} \quad \text{and} \quad \bar{Y}_n = \frac{1}{n} \sum_{i=1}^{n} Y_i. $$

It follows that the joint asymptotic normality of intermediate expectiles and quantiles is equivalent to the joint asymptotic normality of $\hat{E}_n$ and $\hat{F}_n$ at suitably chosen intermediate levels. Our purpose is to finely quantify the joint uncertainty in this limiting relationship and account for the various sources of errors in the asymptotic approximations, as a preliminary to doing so again for the extrapolated version of $\hat{\xi}_{\tau'_n}^{\star,\text{BR}}$. We give a sketch of the most important steps below; full details, following from a careful inspection of the proofs of Theorem 2 in Stupfler and Usseglio-Carleve (2023) and Theorem 2 in Girard et al. (2022b), can be found in Appendix B.1. It is, first of all, straightforward to prove that if $\tau_n$ and $\alpha_n$ are intermediate levels (with $\alpha_n$ to be suitably chosen later), the covariance matrix
<figure><img src="image_3.png" /><figcaption>Figure 2: Left panels: Comparison of the empirical variance of $(\sqrt{n(1-\tau_n)}/\log((1-\tau_n)/(1-\tau'_n))) \log(\tilde{\xi}_{\tau'_n}^{\star,\text{BR}}/\xi_{\tau'_n})$ (red curve) with boxplots of the asymptotic variance estimate $\tilde{\sigma}_n^{2,\text{BR}}$. Middle panels: Empirical coverage probability of the associated 95% confidence interval $\tilde{I}_{\tau'_n}^{(1)}(\alpha)$ with the nominal level $1-\alpha = 0.95$ in red. Right panels: Target value $\xi_{\tau'_n}$ (red line) and medians of the lower and upper bounds of $\tilde{I}_{\tau'_n}^{(1)}(\alpha)$. These results were obtained from $N = 5,000$ replicated samples of $n = 1,000$ independent observations from the following distributions: the Student distribution with $\nu = 1/\gamma = 0.6$ degrees of freedom (first row) and the Generalized Pareto Distribution with shape parameter $\gamma = 0.8$ and unit scale (second row). The target expectile level is $\tau'_n = 1 - 1/n = 0.999$, and the results are represented as functions of the discretized level $k = k_n = n(1-\tau_n)$ throughout.</figcaption></figure>

of $\sqrt{n(1-\tau_n)}(\hat{\varphi}_n^{(1)}(\xi_{\tau_n})/\varphi^{(1)}(\xi_{\tau_n}) - 1, \hat{F}_n(q_{\alpha_n})/\bar{F}(q_{\alpha_n}) - 1)$ is exactly the 2 × 2 symmetric matrix $\mathbf{M}^{\varphi} = \mathbf{M}_n^{\varphi}$ having components:

$$ M_{n,11}^{\varphi} = (1 - \tau_n) \left( \frac{\varphi^{(2)}(\xi_{\tau_n})}{[\varphi^{(1)}(\xi_{\tau_n})]^2} - 1 \right), \quad M_{n,22}^{\varphi} = (1 - \tau_n) \frac{\alpha_n}{1 - \alpha_n} \\ \text{and } M_{n,12}^{\varphi} = (1 - \tau_n) \left( \frac{\varphi^{(1)}(\xi_{\tau_n} \vee q_{\alpha_n}) + (\xi_{\tau_n} \vee q_{\alpha_n} - \xi_{\tau_n}) \bar{F}(\xi_{\tau_n} \vee q_{\alpha_n})}{\varphi^{(1)}(\xi_{\tau_n})(1 - \alpha_n)} - 1 \right). $$

Then, since $\bar{Y}_n - \mathbb{E}(Y)$ converges to 0 at the rate $1/\sqrt{n}$ by the central limit theorem, and thus converges faster than extreme value procedures, we neglect the uncertainty in $\bar{Y}_n$ throughout and in particular in $\hat{E}_n$. A straightforward calculation provides

$$ \frac{\hat{E}_n(\xi_{\tau_n})}{\bar{E}(\xi_{\tau_n})} - 1 \approx \frac{\xi_{\tau_n} - \mathbb{E}(Y)}{2\varphi^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \mathbb{E}(Y)} \left( \frac{\hat{\varphi}_n^{(1)}(\xi_{\tau_n})}{\varphi^{(1)}(\xi_{\tau_n})} - 1 \right). $$

This suggests that the joint distribution of $\sqrt{n(1-\tau_n)}(\hat{E}_n(\xi_{\tau_n})/\bar{E}(\xi_{\tau_n}) - 1, \hat{F}_n(q_{\alpha_n})/\bar{F}(q_{\alpha_n}) - 1)$ has a covariance matrix that can be more accurately approximated by $\mathbf{M}^E = \mathbf{M}_n^E$ hav-
ing components:

$$M_{n,11}^{E} = \frac{(\xi_{\tau_n} - \mathbb{E}(Y))^2 M_{n,11}^{\varphi}}{(2\varphi^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \mathbb{E}(Y))^2}, \quad M_{n,12}^{E} = \frac{(\xi_{\tau_n} - \mathbb{E}(Y)) M_{n,12}^{\varphi}}{2\varphi^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \mathbb{E}(Y)}, \quad M_{n,22}^{E} = M_{n,22}^{\varphi}.$$

The next step is to remark that for $u, v \in \mathbb{R}$,

$$\begin{aligned} & \mathbb{P}\left(\left\{\sqrt{n(1-\tau_n)}\left(\frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1\right) \le u\right\} \cap \left\{\sqrt{n(1-\tau_n)}\left(\frac{\hat{q}_{\alpha_n}}{q_{\alpha_n}} - 1\right) \le v\right\}\right) \\ & \approx \mathbb{P}\left(\left\{\sqrt{n(1-\tau_n)}\left(\frac{\hat{E}_n(\xi_{\tau_n})}{\overline{E}(\xi_{\tau_n})} - 1\right) \le \sqrt{n(1-\tau_n)}\left(\frac{\overline{E}(\xi_{\tau_n})}{\overline{E}(y_n)} - 1\right)\right\} \right. \\ & \qquad \left. \cap \left\{\sqrt{n(1-\tau_n)}\left(\frac{\hat{F}_n(q_{\alpha_n})}{\overline{F}(q_{\alpha_n})} - 1\right) \le \sqrt{n(1-\tau_n)}\left(\frac{\overline{F}(q_{\alpha_n})}{\overline{F}(z_n)} - 1\right)\right\}\right) \end{aligned}$$

where $y_n = y_n(u) = \xi_{\tau_n}(1 + u/\sqrt{n(1-\tau_n)})$ and $z_n = z_n(v) = q_{\alpha_n}(1 + v/\sqrt{n(1-\tau_n)})$. A Taylor expansion yields

$$\frac{\overline{E}(\xi_{\tau_n})}{\overline{E}(y_n)} - 1 = u\xi_{\tau_n} \frac{1}{\sqrt{n(1-\tau_n)}} \frac{\varphi^{(1)}(\xi_{\tau_n}) + \overline{F}(\xi_{\tau_n})(\xi_{\tau_n} - \mathbb{E}(Y))}{\varphi^{(1)}(\xi_{\tau_n})(2\varphi^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \mathbb{E}(Y))(1 + o(1))}$$

and assumption $C_2(\gamma, \rho, A)$ entails $\sqrt{n(1-\tau_n)}(\overline{F}(q_{\alpha_n})/\overline{F}(z_n)-1) \to v/\gamma$. It readily follows that the covariance matrix of $\sqrt{n(1-\tau_n)}(\hat{\xi}_{\tau_n}/\xi_{\tau_n} - 1, \hat{q}_{\alpha_n}/q_{\alpha_n} - 1)$ is approximated by $\mathbf{M}^{\xi} = \mathbf{M}_{n,11}^{\xi}$, whose elements are

$$M_{n,11}^{\xi} = \left( \frac{\varphi^{(1)}(\xi_{\tau_n})(2\varphi^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \mathbb{E}(Y))}{\xi_{\tau_n}(\varphi^{(1)}(\xi_{\tau_n}) + \overline{F}(\xi_{\tau_n})(\xi_{\tau_n} - \mathbb{E}(Y)))} \right)^2 M_{n,11}^{E},$$

$$M_{n,12}^{\xi} = \gamma \left( \frac{\varphi^{(1)}(\xi_{\tau_n})(2\varphi^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \mathbb{E}(Y))}{\xi_{\tau_n}(\varphi^{(1)}(\xi_{\tau_n}) + \overline{F}(\xi_{\tau_n})(\xi_{\tau_n} - \mathbb{E}(Y)))} \right) M_{n,12}^{E}, \text{ and } M_{n,22}^{\xi} = \gamma^2 M_{n,22}^{E}.$$

Now that we have obtained a suitable approximation to the covariance matrix of the Gaussian asymptotic distribution of $\hat{\xi}_{\tau_n}$ and $\hat{q}_{\alpha_n}$, we recall that $\hat{\gamma}_{\tau_n}^E = (1 + \hat{F}_n(\hat{\xi}_{\tau_n})/(1 - \tau_n))^{-1}$ and we note that

$$\left\{ \sqrt{n(1-\tau_n)} \left( \frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1-\tau_n} - \frac{\overline{F}(\xi_{\tau_n})}{1-\tau_n} \right) \le z \right\} = \left\{ \hat{\xi}_{\tau_n} \ge \hat{q}_{1-(1-\tau_n)(\overline{F}(\xi_{\tau_n})/(1-\tau_n)+z/\sqrt{n(1-\tau_n}))} \right\}.$$

The intermediate quantile level $\beta_n = \beta_n(z) = 1 - (1 - \tau_n)(\overline{F}(\xi_{\tau_n})/(1 - \tau_n) + z/\sqrt{n(1 - \tau_n)})$ in the right-hand side is such that $1 - \beta_n \sim 1 - \alpha_n = \overline{F}(\xi_{\tau_n}) \sim (\gamma^{-1} - 1)(1 - \tau_n)$. Then

$$\begin{aligned} & \mathbb{P}\left(\left\{\sqrt{n(1-\tau_n)}\left(\frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1-\tau_n} - \frac{\overline{F}(\xi_{\tau_n})}{1-\tau_n}\right) \le z\right\} \cap \left\{\sqrt{n(1-\tau_n)}\left(\frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1\right) \le z'\right\}\right) \\ & \approx \mathbb{P}\left(\left\{\sqrt{n(1-\tau_n)}\left(\frac{\hat{q}_{\alpha_n}}{q_{\alpha_n}} - 1\right) \frac{q_{\alpha_n}}{\xi_{\tau_n}} - \sqrt{n(1-\tau_n)}\left(\frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1\right) \le -\sqrt{n(1-\tau_n)}\left(\frac{q_{\beta_n}}{\xi_{\tau_n}} - 1\right)\right\} \right. \\ & \qquad \left. \cap \left\{\sqrt{n(1-\tau_n)}\left(\frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1\right) \le z'\right\}\right). \end{aligned}$$
From Proposition 1 in Daouia et al. (2020) and the second-order regular variation assumption,

$$ \sqrt{n(1 - \tau_n)} \left( \frac{q_{\beta_n}}{\xi_{\tau_n}} - 1 \right) = \sqrt{n(1 - \tau_n)} \left( \frac{q_{\beta_n}}{q_{\tau_n}} \times \frac{q_{\tau_n}}{\xi_{\tau_n}} - 1 \right) \approx -\gamma \frac{1 - \tau_n}{\bar{F}(\xi_{\tau_n})} z $$

up to bias terms, which will be taken care of in the bias reduction procedure (and therefore can be neglected at this stage). Conclude that the asymptotic distribution of

$$ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1 - \tau_n} - \frac{\bar{F}(\xi_{\tau_n})}{1 - \tau_n}, \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1 \right) $$

is essentially that of

$$ \left( \begin{matrix} -\frac{1}{\gamma} \frac{\bar{F}(\xi_{\tau_n})}{1 - \tau_n} & \frac{1}{\gamma} \frac{\bar{F}(\xi_{\tau_n})}{1 - \tau_n} \\ 1 & 0 \end{matrix} \right) \sqrt{n(1 - \tau_n)} \left( \begin{matrix} \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1 \\ \frac{\hat{q}_{\alpha_n}}{q_{\alpha_n}} - 1 \end{matrix} \right), $$

and so can be approximated by a centered Gaussian distribution having covariance matrix $M = M_n$ defined by its elements

$$ M_{n,11} = \frac{1}{\gamma^2} \left( \frac{\bar{F}(\xi_{\tau_n})}{1 - \tau_n} \right)^2 \left[ M_{n,11}^{\xi} - 2M_{n,12}^{\xi} + M_{n,22}^{\xi} \right], \\ M_{n,12} = \frac{1}{\gamma} \frac{\bar{F}(\xi_{\tau_n})}{1 - \tau_n} \left[ M_{n,12}^{\xi} - M_{n,11}^{\xi} \right] \text{ and } M_{n,22} = M_{n,11}^{\xi}. $$

It remains to find a precise version of the asymptotic behavior of $\hat{\gamma}_{\tau_n}^{E,BR}$, and for this we start by considering $\hat{\gamma}_{\tau_n}^E$. For this purpose, we set $\psi(x) = 1/(1+x)$ and note the following refinement of the delta-method as a formal power series expansion of $\psi$:

$$ \begin{aligned} & \sqrt{n(1 - \tau_n)}(\hat{\gamma}_{\tau_n}^E - \psi(\bar{F}(\xi_{\tau_n}))/(1 - \tau_n))) \\ &= \sqrt{n(1 - \tau_n)}(\psi(\hat{F}_n(\hat{\xi}_{\tau_n}))/(1 - \tau_n)) - \psi(\bar{F}(\xi_{\tau_n}))/(1 - \tau_n))) \\ &= \sum_{j=1}^{\infty} \frac{\psi^{(j)}(\bar{F}(\xi_{\tau_n}))/(1 - \tau_n))}{j!} [n(1 - \tau_n)]^{(1-j)/2} \left\{ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1 - \tau_n} - \frac{\bar{F}(\xi_{\tau_n})}{1 - \tau_n} \right) \right\}^j. \end{aligned} $$

Evaluating the variance of the right-hand side and its covariance with $\sqrt{n(1 - \tau_n)}(\hat{\xi}_{\tau_n}/\xi_{\tau_n} - 1)$ as formal power series and truncating the resulting series at an order $1/[n(1 - \tau_n)]^J$, for a suitably chosen $J \ge 1$, will lead to an accurate approximation $\mathfrak{M}(J) = \mathfrak{M}_n(J)$ of the covariance matrix of $\sqrt{n(1 - \tau_n)}(\hat{\gamma}_{\tau_n}^E - \gamma, \hat{\xi}_{\tau_n}/\xi_{\tau_n} - 1)$. Of course, first of all, $\mathfrak{M}_{n,22}(J) = M_{n,22}$ since this element is the asymptotic variance of $\sqrt{n(1 - \tau_n)}(\hat{\xi}_{\tau_n}/\xi_{\tau_n} - 1)$ and does not involve the above power series. Then straightforward calculations based on higher order moments and covariances of Gaussian random variables entail

$$ \mathfrak{M}_{n,11}(\infty) = \sum_{j=0}^{\infty} \frac{\mathfrak{m}_{j,11}}{[n(1 - \tau_n)]^j} \text{ with } \mathfrak{m}_{j,11} = (1 + \bar{F}(\xi_{\tau_n}))/(1 - \tau_n)^{-2(j+2)} M_{n,11}^{j+1} \\ \times \left( (2j + 1)!!(2j + 1) - \sum_{i=1}^{j} (2i - 1)!!(2j + 1 - 2i)!! \right) $$
where the double factorial $N!!$ denotes the product of all integers from 1 to $N$ having the same parity as $N$, and

$$ \mathfrak{M}_{n,12}(\infty) = \sum_{j=0}^{\infty} \frac{\mathfrak{m}_{j,12}}{[n(1-\tau_n)]^j} \quad \text{with} \quad \mathfrak{m}_{j,12} = -(1+\overline{F}(\xi_{\tau_n})/(1-\tau_n))^{-2(j+1)} M_{n,11}^j \\ \times M_{n,12}(2j+1)!! $$

The matrix $\mathfrak{M}(J) = \mathfrak{M}_n(J)$, approximating the covariance matrix of $\sqrt{n(1-\tau_n)}(\hat{\gamma}_{\tau_n}^E - \psi(\overline{F}(\xi_{\tau_n})/(1-\tau_n)), \hat{\xi}_{\tau_n}/\xi_{\tau_n} - 1)$, is finally obtained by truncating each of the series defining $\mathfrak{M}_{n,11}(\infty)$ and $\mathfrak{M}_{n,12}(\infty)$ at order $1/[n(1-\tau_n)]^J$. In practice $J=1$ already provides reasonable results, and we adopt this choice in our implementation. We use this in order to approximate the uncertainty about $\hat{\gamma}_{\tau_n}^{\text{E,BR}}$. Straightforward calculations yield

$$ \hat{\gamma}_{\tau_n}^{\text{E,BR}} - \gamma = u_n(\hat{\gamma}_{\tau_n}^E, \hat{\xi}_{\tau_n}/\xi_{\tau_n}) - u_n((1+\overline{F}(\xi_{\tau_n})/(1-\tau_n))^{-1}, 1) + O_P(A((1-\tau_n)^{-1})) $$

where

$$ u_n(x, y) = \left( 1 + (2\tau_n - 1) \left( \frac{1}{x} - 1 \right) \left( 1 - \frac{\mathbb{E}(Y)}{y\xi_{\tau_n}} \right)^{-1} \right)^{-1} $$

Here and throughout the paper, we have neglected the finite-sample variability in $\overline{Y}_n$, and we neglect any term proportional to (or dominated by) $A((1-\tau_n)^{-1})$ for the purpose of approximating the variance of our estimators only. The rationale behind this choice is that, since $A(t) = b\gamma t^\rho$, keeping these terms in this kind of calculation would, because of the bias reduction procedure, entail ultimately having to approximate the correlation of estimators of the second-order parameters $\rho$ and $b$ with estimators of other extreme value parameters, here $\hat{\gamma}_{\tau_n}^E$ and $\hat{\xi}_{\tau_n}$. This is a joint convergence problem which, to the best of our knowledge, remains open, and whose solution deserves a separate in-depth study. Moreover, in many usual cases in extreme value theory (for instance, when $\gamma < 1/2$), one typically has $\gamma < -\rho$, and then $A(t)$ is negligible with respect to $1/q_{1-1/t} \sim (\gamma^{-1}-1)^{-\gamma}/\xi_{1-1/t}$ as $t \to \infty$, which further justifies this modeling choice in our setting. The random vector

$$ \sqrt{n(1-\tau_n)} \left( \hat{\gamma}_{\tau_n}^{\text{E,BR}} - \gamma, \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1 \right) $$

can thus be considered asymptotically Gaussian centered with a covariance matrix approximated by $\mathfrak{M}^{\text{BR}}(J) = \mathfrak{M}_n^{\text{BR}}(J)$ whose elements are

$$ \mathfrak{M}_{n,11}^{\text{BR}}(J) = (\partial_1 u_n((1 + \overline{F}(\xi_{\tau_n})/(1-\tau_n))^{-1}, 1))^2 \mathfrak{M}_{n,11}(J) \\ + (\partial_2 u_n((1 + \overline{F}(\xi_{\tau_n})/(1-\tau_n))^{-1}, 1))^2 \mathfrak{M}_{n,22}(J) \\ + 2\partial_1 u_n((1 + \overline{F}(\xi_{\tau_n})/(1-\tau_n))^{-1}, 1)\partial_2 u_n((1 + \overline{F}(\xi_{\tau_n})/(1-\tau_n))^{-1}, 1)\mathfrak{M}_{n,12}(J), $$

$$ \mathfrak{M}_{n,12}^{\text{BR}}(J) = \partial_1 u_n((1 + \overline{F}(\xi_{\tau_n})/(1 - \tau_n))^{-1}, 1)\mathfrak{M}_{n,12}(J) \\ + \partial_2 u_n((1 + \overline{F}(\xi_{\tau_n})/(1 - \tau_n))^{-1}, 1)\mathfrak{M}_{n,22}(J), $$

$$ \mathfrak{M}_{n,22}^{\text{BR}}(J) = \mathfrak{M}_{n,22}(J). $$

Our final step is to combine all these elements in order to accurately quantify the uncertainty in $\log(\hat{\xi}_{\tau_n}^{\text{BR}}/\xi_{\tau_n})$. Starting with (2.5), and neglecting any term proportional to
$A((1 - \tau_n)^{-1})$ or $A((1 - \tau'_n)^{-1})$ and approximating $\bar{Y}_n$ by $\mathbb{E}(Y)$, we find

$$
\begin{aligned}
\log \frac{\hat{\xi}_{\tau'_n}^{\star, \mathrm{BR}}}{\xi_{\tau'_n}} &\approx \left[ \log \left( \frac{1 - \tau_n}{1 - \tau'_n} \right) + \log \left( \frac{2\tau'_n - 1}{2\tau_n - 1} \right) \right] (\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}} - \gamma) + \log \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} \\
&\quad + \hat{\gamma}_{\tau_n}^{\mathrm{E,BR}} \log \left( 1 - \frac{\mathbb{E}(Y)}{\hat{\xi}_{\tau_n}} \right) - \gamma \log \left( 1 - \frac{\mathbb{E}(Y)}{\xi_{\tau_n}} \right) \\
&\quad - \left[ \hat{\gamma}_{\tau_n}^{\mathrm{E,BR}} \log \left( 1 - \frac{\mathbb{E}(Y)}{\left( \frac{1-\tau'_n}{1-\tau_n} \right)^{-\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}} \hat{\xi}_{\tau_n}} \right) - \gamma \log \left( 1 - \frac{\mathbb{E}(Y)}{\left( \frac{1-\tau'_n}{1-\tau_n} \right)^{-\gamma} \xi_{\tau_n}} \right) \right].
\end{aligned}
$$

This eventually leads to the approximation

$$
\frac{\sqrt{n(1-\tau_n)}}{\log((1-\tau_n)/(1-\tau'_n))} \log \frac{\hat{\xi}_{\tau'_n}^{\star, \mathrm{BR}}}{\xi_{\tau'_n}} \approx \sqrt{n(1-\tau_n)} \left\{ g_n(\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}, \hat{\xi}_{\tau_n}/\xi_{\tau_n}) - g_n(\gamma, 1) \right\}
$$

where

$$
\begin{aligned}
g_n(x, y) = & x \left( 1 + \frac{\log\left(\frac{2\tau'_n - 1}{2\tau_n - 1}\right)}{\log\left(\frac{1 - \tau_n}{1 - \tau'_n}\right)} \right) + \frac{\log(y)}{\log\left(\frac{1 - \tau_n}{1 - \tau'_n}\right)} \\
& + \frac{x}{\log\left(\frac{1 - \tau_n}{1 - \tau'_n}\right)} \log\left(1 - \frac{\mathbb{E}(Y)}{\xi_{\tau_n} y}\right) \\
& - \frac{x}{\log\left(\frac{1 - \tau_n}{1 - \tau'_n}\right)} \log\left(1 - \frac{\mathbb{E}(Y)}{\left(\frac{1 - \tau'_n}{1 - \tau_n}\right)^{-x} \xi_{\tau_n} y}\right).
\end{aligned}
$$

The variance of the asymptotic Gaussian distribution of $(\sqrt{n(1-\tau_n)}/\log((1-\tau_n)/(1-\tau'_n))) \log(\hat{\xi}_{\tau'_n}^{\star,\mathrm{BR}}/\xi_{\tau'_n})$ is then well approximated by

$$
s_n^2(J) = (\partial_1 g_n(\gamma, 1))^2 \mathfrak{M}_{n,11}^{\mathrm{BR}}(J) + 2\partial_1 g_n(\gamma, 1)\partial_2 g_n(\gamma, 1) \mathfrak{M}_{n,12}^{\mathrm{BR}}(J) + (\partial_2 g_n(\gamma, 1))^2 \mathfrak{M}_{n,22}^{\mathrm{BR}}(J).
$$

In order to produce confidence intervals, we estimate $\gamma$ by $\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}$, $\xi_{\tau_n}$ by $\hat{\xi}_{\tau_n}$, $\mathbb{E}(Y)$ by $\bar{Y}_n$,
$\bar{F}(\xi_{\tau_n})$ by $\bar{\hat{F}}_n(\hat{\xi}_{\tau_n})$, $\varphi^{(1)}(\xi_{\tau_n})$ by

$$
\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n}) = \frac{1}{n} \sum_{i=1}^{n} (Y_i - \hat{\xi}_{\tau_n}) \mathbb{1}\{Y_i > \hat{\xi}_{\tau_n}\}.
$$

The estimation of $\varphi^{(2)}(\xi_{\tau_n})$ is more complex, because its naive empirical counterpart
$\hat{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n})$ is unbiased but highly skewed, and therefore tends to vastly underestimate
$\varphi^{(2)}(\xi_{\tau_n})$. A second-order approximation of the underlying distribution function $\bar{F}$ in
a neighborhood of infinity suggests the more accurate estimator

$$
\tilde{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n}) = 2\hat{F}_n(\hat{\xi}_{\tau_n})\hat{\xi}_{\tau_n}^2(\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}})^2 \left( \frac{1}{(1-\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}})(1-2\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}})} + \frac{\bar{b}(\hat{F}_n(\hat{\xi}_{\tau_n}))^{-\bar{\rho}}}{\bar{\rho}} \left\{ \frac{1}{(1-\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}-\bar{\rho})(1-2\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}-\bar{\rho})} - \frac{1}{(1-\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}})(1-2\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}})} \right\} \right).
$$
We deduce from these calculations a corrected asymptotic Gaussian confidence interval for $\xi_{\tau_n'}$ at level $1 - \alpha$ as

$$ \hat{I}_{\tau_n'}^{(2)}(\alpha) = \hat{I}_{\tau_n'}^{(2)}(\alpha; J) = \left[ \hat{\xi}_{\tau_n'}^{\star, \text{BR}} \exp \left( \pm \frac{\log((1 - \tau_n)/(1 - \tau_n'))}{\sqrt{n(1 - \tau_n)}} \sqrt{\hat{s}_n^2(J)} \times z_{1-\alpha/2} \right) \right] $$

where $\hat{s}_n^2(J)$ is obtained from $s_n^2(J)$ by replacing $\gamma$ by $\hat{\gamma}_{\tau_n}^{\text{E,BR}}$, $\xi_{\tau_n}$ by $\hat{\xi}_{\tau_n}$, $\mathbb{E}(Y)$ by $\bar{Y}_n$, $\bar{F}(\xi_{\tau_n})$ by $\bar{F}_n(\hat{\xi}_{\tau_n})$, $\varphi^{(1)}(\xi_{\tau_n})$ by $\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n})$ and $\varphi^{(2)}(\xi_{\tau_n})$ by $\hat{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n})$. The R function `CIextExpect` (with `method="direct"`), available as part of our package `Expectrem`, computes this confidence interval<fnref n="2" />.

We finally state and prove that this corrected confidence interval has asymptotically correct coverage and is always (asymptotically) longer than the naive interval $\hat{I}_{\tau_n'}^{(1)}(\alpha)$.

**Theorem 1.** Assume that $\mathbb{E}(|\min(Y, 0)|^{2+\delta}) < \infty$ for some $\delta > 0$, and condition $C_2(\gamma, \rho, A)$ holds with $0 < \gamma < 1/2$, $\rho < 0$ and $A(t) = b\gamma t^\rho$. Let $\tau_n, \tau_n'$ be two sequences such that $\tau_n, \tau_n' \uparrow 1$ as $n \to \infty$, $n(1 - \tau_n) \to \infty$, $(1 - \tau_n')/(1 - \tau_n) \to 0$ and $\log((1 - \tau_n)/(1 - \tau_n'))/\sqrt{n(1 - \tau_n)} \to 0$ as $n \to \infty$. Assume further that $\sqrt{n(1 - \tau_n)}A((1 - \tau_n)^{-1}) \to 0$ and $\sqrt{n(1 - \tau_n)}/q_{\tau_n} \to \lambda \in \mathbb{R}$, and $\bar{\rho}$ and $\bar{b}$ are consistent estimators of $\rho$ and $b$ such that $(\bar{\rho} - \rho)\log(n) = o_P(1)$. Then:

(i) For any $J \ge 1$, $\hat{s}_n^2(J) \to \gamma^3(1 - \gamma)/(1 - 2\gamma)$ in probability, and

$$ \forall \alpha \in (0, 1), \lim_{n \to \infty} \mathbb{P} \left( \xi_{\tau_n'} \in \hat{I}_{\tau_n'}^{(2)}(\alpha; J) \right) = 1 - \alpha. $$

(ii) One has

$$ \frac{\xi_{\tau_n'}}{\sqrt{n(1 - \tau_n)}} (\text{length}(\hat{I}_{\tau_n'}^{(2)}(\alpha; J)) - \text{length}(\hat{I}_{\tau_n'}^{(1)}(\alpha))) \to \sqrt{\frac{\gamma^3}{(1 - \gamma)(1 - 2\gamma)}} z_{1-\alpha/2} $$

in probability.

# 4 Bias and variance-corrected quantile-based inference

For the purpose of carrying out inference using the quantile-based extrapolated estimator $\tilde{\xi}_{\tau_n'}^{\star, \text{RB}}$, the key is to obtain an accurate representation of the uncertainty in the pair $(\hat{\gamma}_{\tau_n}^{\text{H}}, \tilde{\xi}_{\tau_n})$, where the quantile-based intermediate expectile estimator is $\tilde{\xi}_{\tau_n} = (1/\hat{\gamma}_{\tau_n}^{\text{H}} - 1)^{-\hat{\gamma}_{\tau_n}^{\text{H}}} \hat{q}_{\tau_n}$. Throughout this section we take $\tau_n = 1 - k_n/n$, where $k_n$ is a sequence of integers, $\hat{q}_{\tau_n} = Y_{n-k_n,n}$ is the corresponding intermediate order statistic, and $\hat{\gamma}_{\tau_n}^{\text{H}} = \hat{\gamma}_{1-k_n/n}^{\text{H}}$ is the usual Hill estimator of $\gamma$ calculated upon the top $k_n$ log-spacings. Again, we only give the main steps of our construction, with full technical details reported in Appendix B.2.

The idea is to write

$$ \sqrt{k_n} \log \left( \frac{\tilde{\xi}_{1-k_n/n}}{(\gamma^{-1} - 1)^{-\gamma} q_{1-k_n/n}} \right) = \sqrt{k_n} \{ \phi(\hat{\gamma}_{1-k_n/n}^{\text{H}}) - \phi(\gamma) \} + \sqrt{k_n} \log \left( \frac{\hat{q}_{1-k_n/n}}{q_{1-k_n/n}} \right) $$

<footnote n="2">²In small samples we compute a very slightly different version based on a modification $\tilde{\varphi}_n^{(2)}$ of $\tilde{\varphi}_n^{(2)}$ for added stability; see Section 5.1 for all necessary details.</footnote>
where $\phi(x) = -x \log(x^{-1} - 1)$ and, analogously to Section 3, to construct a formal power series expansion of $\phi$ in a neighborhood of $\gamma$:

$$ \sqrt{k_n} \log \left( \frac{\tilde{\xi}_{1-k_n/n}}{(\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n}} \right) = \sqrt{k_n} \log \left( \frac{\hat{q}_{1-k_n/n}}{q_{1-k_n/n}} \right) + \sum_{j=1}^{\infty} \frac{\phi^{(j)}(\gamma)}{j!} k_n^{\frac{1-j}{2}} \left\{ \sqrt{k_n} (\hat{\gamma}_{1-k_n/n}^{H} - \gamma) \right\}^j . $$

We approximate the variance of $\sqrt{k_n} \log(\tilde{\xi}_{1-k_n/n}/((\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n}))$ and its covariance with $\sqrt{k_n}(\hat{\gamma}_{1-k_n/n}^H - \gamma)$ through a truncated formal power series at an order $1/k_n^J$, for a suitably chosen $J \ge 1$. This provides a finite-sample correction for the asymptotic covariance matrix of $\sqrt{k_n}(\hat{\gamma}_{1-k_n/n}^H - \gamma, \log(\tilde{\xi}_{1-k_n/n}/((\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n})))$, and hence a correction for the asymptotic variance of the extrapolated quantile-based extreme expectile estimator. The resulting matrix $V(J) = V_n(J)$ has the following coefficients: the first variance term $V_{11}(J) = V_{11}(0) = \gamma^2$, obtained from the asymptotic variance of $\sqrt{k_n}(\hat{\gamma}_{1-k_n/n}^H - \gamma)$, does not depend on $J$. The second variance term $V_{22}(J)$, found by evaluating the variance of the right-hand side in the above expansion, is obtained by truncating the series

$$ V_{22}(\infty) = \gamma^2 + \mathrm{Var} \left( \sum_{j=1}^{\infty} \frac{\phi^{(j)}(\gamma)}{j!} k_n^{\frac{1-j}{2}} \gamma^j Z^j \right) $$

where $Z$ is standard Gaussian. [This is justified because, by Theorem 2.4.8 p.52, Lemma 3.2.3 p.71 and Theorem 3.2.5 p.74 in de Haan and Ferreira (2006), $\sqrt{k_n}(\hat{\gamma}_{1-k_n/n}^H - \gamma)$ and $\sqrt{k_n}(\hat{q}_{1-k_n/n}/q_{1-k_n/n} - 1)$ are asymptotically Gaussian and independent with asymptotic variance $\gamma^2$.] After straightforward calculations, gathering powers of $k_n$ together, we find

$$ V_{22}(\infty) = \gamma^2 \left( 1 + (m(\gamma))^2 + \sum_{j=1}^{\infty} \frac{\gamma^{2j}}{k_n^j} \left( (2j+1)!! \sum_{i=1}^{2j+1} \frac{\phi^{(i)}(\gamma)\phi^{(2j+2-i)}(\gamma)}{i!(2j+2-i)!} \right. \right. \\ \left. \left. - \frac{1}{2^{j+1}} \sum_{i=1}^{j} \frac{\phi^{(2i)}(\gamma)\phi^{(2j+2-2i)}(\gamma)}{i!(j+1-i)!} \right) \right) =: \sum_{j=0}^{\infty} \frac{v_{j,22}}{k_n^j} $$

where $m(\gamma) = \phi'(\gamma) = (1 - \gamma)^{-1} - \log(\gamma^{-1} - 1)$. An expression of the formal covariance term $V_{12}(\infty)$ is similarly easily derived: if $Z$ is standard Gaussian,

$$ V_{12}(\infty) = \sum_{j=0}^{\infty} \frac{\phi^{(j)}(\gamma)}{j!} k_n^{\frac{1-j}{2}} \gamma^{j+1} \mathbb{E}(Z^{j+1}) = \sum_{j=0}^{\infty} \frac{\phi^{(2j+1)}(\gamma)}{2^j j! k_n^j} \gamma^{2j+2} =: \sum_{j=0}^{\infty} \frac{v_{j,12}}{k_n^j} $$

In practice we truncate these power series at order $1/k_n^J$, for a given $J \ge 1$, so that a refined approximation to the covariance matrix of $\sqrt{k_n}(\hat{\gamma}_{1-k_n/n}^H - \gamma, \log(\tilde{\xi}_{1-k_n/n}/((\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n})))$ is the symmetric matrix $V(J)$ whose elements are $V_{11}(J), V_{12}(J)$ and $V_{22}(J)$. In our implementation we stop the approximation at order $1/k_n^2$, that is, $J = 2$. Simple expressions of the coefficients $v_{j,12}$ and $v_{j,22}$, for $j = 0, 1, 2$, are provided in Appendix B.2.

We may now proceed with the construction of a confidence interval for $\xi_{\tau_n'}$. We again neglect any term proportional to $A((1 - \tau_n)^{-1})$ or $A((1 - \tau_n')^{-1})$, and we neglect the
finite-sample variability in $\bar{Y}_n$ as an estimator of $\mathbb{E}(Y)$. Starting from Equation (2.6), straightforward calculations lead to the approximation

$$ \log \frac{\tilde{\xi}_{\tau_n'}^{*,\text{BR}}}{\xi_{\tau_n'}} \approx \left( \log \left( \frac{k_n}{n(1-\tau_n')} \right) + \log(2\tau_n' - 1) \right) \left( \hat{\gamma}_{1-k_n/n}^H - \gamma \right) + \log \frac{\tilde{\xi}_{1-k_n/n}}{(\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n}} \\ - \left( \hat{\gamma}_{1-k_n/n}^H \log \left( 1 - \frac{\mathbb{E}(Y)}{\tilde{\xi}_{\tau_n'}^*} \right) - \gamma \log \left( 1 - \frac{\mathbb{E}(Y)}{\left( \frac{n(1-\tau_n')}{k_n} \right)^{-\gamma} (\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n}} \right) \right). $$

This can be reformulated as

$$ \frac{\sqrt{k_n}}{\log(k_n/(n(1-\tau_n')))} \log \frac{\tilde{\xi}_{\tau_n'}^{*,\text{BR}}}{\xi_{\tau_n'}} \approx \sqrt{k_n} \left\{ h_n \left( \hat{\gamma}_{1-k_n/n}^H, \frac{\tilde{\xi}_{1-k_n/n}}{(\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n}} \right) - h_n(\gamma, 1) \right\} $$

where

$$ h_n(x, y) = \left( 1 + \frac{\log(2\tau_n' - 1)}{\log(k_n/(n(1-\tau_n')))} \right) x + \frac{\log(y)}{\log(k_n/(n(1-\tau_n')))} \\ - \frac{x}{\log(k_n/(n(1-\tau_n')))} \log \left( 1 - \frac{\mathbb{E}(Y)}{\left( \frac{n(1-\tau_n')}{k_n} \right)^{-x} (\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n} y} \right). $$

This suggests that the variance of $(\sqrt{k_n}/\log(k_n/(n(1-\tau_n'))))\log(\tilde{\xi}_{\tau_n'}^{*,\text{BR}}/\xi_{\tau_n'})$ is well approximated by

$$ \sigma_n^2(J) = (\partial_1 h_n(\gamma, 1))^2 V_{11}(J) + 2\partial_1 h_n(\gamma, 1)\partial_2 h_n(\gamma, 1)V_{12}(J) + (\partial_2 h_n(\gamma, 1))^2 V_{22}(J). $$

We deduce from these calculations a corrected asymptotic Gaussian confidence interval for $\xi_{\tau_n'}$ at level $1 - \alpha$ as

$$ \tilde{I}_{\tau_n'}^{(2)}(\alpha) = \tilde{I}_{\tau_n'}^{(2)}(\alpha; J) = \left[ \tilde{\xi}_{\tau_n'}^{*,\text{BR}} \exp \left( \pm \frac{\log(k_n/(n(1-\tau_n')))}{\sqrt{k_n}} \sqrt{\tilde{\sigma}_n^2(J)} \times z_{1-\alpha/2} \right) \right] $$

where $\tilde{\sigma}_n^2(J)$ is obtained from $\sigma_n^2(J)$ by replacing $\gamma$ by $\hat{\gamma}_{1-k_n/n}^{H,\text{BR}}$, $q_{1-k_n/n}$ by $\hat{q}_{1-k_n/n}$ and $\mathbb{E}(Y)$ by $\bar{Y}_n$ in the quantities $V_{11}(J) = \gamma^2$, $V_{12}(J)$, $V_{22}(J)$, $\partial_1 h_n(\gamma, 1)$ and $\partial_2 h_n(\gamma, 1)$, thus producing their respective estimators $\hat{V}_{11}(J) = (\hat{\gamma}_{1-k_n/n}^{H,\text{BR}})^2$, $\hat{V}_{12}(J)$, $\hat{V}_{22}(J)$, $\partial_1 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{H,\text{BR}})$, $1$ and $\partial_2 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{H,\text{BR}})$, $1$.

Next, we examine the asymptotic properties of this confidence interval that can be calculated by the R function `CItextExpect` (with `method="indirect"`) available as part of our package `Expectrem`.

**Theorem 2.** Assume that $\mathbb{E}|\min(Y, 0)| < \infty$ and condition $C_2(\gamma, \rho, A)$ holds with $0 < \gamma < 1$, $\rho < 0$ and $A(t) = b\gamma t^\rho$. Let $k_n$, $\tau_n'$ be two sequences such that $k_n \to \infty$, $k_n/n \to 0$, $n(1-\tau_n')/k_n \to 0$ and $\log(k_n/(n(1-\tau_n')))/\sqrt{k_n} \to 0$ as $n \to \infty$. Assume further that $\sqrt{k_n}A(n/k_n) \to 0$ and $\sqrt{k_n}/q_{1-k_n/n} \to \lambda \in \mathbb{R}$, and $\bar{\rho}$ and $\bar{b}$ are consistent estimators of $\rho$ and $b$ such that $(\bar{\rho} - \rho)\log(n) = o_P(1)$. Then:

(i) For any $J \ge 1$, $\tilde{\sigma}_n^2(J) \to \gamma^2$ in probability, and

$$ \forall \alpha \in (0, 1), \lim_{n \to \infty} \mathbb{P} \left( \xi_{\tau_n'} \in \tilde{I}_{\tau_n'}^{(2)}(\alpha; J) \right) = 1 - \alpha. $$
(ii) One has

$$ \frac{\xi_{\tau'_n}}{\sqrt{k_n}}(\text{length}(\tilde{I}_{\tau'_n}^{(2)}(\alpha; J)) - \text{length}(\tilde{I}_{\tau'_n}^{(1)}(\alpha))) \to 2\gamma m(\gamma)z_{1-\alpha/2} $$

in probability, where $m(\gamma) = (1 - \gamma)^{-1} - \log(\gamma^{-1} - 1)$.

The corrected interval $\tilde{I}_{\tau'_n}^{(2)}(\alpha; J)$ is not always longer than the naive Gaussian interval $\tilde{I}_{\tau'_n}^{(1)}(\alpha)$, since $m(\gamma) > 0$ if and only if $\gamma > \gamma_0 \approx 0.218$. Interestingly, even when $\tilde{I}_{\tau'_n}^{(2)}(\alpha, J)$ is shorter than $\tilde{I}_{\tau'_n}^{(1)}(\alpha)$, the former tends to have better coverage than the latter in finite samples; this is due to the fact that coverage of $\tilde{I}_{\tau'_n}^{(1)}(\alpha)$ appears to be higher than the nominal level when $\gamma$ is small. We note moreover that, since $\gamma^3(1 - \gamma)/(1 - 2\gamma) > \gamma^2$ if and only if $\gamma \in ((3 - \sqrt{5})/2, 1/2) \approx (0.382, 0.5)$, the corrected interval $\tilde{I}_{\tau'_n}^{(2)}(\alpha; J)$ will tend to be shorter than $\tilde{I}_{\tau'_n}^{(2)}(\alpha; J)$ when the underlying distribution has a finite third moment. Of course, this rule-of-thumb ignores the variability of the LAWS extreme expectile estimator, which may be considerable in finite samples. We shall check these insights on simulated data in Section 5.

# 5 Finite-sample study

## 5.1 The Expectrem package

The R package Expectrem, freely available at https://github.com/AntoineUC/Expectrem, has been updated to include the methodology developed in this article and fixes have been made to ensure correctness of the implementation. Aside from containing the functions tindexp (to compute either $\hat{\gamma}_{\tau_n}^E$ if argument br=FALSE or $\hat{\gamma}_{\tau_n}^{E,BR}$ if br=TRUE) and extExpect to compute LAWS (if argument method="direct") or quantile-based (if method="indirect") extreme expectile estimators $\hat{\xi}_{\tau_n}^*$ and $\hat{\xi}_{\tau_n}^{*BR}$ (the bias-reduced versions $\hat{\xi}_{\tau_n}^{*BR}$ and $\hat{\xi}_{\tau_n}^{*BR}$ are obtained with br=TRUE), we added two functions in order to make the simulation study and applications introduced in this paper easily reproducible:

* CIExpect: an improved version of extExpect, which returns the extreme expectile estimate, and its associated bias-reduced and variance-corrected confidence interval $\tilde{I}_{\tau'_n}^{(2)}(\alpha; J)$ (for LAWS estimates) and $\tilde{I}_{\tau'_n}^{(2)}(\alpha; J)$ (for quantile-based estimates). The required inputs are the dataset X, the value $k_n = n(1 - \tau_n)$ as k and the target expectile level $\tau'_n$ as tau. The LAWS (method="direct", default) and quantile-based (method="indirect") Weissman estimators can be computed, and the nominal level (ci.level, default 0.95) may be chosen. For comparison purposes, one may also compute the naive, uncorrected confidence intervals $\tilde{I}_{\tau'_n}^{(1)}(\alpha)$ and $\tilde{I}_{\tau'_n}^{(1)}(\alpha)$, available respectively as method="direct_naive" and method="indirect_naive", and univariate versions of the corrected confidence regions of Padoan and Stupfler (2022), available as method="direct_PS" and method="indirect_PS" for LAWS and quantile-based estimators, respectively.

* logspacqqplot: for a dataset X given in input, this function returns an exponential quantile-quantile plot of the $k_n$ (as k) top log-spacings. If argument weighted=TRUE (default), the weighted log-spacings $i \log(X_{n-i+1,n}/X_{n-i,n})$ are computed ($1 \le i \le k$). If weighted=FALSE, they are replaced by $\log(X_{n-i+1,n}/X_{n-k_n,n})$. The straight line with slope $\hat{\gamma}_{1-k_n/n}^{H,BR}$ may be added if add line=TRUE. This is useful for the purpose of diagnosing tail heaviness a posteriori once a choice of $k_n$ has been made.
Note that the expression of $\tilde{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n})$, introduced in Section 3 in the calculation of $\hat{I}_{\tau_n}^{(2)}(\alpha; J)$, can be somewhat unstable in small samples due to the estimation of $\gamma$ and $\rho$. The Cauchy-Schwarz inequality yields $(\varphi^{(1)}(\xi_{\tau_n}))^2 \le \varphi^{(2)}(\xi_{\tau_n}) \le \sqrt{\varphi^{(4)}(\xi_{\tau_n})}$ and motivates the constrained estimator

$$ \tilde{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n}) = \min \left\{ \max \left\{ \tilde{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n}), (\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n}))^2 \right\}, \sqrt{\frac{1}{n} \sum_{i=1}^{n} (Y_i - \hat{\xi}_{\tau_n})^4} \{Y_i > \hat{\xi}_{\tau_n}\} \right\} $$

which we found to work reasonably well in practice in a wide array of situations. This constrained estimator is used in the CItextExpect function (if method="direct") in place of $\tilde{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n})$.

## 5.2 Simulation study

In order to get a comprehensive overview of the practical performance of the suggested confidence intervals, we consider the following test cases for the distribution of Y (see Table 4 for more details):

* A Fréchet distribution with tail index $\gamma > 0$, i.e. $\bar{F}(y) = 1 - \exp(-y^{-1/\gamma})$ for $y > 0$. Here $A(t) = \gamma/(2t)$, and in particular $\rho = -1$. The Fréchet distribution is ubiquitous in extreme value analysis, especially in multivariate problems where it is the standard marginal distribution when estimating multivariate extreme dependence structures, see for instance Chapter 8 of Beirlant et al. (2004).

* A Burr distribution with tail index $\gamma > 0$ and second-order parameter $\rho < 0$, i.e. $\bar{F}(y) = (1 + y^{-\rho/\gamma})^{1/\rho}$ for $y > 0$. Here $A(t) = \gamma t^{\rho}$, so we can make $\rho$ vary in order to generate scenarios with various degrees of difficulty. We consider here $\rho = -5, -1, -0.5$, representing easy, standard and hard scenarios, respectively.

For each of these families of distributions, we consider $\gamma = 0.1, 0.2, 0.3$ and $0.4$ for the LAWS estimator, and $\gamma = 0.1, 0.3, 0.5$ and $0.7$ for the quantile-based estimator: recall that the latter is valid on the larger range of values $\gamma \in (0, 1)$, compared to the former which is restricted to $\gamma \in (0, 1/2)$. In an attempt to evaluate further the influence of finite-sample bias upon the performance of each inference technique, we also consider the following two models:

* A Student distribution with number of degrees of freedom $\nu = 10/3$, i.e. $\gamma = 1/\nu = 0.3$. Here

$$ A(t) = \frac{\nu + 1}{\nu + 2} \left( \frac{\Gamma\left(\frac{\nu+1}{2}\right) \nu^{(\nu-1)/2}}{\sqrt{\nu \pi} \Gamma\left(\frac{\nu}{2}\right)} \right)^{-2/\nu} t^{-2/\nu} \\ \approx 0.635 t^{-0.6}, $$

and in particular $\rho = -2/\nu = -0.6$.

* A Student distribution with number of degrees of freedom $\nu = 10/3$ and location parameter 1, i.e. obtained by adding 1 to a random draw from the previous model. Here again $\gamma = 1/\nu = 0.3$, but by Lemma C.1 in Appendix C.1,

$$ A(t) = -\frac{1}{\nu} \left( \frac{\Gamma\left(\frac{\nu+1}{2}\right) \nu^{(\nu-1)/2}}{\sqrt{\nu \pi} \Gamma\left(\frac{\nu}{2}\right)} \right) t^{-1/\nu} \\ \approx -0.265 t^{-0.3}, $$

and therefore $\rho = -1/\nu = -0.3$.
The rationale for considering these last two models is that while half of the Student realizations will be positive on average, we observed that 80% of the observations generated by the shifted Student model were positive, so that theoretically the techniques considered in this paper will be applicable on a larger range of values of k in the last model compared to the vanilla Student model; this should, however, be at the expense of a (theoretically) larger bias, since the decay of the function A to 0 is slower in the shifted Student model. In each of the considered models, we simulate $N = 5,000$ datasets $\{Y_1, \dots, Y_n\}$ of $n \in \{200, 1,000, 5,000\}$ independent realizations of Y, with distribution function $F = 1 - \bar{F}$. We set $\tau'_n = 1 - 1/n$ and compare the finite-sample performance of the following inference methods for $\xi_{\tau'_n}$:

* The standard confidence intervals $\hat{I}_{\tau'_n}^{(1)}(\alpha)$ and $\tilde{I}_{\tau'_n}^{(1)}(\alpha)$. In order to show the impact of bias reduction, the versions of these intervals not featuring bias correction at all, i.e. where $\hat{\xi}_{\tau'_n}^{*,BR}$, $\tilde{\xi}_{\tau'_n}^{*,BR}$, $\hat{\gamma}_{\tau_n}^{E,BR}$ and $\hat{\gamma}_{\tau_n}^{H,BR}$ are replaced by $\hat{\xi}_{\tau'_n}^{*}$, $\tilde{\xi}_{\tau'_n}^{*}$, $\hat{\gamma}_{\tau_n}^{E}$ and $\hat{\gamma}_{\tau_n}^{H}$, respectively, are also considered, namely

$$ \hat{I}_{\tau'_n}^{(3)}(\alpha) = \left[ \hat{\xi}_{\tau'_n}^{*} \exp \left( \pm \frac{\log((1-\tau_n)/(1-\tau'_n))}{\sqrt{n(1-\tau_n)}} \sqrt{\hat{s}_n^2} \times z_{1-\alpha/2} \right) \right] \\ \text{with } \hat{s}_n^2 = \frac{(\hat{\gamma}_{\tau_n}^E)^3 (1 - \hat{\gamma}_{\tau_n}^E)}{1 - 2\hat{\gamma}_{\tau_n}^E}, \qquad (5.1) $$

and

$$ \tilde{I}_{\tau'_n}^{(3)}(\alpha) = \left[ \tilde{\xi}_{\tau'_n}^{*} \exp \left( \pm \frac{\log((1-\tau_n)/(1-\tau'_n))}{\sqrt{n(1-\tau_n)}} \sqrt{\tilde{\sigma}_n^2} \times z_{1-\alpha/2} \right) \right] \\ \text{with } \tilde{\sigma}_n^2 = (\hat{\gamma}_{\tau_n}^H)^2. \qquad (5.2) $$

* Our proposed confidence intervals $\hat{I}_{\tau'_n}^{(2)}(\alpha; J)$ and $\tilde{I}_{\tau'_n}^{(2)}(\alpha; J)$, with $J = 1$ and $J = 2$, respectively.

* The corrected confidence intervals introduced in Padoan and Stupfler (2022), i.e. univariate versions of their confidence regions $\hat{\mathcal{E}}_{\tau'_n, \alpha}^*$ and $\tilde{\mathcal{E}}_{\tau'_n, \alpha}^*$ for LAWS and quantile-based estimators, respectively. These take the form

$$ \hat{I}_{\tau'_n}^{(4)}(\alpha) = \left[ \hat{\xi}_{\tau'_n}^{*,PS} \exp \left( \pm \frac{\log((1-\tau_n)/(1-\tau'_n))}{\sqrt{n(1-\tau_n)}} \sqrt{\hat{s}_n^2} \times z_{1-\alpha/2} \right) \right] \qquad (5.3) $$

$$ \text{and } \tilde{I}_{\tau'_n}^{(4)}(\alpha) = \left[ \tilde{\xi}_{\tau'_n}^{*,PS} \exp \left( \pm \frac{\log((1-\tau_n)/(1-\tau'_n))}{\sqrt{n(1-\tau_n)}} \sqrt{\tilde{\sigma}_n^2} \times z_{1-\alpha/2} \right) \right] \qquad (5.4) $$

where $\hat{\xi}_{\tau'_n}^{*,PS}$ (resp. $\tilde{\xi}_{\tau'_n}^{*,PS}$) and $\hat{s}_n^2$ (resp. $\tilde{\sigma}_n^2$) are suitably chosen bias-corrected LAWS (resp. quantile-based) estimates of $\xi_{\tau'_n}$ and asymptotic variance estimates for these quantities. We refer to Padoan and Stupfler (2022) for full details.

* Our proposed confidence intervals $\hat{I}_{\tau'_n}^{(2)}(\alpha; J)$ and $\tilde{I}_{\tau'_n}^{(2)}(\alpha; J)$, with $J = 1$ and $J = 2$, respectively.

These are all calculated using the R function `CIExpect` from the `Expectrem` package, described in Section 5.1, at the default $1 - \alpha = 0.95$ confidence level. To assess the finite-sample performance of these inference methods, we provide three types of graphs:
* A graph comparing the empirical variances (across the $N = 5,000$ replicated samples and on the log-scale) of $\hat{\xi}_{\tau_n'}^{\star,\text{BR}}$ (resp. $\hat{\xi}_{\tau_n'}^{\star,\text{BR}}$) with the average of the estimated variances $\hat{s}_n^{2,\text{BR}}$, $\hat{s}_n^2(J)$, $\hat{s}_n^2$ and $\bar{s}_n^2$ (resp. $\tilde{\sigma}_n^{2,\text{BR}}$, $\tilde{\sigma}_n^2(2)$, $\tilde{\sigma}_n^2$ and $\bar{\sigma}_n^2$) giving rise to the considered confidence intervals,

* A graph comparing the empirical coverage probability (across the $N = 5,000$ replicated samples) of each of the confidence intervals $\hat{I}_{\tau_n'}^{(1)}(\alpha)$, $\hat{I}_{\tau_n'}^{(2)}(\alpha; 1)$, $\hat{I}_{\tau_n'}^{(3)}(\alpha)$ and $\hat{I}_{\tau_n'}^{(4)}(\alpha)$ (resp. $\tilde{I}_{\tau_n'}^{(1)}(\alpha)$, $\tilde{I}_{\tau_n'}^{(2)}(\alpha; 2)$, $\tilde{I}_{\tau_n'}^{(3)}(\alpha)$ and $\tilde{I}_{\tau_n'}^{(4)}(\alpha)$) to the nominal level $1-\alpha = 0.95$,

* A graph comparing the median length (across the $N = 5,000$ replicated samples) of each of these confidence intervals.

This results in 96 series of 3 graphs, for the Fréchet and Burr models, and 12 series of 3 graphs, for the (location) Student model. These are deferred to Appendix D.1, as Figures D.1–D.28, for the sake of brevity. We give representative examples of our results in Figures 3 and 4.

As the left panels make clear, the confidence intervals $\hat{I}_{\tau_n'}^{(2)}(\alpha; 1)$ and $\tilde{I}_{\tau_n'}^{(2)}(\alpha; 2)$ are constructed based on variance estimates whose median tends to be much closer to the empirical variance of $\hat{\xi}_{\tau_n'}^{\star,\text{BR}}$ and $\tilde{\xi}_{\tau_n'}^{\star,\text{BR}}$ than those on which the examined competitors are built. That the average value of our variance estimates may be poor should not be surprising since we are working with heavy-tailed data, whose high values are very disruptive to an accurate estimation of the variance. The middle panels indicate that overall there is substantial gain in preferring $\hat{I}_{\tau_n'}^{(2)}(\alpha; 1)$ (resp. $\tilde{I}_{\tau_n'}^{(2)}(\alpha; 2)$) to $\hat{I}_{\tau_n'}^{(1)}(\alpha)$ and $\tilde{I}_{\tau_n'}^{(3)}(\alpha)$ (resp. $\hat{I}_{\tau_n'}^{(1)}(\alpha)$ and $\tilde{I}_{\tau_n'}^{(3)}(\alpha)$). Empirical coverage probabilities of $\hat{I}_{\tau_n'}^{(2)}(\alpha; 1)$ and $\tilde{I}_{\tau_n'}^{(2)}(\alpha; 2)$ are also typically more stable than those of $\hat{I}_{\tau_n'}^{(4)}(\alpha)$ and $\tilde{I}_{\tau_n'}^{(4)}(\alpha)$ as a function of $k$, while generally being as close or even closer to the nominal level (see e.g. the example of the Burr distribution with $\gamma = 0.3$ and $\rho = -1$, with the LAWS approach, or the difficult case of the Burr distribution with $\gamma = 0.1$ and $\rho = -0.5$, with the quantile-based method). The right panels confirm the observations made following Theorems 1 and 2: the LAWS (resp. quantile-based) interval $\hat{I}_{\tau_n'}^{(2)}(\alpha; 1)$ (resp. $\tilde{I}_{\tau_n'}^{(2)}(\alpha; 2)$) is longer than its counterpart $\hat{I}_{\tau_n'}^{(1)}(\alpha)$ (resp. longer than $\tilde{I}_{\tau_n'}^{(1)}(\alpha)$ for values of $\gamma \ge 0.3$), while being kept to a reasonable length. There does not seem to be a rule for comparing these lengths with those of $\hat{I}_{\tau_n'}^{(4)}(\alpha)$ and $\tilde{I}_{\tau_n'}^{(4)}(\alpha)$. Finally, in the Student case, which should be considered difficult because the second-order parameter $\rho$ is then close to 0, it is interesting to note that performance is not worse on the shifted distribution compared to the vanilla Student setting. In fact, at least for small $k$, our variance estimation approach seems to be excellent, while accurately estimating the variance of the (bias-reduced) LAWS estimator for a standard Student sample appears to be very difficult. This is quite surprising, since the second-order auxiliary function $A : t \mapsto -0.265 t^{-0.3}$ of the shifted Student distribution decays more slowly to 0 than the corresponding function $A : t \mapsto 0.635 t^{-0.6}$ for the standard Student distribution. We believe that this is likely due to the constant 0.635, in front of the leading power, being greater in absolute value than -0.265, meaning that while the bias should be theoretically expected to be lower for the Student distribution in very large samples, this is not necessarily the case in smaller samples. This last set of results should be viewed as an illustration of the fact that our bias and variance corrections, and indeed our estimation methods generally, are not location-invariant.
<figure><img src="image_4.png" /><figcaption>Figure 3: Left panels: Comparison of the empirical variance of $(\sqrt{n(1-\tau_n)}/\log((1-\tau_n)/(1-\tau'_n))) \log(\xi_{\tau'_n}^{\chi,\text{BR}}/\xi_{\tau'_n})$ (red curve) with boxplots of the asymptotic variance estimate $\hat{s}_n^2(1)$ and the average values of $\hat{s}_n^{2,\text{BR}}$ (dashed curve), $\hat{s}_n^2(1)$ (green curve), $\hat{s}_n^2$ (dashed-dotted curve) and $\bar{s}_n^2$ (blue curve) as functions of $k=k_n=n(1-\tau_n)$. Middle panels: coverage probabilities of the associated 95% confidence intervals $\hat{I}_{\tau'_n}^{(1)}(\alpha)$, $\hat{I}_{\tau'_n}^{(2)}(\alpha; 1)$, $\hat{I}_{\tau'_n}^{(3)}(\alpha)$ and $\hat{I}_{\tau'_n}^{(4)}(\alpha)$ (respectively in black dashed line, green line, black dotted line and blue line) with the nominal level 0.95 in red. Right panels: Median length of these confidence intervals (same color code as in the middle panels). First row: Burr distribution with $\gamma = 0.1$ and $\rho = -0.5$, second row: Fréchet distribution with $\gamma = 0.2$, third row: Burr distribution with $\gamma = 0.3$ and $\rho = -1$, fourth row: Burr distribution with $\gamma = 0.4$ and $\rho = -5$. The sample size is $n = 1,000$ and the target expectile level is $\tau'_n = 1-1/n = 0.999$.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 4: Left panels: Comparison of the empirical variance of $(\sqrt{n(1-\tau_n)}/\log((1-\tau_n)/(1-\tau'_n))) \log(\xi_{\tau'_n}^{\star,\text{BR}}/\xi_{\tau'_n})$ (red curve) with boxplots of the asymptotic variance estimate $\tilde{\sigma}_n^2(2)$ and the average values of $\tilde{\sigma}_n^{2,\text{BR}}$ (dashed curve), $\tilde{\sigma}_n^2(2)$ (green curve), $\tilde{\sigma}_n^2$ (dashed-dotted curve) and $\bar{\sigma}_n^2$ (blue curve) as functions of $k=k_n=n(1-\tau_n)$. Middle panels: coverage probabilities of the associated 95% confidence intervals $\tilde{I}_{\tau'_n}^{(1)}(\alpha)$, $\tilde{I}_{\tau'_n}^{(2)}(\alpha; 2)$, $\tilde{I}_{\tau'_n}^{(3)}(\alpha)$ and $\tilde{I}_{\tau'_n}^{(4)}(\alpha)$ (respectively in black dashed line, green line, black dotted line and blue line) with the nominal level 0.95 in red. Right panels: Median length of these confidence intervals (same color code as in the middle panels). First row: Burr distribution with $\gamma = 0.1$ and $\rho = -0.5$, second row: Fréchet distribution with $\gamma = 0.3$, third row: Burr distribution with $\gamma = 0.5$ and $\rho = -1$, fourth row: Burr distribution with $\gamma = 0.7$ and $\rho = -5$. The sample size is $n = 1,000$ and the target expectile level is $\tau'_n = 1-1/n = 0.999$.</figcaption></figure>
## 5.3 Real data analyses

We showcase the proposed inference methods on three real data examples: the first two revisit the applications of Daouia et al. (2018) and Girard et al. (2022b) to medical insurance data and stock price data, respectively, and the third one addresses cyber insurance of data breaches through the assessment of, and inference about, their associated extreme risk.

### 5.3.1 Example 1: Society of Actuaries medical insurance claims data

The Society of Actuaries (SOA Group) Medical Insurance Large Claims Database records all the claim amounts exceeding 25,000 USD over the timeframe 1991-92. We only deal here with the $n = 75,789$ claims for 1991. This dataset (available in the R package `Expectrem` as the dataset `soa`) has been analyzed by a number of authors, including Beirlant et al. (2004) who find evidence of a heavy right tail (see also the histogram and rug plot in the top left panel of Figure 5) with tail index $\gamma \approx 1/3$, and Daouia et al. (2018), who estimate extreme quantiles at level $\alpha_n = 1 - 1/100,000$ by making use of extrapolated extreme expectile estimates.

Our objective here is to estimate and infer extreme expectiles themselves in this dataset using our bias-reduced and variance-corrected techniques. The top right panel of Figure 5 compares the point estimates $\hat{\xi}_{\tau_n'}^*$ in (2.1) and $\tilde{\xi}_{\tau_n'}^*$ in (2.2) with their bias-reduced counterparts in (2.5) and (2.6), as functions of the intermediate anchor level $k = k_n = n(1 - \tau_n)$ and at level $\tau_n' = 1 - 1/100,000$. It is readily seen that the bias reduction scheme is highly effective, which here stems from the fact that the estimated value of $\rho$ hovers around $\bar{\rho} = -0.2$ and thus the bias component is expected to be sizeable. We then carry out inference on the basis of the bias-reduced estimates, see the bottom panels of Figure 5, where we compare the straightforward intervals $\hat{I}_{\tau_n'}^{(1)}(\alpha)$ and $\tilde{I}_{\tau_n'}^{(1)}(\alpha)$ to the variance-corrected versions $\hat{I}_{\tau_n'}^{(2)}(\alpha; 1)$ and $\tilde{I}_{\tau_n'}^{(2)}(\alpha; 2)$, respectively. We first observe that the variance-corrected interval is very similar to its naive counterpart when the quantile-based method is used. This is sensible in view of the results of our simulation study, where one sees that the naive and corrected quantile-based extreme expectile confidence intervals are often close when $\gamma \approx 0.3$, which is precisely the case here. When using the LAWS-based confidence intervals, however, the bottom left panel of Figure 5 suggests that the finite-sample variability of the LAWS extreme expectile estimate is more important than what is suggested by the straightforward interval $\hat{I}_{\tau_n'}^{(1)}(\alpha)$. Using the latter would lead to a substantial underestimation of the uncertainty on the target extreme expectile and thus to an incorrect assessment of tail risk carried out from large medical claims with potentially detrimental consequences on the insurance companies involved.

### 5.3.2 Example 2: Commerzbank stock price data

We consider the daily negative log-returns of the Commerzbank stock price on the DAX30 stock exchange between March 6, 2012 and July 28, 2016, resulting in a sample $Y_1, \dots, Y_n$ of size $n = 1,048$ (available in the R package `Expectrem` as the dataset `commerzbank`). This dataset is naturally a time series with substantial serial dependence, so that we cannot apply the bias and variance corrections we develop to the raw data; nevertheless, and as already pointed out by Girard et al. (2022b) for this dataset, the temporal dependence can be handled in a dynamic extreme value estimation setup by, first, filtering the time series with an ARMA(1,1)-GARCH(1,1) model:

$$Y_t = \mu + \phi Y_{t-1} + u_t + \theta u_{t-1},$$
<figure><img src="image_6.png" /><figcaption>Figure 5: Society of Actuaries medical insurance claims data (sample size $n = 75,789$). Top left: Histogram and rug plot of the log-claims $\log Y_i$, $1 \le i \le n$. Top right: Expectile point estimates $\tilde{\xi}_{\tau_n'}^*$ (solid black), $\tilde{\xi}_{\tau_n'}^*$ (dotted black), $\tilde{\xi}_{\tau_n'}^{\star,\text{BR}}$ (solid green) and $\tilde{\xi}_{\tau_n'}^{\star,\text{BR}}$ (dotted green), at level $\tau_n' = 1 - 1/100,000$, as a function of $k = n$ (dashed black) and $k = n(1 - \tau_n)$ (solid green). Bottom left: Bias-reduced LAWS expectile point estimate $\tilde{\xi}_{\tau_n'}^{\star,\text{BR}}$ (solid black) with associated confidence intervals $\hat{I}_{\tau_n'}^{(1)}(\alpha)$ (dashed black) and $\hat{I}_{\tau_n'}^{(2)}(\alpha; 1)$ (dashed green) at the 95% confidence level. Bottom right: Bias-reduced quantile-based expectile point estimate $\tilde{\xi}_{\tau_n'}^{\star,\text{BR}}$ (solid black) with associated confidence intervals $\hat{I}_{\tau_n'}^{(1)}(\alpha)$ (dashed black) and $\hat{I}_{\tau_n'}^{(2)}(\alpha; 2)$ (dashed green) at the 95% confidence level.</figcaption></figure>
where $u_t = \sigma_t \varepsilon_t$ is such that $\sigma_t^2 = c + au_{t-1}^2 + b\sigma_{t-1}^2$, with $(\varepsilon_t)$ being an unobserved independent nondegenerate white noise sequence, i.e. copies of a random variable $\varepsilon$ such that $\mathbb{E}(\varepsilon) = 0$, $\mathbb{E}(\varepsilon^2) = 1$ and $\mathbb{P}(\varepsilon^2 = 1) < 1$, and the constants $\mu, \phi, \theta, a, b$ and $c$ being the model parameters. Denote by $\mathcal{F}_n$ the algebra generated by the ARMA-GARCH process up to time $n$, and recall that the expectile is positive homogeneous and location equivariant, so that the conditional expectile for the next day given the data up to time $n$ can be modeled as

$$ \xi_{\tau}(Y_{n+1}|\mathcal{F}_n) = \mu + \phi Y_n + \sigma_{n+1} \xi_{\tau}(\varepsilon) + \theta u_n. $$

The estimation of expectiles in practice requires an informed choice of the target level. One possible strategy, studied among others by Bellini and Di Bernardino (2017), is to assume that the underlying distribution function $\Phi$ of the innovations $\varepsilon_t$ is standard Gaussian. For example, in this setting, an expectile of level $\tau'_n = 0.995$ of the $\varepsilon_t$ coincides with their quantile at level $\beta_n = \Phi(\xi_{\tau'_n}(\varepsilon)) \approx 0.974$. As such, if the predicted expectiles and quantiles at the levels $\tau'_n = 0.995$ and $0.974$, respectively, are significantly different, then there is evidence that the distribution of the errors is not Gaussian.

We first estimate the ARMA-GARCH model parameters and predict the residuals $\hat{\varepsilon}_i$ and $\hat{u}_i$ using the methodology introduced in Girard et al. (2022b), i.e. by, first, using the function garchFit in the R package fGarch and the option @residuals in order to retrieve the $\hat{u}_i$, and then, by obtaining the predictions $\hat{\varepsilon}_i$ of the innovations by fitting a pure GARCH model directly to the $\hat{u}_i$ and applying garchFit(...)@residuals. Following the theory developed in Girard et al. (2021), we treat the residuals $\hat{\varepsilon}_i$ from the model as independent and identically distributed copies of $\varepsilon$ for the estimation of the tail index $\gamma$ of $\varepsilon$ and its extreme expectile $\xi_{\tau'_n}(\varepsilon)$. The independence assumption was checked using a series of Ljung-Box independence tests on residuals and their squares, with the lowest p-value being 0.39. Evidence that $\varepsilon$ is indeed heavy-tailed is gathered in the top-right panel of Figure 6 using exponential QQ-plots of the log-spacings. The extreme ($\tau'_n = 0.995$) expectile estimates with their associated 90% and 95% confidence intervals (as functions of $k = k_n = n(1 - \tau_n)$) are reported in the middle panels of Figure 6. We eyeball the threshold $k_n$ for stability of the estimates and take in the sequel $k_n = 32$, leading to the tail index estimates $\hat{\gamma}_{k_n}^{\text{H,BR}}(\varepsilon) = 0.329$ and $\hat{\gamma}_{k_n}^{\text{E,BR}}(\varepsilon) = 0.315$. The extreme quantile (of level $\beta_n$, extrapolated and non-extrapolated) estimators are also reported in the middle panels of Figure 6, where

$$ \hat{q}_{\beta_n}(\varepsilon) = \hat{\varepsilon}_{n-|n(1-\beta_n)|,n} \quad \text{and} \\ \hat{q}_{\beta_n}^{\star, \text{BR}}(\varepsilon) = \left( \frac{k_n}{n(1-\beta_n)} \right)^{\hat{\gamma}_{k_n}^{\text{H,BR}}(\varepsilon)} \hat{\varepsilon}_{n-k_n, n} \left( 1 + \frac{1}{\bar{\rho}} \left\{ \left( \frac{k_n}{n(1-\beta_n)} \right)^{\bar{\rho}} - 1 \right\} \times \bar{b} \hat{\gamma}_{k_n}^{\text{H,BR}}(\varepsilon) \left( \frac{n}{k_n} \right)^{\bar{\rho}} \right) $$

are the empirical and bias-reduced extrapolated Weissman quantile estimates of $q_{\beta_n}(\varepsilon)$. The resulting quantile and expectile estimates are reported, with confidence intervals, in Table 1 for $k_n = 32$.

Dynamic predictions of extreme expectiles of $Y_{n+1}$ given $Y_n$ are then obtained as

$$ \hat{\xi}_{\tau'_n}^{\star, \text{BR}}(Y_{n+1}|\mathcal{F}_n) = \hat{\mu} + \hat{\phi} Y_n + \hat{\sigma}_{n+1} \hat{\xi}_{\tau'_n}^{\star, \text{BR}}(\varepsilon) + \hat{\theta} \hat{u}_n \\ \text{and} \\ \tilde{\xi}_{\tau'_n}^{\star, \text{BR}}(Y_{n+1}|\mathcal{F}_n) = \hat{\mu} + \hat{\phi} Y_n + \hat{\sigma}_{n+1} \tilde{\xi}_{\tau'_n}^{\star, \text{BR}}(\varepsilon) + \hat{\theta} \hat{u}_n. $$

In order to visualize these predictions, the key point is that, if $\hat{\mu}, \hat{\phi}, \hat{\theta}$ and the past $(Y_t)_{t \le n-1}$ of the process are considered as fixed, then $\hat{u}_n = Y_n - \hat{\phi}Y_{n-1} - \hat{\theta}\hat{u}_{n-1} - \hat{\mu}$
<figure><img src="image_7.png" /><figcaption>Figure 6: Commerzbank stock price data (sample size *n* = 1,048). Top left: Daily log-returns (black curve) and ARMA-GARCH log-volatility estimates (red bold curve, smoothed using the R function *smooth.spline* with smoothing parameter *λ* = 5 × 10⁻⁷). Top right: Exponential QQ-plot of the log-spacings log(ε̂<sub>n-i+1,n</sub>/ε̂<sub>n-k_n,n</sub>), 1 ≤ *i* ≤ *k*_n = 32. The straight line has slope γ̂<sub>k_n</sub><sup>H,BR</sup>(ε) = 0.329. Middle panels: Expectile estimates ξ<sub>τ'_n</sub><sup>*,BR</sup>(ε) and ξ̃<sub>τ'_n</sub><sup>*,BR</sup>(ε) at τ'_n = 0.995 (solid blue and green curves) with their associated confidence intervals I<sub>τ'_n</sub><sup>(2)</sup>(α; 1) and I<sub>τ'_n</sub><sup>(2)</sup>(α; 2) (dashed blue and green curves) at the 90% (left) and 95% (right) confidence levels, along with quantile estimates q̂<sub>β_n</sub>(ε) and q<sub>β_n</sub><sup>*,BR</sup>(ε) at level *β*_n ≈ 0.974 (dashed-dotted and solid red curves). Bottom panels: Dynamic estimates and their associated 90% (left) and 95% (right) confidence intervals, with the same color code.</figcaption></figure>
<table><caption>Table 1: 90% (first column) and 95% (second column) confidence intervals for the extreme expectile $\xi_{\tau'_n}$ of the innovations $\varepsilon_t$ at $\tau'_n = 0.995$, along with their associated LAWS and quantile-based point estimates (third column), and empirical and Weissman quantile estimates (fourth and fifth columns) calculated at the level $\beta_n \approx 0.974$.</caption><thead><tr><th>Interval</th><th>90% CI</th><th>95% CI</th><th>ξ<sub>τ'<sub>n</sub></sub>(ε) estimate</th><th>q̂<sub>β<sub>n</sub></sub>(ε)</th><th>q̂<sup>*</sup><sub>β<sub>n</sub></sub>(ε)</th></tr></thead><tbody><tr><td>Î<sub>τ'<sub>n</sub></sub><sup>(1)</sup>(α)</td><td>[2.435, 3.136]</td><td>[2.377, 3.213]</td><td rowspan="2">ξ<sub>τ'<sub>n</sub></sub><sup>BR</sup> ≈ 2.763</td><td rowspan="2">2.017</td><td rowspan="2">1.999</td></tr><tr><td>Î<sub>τ'<sub>n</sub></sub><sup>(2)</sup>(α)</td><td>[2.091, 3.652]</td><td>[1.982, 3.852]</td></tr><tr><td>Î<sub>τ'<sub>n</sub></sub><sup>(1)</sup>(α)</td><td>[2.337, 3.305]</td><td>[2.261, 3.416]</td><td rowspan="2">ξ<sub>τ'<sub>n</sub></sub><sup>BR</sup> ≈ 2.779</td><td rowspan="2"></td><td rowspan="2"></td></tr><tr><td>Î<sub>τ'<sub>n</sub></sub><sup>(2)</sup>(α)</td><td>[2.129, 3.627]</td><td>[2.023, 3.817]</td></tr></tbody></table>

(resp. $\hat{\sigma}_{n+1}^2$) is an affine function of $Y_n$ (resp. is quadratic in $\hat{u}_n$ and hence in $Y_n$). The dynamic extreme expectile estimates can thus be seen as functions of the observation $Y_n = y_n$. These functions are represented in the bottom panels of Figure 6, together with confidence intervals deduced from those of $\xi_{\tau'_n}(\varepsilon)$ by ignoring the uncertainty in $\hat{\mu}$, $\hat{\phi}$, $\hat{\theta}$, $\hat{\alpha}$, $\hat{\beta}$ and $\hat{\gamma}$: this is justified theoretically by the fact that these estimators converge at the rate $1/\sqrt{n}$ and thus much faster than those obtained via our extreme value procedures. They are compared with dynamic bias-reduced Weissman quantile estimates of level $\beta_n$, calculated as

$$ \hat{q}_{\beta_n}^{*, \text{BR}}(Y_{n+1} | \mathcal{F}_n) = \hat{\mu} + \hat{\phi} Y_n + \hat{\sigma}_{n+1} \hat{q}_{\beta_n}^{*, \text{BR}}(\varepsilon) + \hat{\theta} \hat{u}_n. $$

It is readily seen that the quantile estimates do not belong to the 90% confidence intervals of the expectile estimates, and belong to the LAWS-based 95%–confidence interval but not to the quantile-based 95%–confidence interval. This constitutes further evidence that the residuals of the ARMA-GARCH model are not Gaussian and therefore that the standard selection rules of expectile levels on light-tailed data do not apply to stock price data in general. A financial consequence of this analysis is that a risk assessment based on expectiles using the guidelines provided by Bellini and Di Bernardino (2017) in the Gaussian case would be more conservative than if it were based on quantiles.

### 5.3.3 Example 3: Privacy Rights Clearinghouse cyber risk data

Cyber risks, especially the subclass of personal data breaches from firms/organisations, defined by Article 4(12) of the European Union General Data Protection Regulation as breaches of security leading to mass identity fraud, are a rapidly increasing threat to individuals, companies, public services, and governments. Insurance against data breaches constitutes the bulk of coverage of losses related to cyber events. The financial loss due to a single breached piece of private data is estimated to be 213 USD on average, which translates into a global cost estimated at hundreds of billions of USD per year, see Wheatley et al. (2016). Early work on data breach risks from an extreme value perspective includes quantifying the heavy-tailed nature of breach sizes in Maillart and Sornette (2010) and Wheatley et al. (2016). Very recently, Farkas et al. (2021) have computed loss quantiles for a cyber portfolio by combining Generalized Pareto modeling and regression trees. We further investigate such questions here using expectiles. As described below in our exploratory analysis, tail index estimates are found to be mostly larger than 0.5 in the different studied classes of breached data, thus leading to more conservative (i.e. prudent) risk measurements against catastrophic breach events than if quantiles were used, according to Bellini et al. (2014).
We use the PRC (Privacy Rights Clearinghouse) database<fnref n="3" /> which is the most comprehensive open scientific dataset for breaches occurring in businesses, educational institutions, government and military, healthcare institutions, and other media and nonprofit organisations. The PRC database gathers cyber events from different sources, which introduces heterogeneity among the reported events. In order to reduce the heterogeneity effect when inferring the tail risk, we follow the setup of Farkas et al. (2021) and let the sources of information be grouped into four types: (1) US federal government agencies; (2) US state government agencies; (3) Nonprofit organizations; and (4) Media. To strive towards checking the fundamental assumption of identically distributed data, we stratify the breaches occurring in healthcare, medical providers and medical insurance services (represented by the MED acronym) and educational institutions (represented by the EDU acronym) along the source of information and breach event type. This differentiation of the cyber events (after eliminating duplicates and NA values) occurring between January 1st, 2005 and October 25th, 2019, results in 6 MED (respectively, 3 EDU) classes of sizes $n > 100$; we exclude small clusters of less than 100 observations. As summarized in Table 2, the classes we obtain are determined by the following five types of breaches: HACK (hacking by an outside party or infection by malware), INSD (similar to HACK but performed by an insider), PHYS (physical breach, including paper documents that are lost, discarded or stolen), PORT (portable device lost, discarded or stolen) and DISC (unintended disclosure not involving hacking, intentional breach or physical loss), and by Federal, State and Nonprofit sources of information (the Media source appears only in very small clusters of size $n < 50$).

A more global and automated clustering is adopted by Farkas et al. (2021) based on their GPD regression tree analysis of the severity distribution, but nothing guarantees that our model assumptions (e.g. stationarity) are satisfied when using the leaves of their trees. It should also be noted that the PRC database only reports breach sizes, that is, the number of records affected by each breach event. The related financial loss is not provided directly by the database, but can be quantified approximately from a claims-driven loss formula, calibrated by Farkas et al. (2021) to the specific PRC database. We refer to Section 5.1 in Farkas et al. (2021) for a detailed discussion and justification of the rationale behind their calibration. The resulting losses provide a rough approximation of the real claim data. Lack of this data still constitutes a bottleneck to precise risk assessment and insurance pricing, but the analysis presented here might already help “risk thinkers” to better grasp these new cyber risks from the extreme value perspective.

Numerical results obtained for the nine studied classes are displayed in Table 2. The tuning parameter $k = n(1 - \tau_n)$ appearing in both $\hat{\gamma}_{\tau_n}^{H,BR}$ and $\tilde{\xi}_{1-1/n}^{*,BR}$ is chosen following the path stability algorithm developed by El Methni and Stupfler (2017) which consists in computing the standard deviations of the estimators over a “moving window” of successive values of $k$, and then by selecting the first value of $k$ where the standard deviation is minimal and sufficiently low for each estimator. Graphs of $\hat{\gamma}_{1-k/n}^{H,BR}$ versus $k$, along with asymptotic Gaussian 95% confidence intervals and the value of $k$ selected by the path stability procedure as well as the final pointwise estimate, and their analogs for $\tilde{\xi}_{1-1/n}^{*,BR}$ can be found in Appendix D.2.

The results gathered in Table 2 show a substantial difference in the bias-corrected Hill estimates $\hat{\gamma}_{\tau_n}^{H,BR} \in [0.40, 0.78]$ of the tail index across all clusters. We also observe 7 extremely heavy-tailed cyber clusters with tail index estimates exceeding 0.5 (as indicated in green), in which the estimated expectile risk measure at level $\tau'_n = 1 - 1/n$ exceeds the maximum historical loss (as indicated in purple). Our asymptotic corrected Gaussian

<footnote n="3">³See https://privacyrights.org/data-breaches. Data on file with the authors.</footnote>
confidence intervals point towards upper risk bounds whose value may even exceed three times the maximum loss. An accurate quantification of this uncertainty is of course crucial, especially in this challenging context with small datasets. We provide a further point of comparison by estimating a quantile-based Value-at-Risk using the bias-corrected extrapolated Weissman estimator $\hat{q}_{1-1/n}^{*,BR}$. It is precisely on the 7 extremely heavy-tailed cyber risk classes where the difference between extreme expectiles and quantiles is strongest. In particular, the estimated expectile loss measures can be as much as three times higher than their quantile counterparts, and the quantile estimates do not even exceed the maximum historical loss for five of these classes (as indicated in orange). It is here more prudent to measure tail risk based on expectiles rather than on quantiles, which solely depend on the frequency of tail losses and not on their severity.

<table><caption>Table 2: Privacy Rights Clearinghouse cyber risk data. First column: Sector, level and type of breach. Second column: Sample size *n* of each class. Third column: Bias-reduced Hill estimate of $\gamma$ along with the selected *k* and asymptotic Gaussian 95% confidence interval. Fourth column: Quantile-based bias-reduced extreme expectile estimator $\tilde{\xi}_{1-1/n}^{*,BR}$ along with the selected *k* and bias-reduced and variance-corrected asymptotic Gaussian 95% confidence interval $\tilde{I}_{1-1/n}^{(2)}(\alpha; 2)$. Fifth column: Maximum observation. Sixth column: Weissman extreme quantile estimator $\hat{q}_{1-1/n}^{*,BR}$ obtained by making use of the bias-reduced Hill estimator, along with the selected *k* and asymptotic Gaussian 95% confidence interval.</caption><thead><tr><th>Sector, level, type of breach</th><th>n</th><th>γ̂<sup>H,BR</sup><sub>τ<sub>n</sub></sub> (k̂<sub>opt</sub>) [95% CI]</th><th>ξ̂<sup>,BR</sup><sub>1-1/n</sub> (k̂<sub>opt</sub>) (million USD) [95% CI]</th><th>max<sub>1≤i≤n</sub> Y<sub>i</sub> (million USD)</th><th>q̂<sup>,BR</sup><sub>1-1/n</sub> (k̂<sub>opt</sub>) (million USD) [95% CI]</th></tr></thead><tbody><tr><td>MED, federal, PHYS</td><td>1140</td><td>0.654 (320) [0.582, 0.726]</td><td>156.04 (322) [79.91, 304.70]</td><td>95.12</td><td>84.17 (316) [55.49, 127.68]</td></tr><tr><td>MED, federal, DISC</td><td>766</td><td>0.582 (299) [0.516, 0.648]</td><td>61.54 (323) [36.26, 104.44]</td><td>57.07</td><td>33.51 (314) [23.11, 48.58]</td></tr><tr><td>MED, federal, HACK</td><td>545</td><td>0.782 (263) [0.688, 0.877]</td><td>545.01 (265) [181.26, 1638.73]</td><td>463.34</td><td>157.36 (266) [92.96, 266.36]</td></tr><tr><td>MED, state, HACK</td><td>178</td><td>0.497 (67) [0.378, 0.616]</td><td>17.49 (69) [9.21, 33.22]</td><td>60.26</td><td>13.07 (42) [7.21, 23.69]</td></tr><tr><td>MED, nonprofit, PORT</td><td>223</td><td>0.714 (54) [0.523, 0.904]</td><td>169.99 (42) [37.90, 762.35]</td><td>60.26</td><td>99.86 (56) [47.12, 211.64]</td></tr><tr><td>MED, nonprofit, INSD</td><td>106</td><td>0.752 (42) [0.525, 0.980]</td><td>38.03 (43) [4.40, 328.57]</td><td>16.57</td><td>12.18 (40) [5.13, 28.92]</td></tr><tr><td>EDU, nonprofit, HACK</td><td>178</td><td>0.402 (36) [0.270, 0.533]</td><td>29.17 (43) [16.57, 51.34]</td><td>32.63</td><td>28.54 (68) [17.49, 46.56]</td></tr><tr><td>EDU, nonprofit, DISC</td><td>174</td><td>0.529 (27) [0.329, 0.728]</td><td>21.48 (27) [7.10, 64.94]</td><td>21.13</td><td>15.98 (28) [8.53, 29.91]</td></tr><tr><td>EDU, nonprofit, PORT</td><td>103</td><td>0.578 (46) [0.411, 0.745]</td><td>25.03 (46) [8.80, 71.17]</td><td>10.52</td><td>11.78 (58) [6.36, 2181]</td></tr></tbody></table>

## Acknowledgments

The authors acknowledge an anonymous reviewer for helpful comments that led to an improved version of this article. This research was supported by the French National Research Agency under the grants ANR-19-CE40-0013 (ExtremReg project), ANR-17-EURE-0010 (EUR CHESS) and ANR-11-LABX-0020-01 (Centre Henri Lebesgue). A. Daouia and G. Stupfler acknowledge financial support from the TSE-HEC ACPR Chair and from an AXA Research Fund Award on “Mitigating risk in the wake of the COVID-19 pandemic”.
References

Angus, J. (1993). Asymptotic theory for bootstrapping the extremes. *Communications in Statistics - Theory and Methods*, 22(1):15-30.

Artzner, P., Delbaen, F., Eber, J.-M., and Heath, D. (1999). Coherent measures of risk. *Mathematical Finance*, 9(3):203-228.

Athreya, K. B. (1987). Bootstrap of the mean in the infinite variance case. *Annals of Statistics*, 15(2):724-731.

Beirlant, J., Goegebeur, Y., Segers, J., and Teugels, J. (2004). *Statistics of Extremes: Theory and Applications*. Wiley, Chichester, United Kingdom.

Bellini, F. and Di Bernardino, E. (2017). Risk management with expectiles. *The European Journal of Finance*, 23(6):487-506.

Bellini, F., Klar, B., Müller, A., and Gianin, E. R. (2014). Generalized quantiles as risk measures. *Insurance: Mathematics and Economics*, 54:41-48.

Caeiro, F., Gomes, M. I., and Pestana, D. (2005). Direct reduction of bias of the classical Hill estimator. *Revstat*, 3(2):113-136.

Cornea-Madeira, A. and Davidson, R. (2015). A parametric bootstrap for heavy-tailed distributions. *Econometric Theory*, 31(3):449-470.

Daouia, A., Girard, S., and Stupfler, G. (2018). Estimation of tail risk based on extreme expectiles. *Journal of the Royal Statistical Society: Series B*, 80(2):263-292.

Daouia, A., Girard, S., and Stupfler, G. (2019). Extreme M-quantiles as risk measures: from $L^1$ to $L^p$ optimization. *Bernoulli*, 25(1):264-309.

Daouia, A., Girard, S., and Stupfler, G. (2020). Tail expectile process and risk assessment. *Bernoulli*, 26(1):531-556.

de Haan, L. and Ferreira, A. (2006). *Extreme Value Theory: An Introduction*. Springer-Verlag, New York, New York.

El Methni, J. and Stupfler, G. (2017). Extreme versions of Wang risk measures and their estimation for heavy-tailed distributions. *Statistica Sinica*, 27(2):907-930.

Embrechts, P., Klüppelberg, C., and Mikosch, T. (1997). *Modelling Extremal Events for Insurance and Finance*. Springer, Berlin-Heidelberg.

Farkas, S., Lopez, O., and Thomas, M. (2021). Cyber claim analysis using Generalized Pareto regression trees with applications to insurance. *Insurance: Mathematics and Economics*, 98:92-105.

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

Girard, S., Stupfler, G., and Usseglio-Carleve, A. (2022a). Nonparametric extreme conditional expectile estimation. *Scandinavian Journal of Statistics*, 49(1):78-115.
Girard, S., Stupfler, G., and Usseglio-Carleve, A. (2022b). On automatic bias reduction for extreme expectile estimation. *Statistics & Computing*, 32(4):64.

Gneiting, T. (2011). Making and evaluating point forecasts. *Journal of the American Statistical Association*, 106(494):746–762.

Gomes, M. I., Brilhante, M. F., and Pestana, D. (2016). New reduced-bias estimators of a positive extreme value index. *Communications in Statistics-Simulation and Computation*, 45(3):833–862.

Guillou, A. (2000). Bootstrap confidence intervals for the Pareto index. *Communications in Statistics - Theory and Methods*, 29(1):211–226.

Hall, P. (1990). Asymptotic properties of the bootstrap for heavy-tailed distributions. *Annals of Probability*, 18(3):1342–1360.

Hall, P. and Jing, B.-Y. (1998). Comparison of bootstrap and asymptotic approximations to the distribution of a heavy-tailed mean. *Statistica Sinica*, 8(3):887–906.

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

Holzmann, H. and Klar, B. (2016). Expectile asymptotics. *Electronic Journal of Statistics*, 10(2):2355–2371.

Jones, M. C. (1994). Expectiles and M-quantiles are quantiles. *Statistics & Probability Letters*, 20:149–153.

Knight, K. (1989). On the bootstrap of the sample mean in the infinite variance case. *Annals of Statistics*, 17(3):1168–1175.

Krätschmer, V. and Zähle, H. (2017). Statistical inference for expectile-based risk measures. *Scandinavian Journal of Statistics*, 44(2):425–454.

Maillart, T. and Sornette, D. (2010). Heavy-tailed distribution of cyber-risks. *European Physical Journal B*, 75(3):357–364.

Newey, W. K. and Powell, J. L. (1987). Asymmetric least squares estimation and testing. *Econometrica*, 55(4):819–847.

Padoan, S. A. and Stupfler, G. (2022). Joint inference on extreme expectiles for multivariate heavy-tailed distributions. *Bernoulli*, 28(2):1021–1048.

Stupfler, G. and Usseglio-Carleve, A. (2023). Composite bias-reduced $L^p$-quantile-based estimators of extreme quantiles and expectiles. *Canadian Journal of Statistics*, 51(2):704–742.

Wang, J. and Genton, M. G. (2006). The multivariate skew-slash distribution. *Journal of Statistical Planning and Inference*, 136(1):209–220.

Wheatley, S., Maillart, T., and Sornette, D. (2016). The extreme risk of personal data breaches and the erosion of privacy. *European Physical Journal B*, 89(7):1–12.

Ziegel, J. F. (2016). Coherence and elicitability. *Mathematical Finance*, 26(4):901–918.
# Appendix to the paper “Bias-reduced and variance-corrected asymptotic Gaussian inference about extreme expectiles”

Abdelaati Daouia, Gilles Stupfler & Antoine Usseglio-Carleve

This appendix is organized as follows. Section A describes the estimators of the second-order parameters used in our implementation. Section B gives further details as to how our corrections are calculated. Section C contains all necessary proofs and Section D provides extra finite-sample results about our simulation study and real data analyses.

## A Estimation of the second-order parameters b and ρ

Under the working assumption that $A(t) = b\gamma t^{\rho}$, the estimators $\bar{b} = \bar{b}_n$ and $\bar{\rho} = \bar{\rho}_n$ of the b and ρ parameters are calculated as follows: for a given n-sample $(Y_1, \dots, Y_n)$, let $Y_{1,n} \le Y_{2,n} \le \dots \le Y_{n,n}$ be the associated order statistics and set

$$M_{\kappa_n}^{(j)} = \frac{1}{\kappa_n} \sum_{i=1}^{\kappa_n} (\log Y_{n-i+1,n} - \log Y_{n-\kappa_n,n})^j, \text{ for } j = 1, 2, 3.$$

For $j = 1$, this is just the Hill estimator. These quantities are the basic building blocks for the quantity $T_{\kappa_n}^{(\tau)}$ defined as

$$T_{\kappa_n}^{(\tau)} = \begin{cases} \frac{(M_{\kappa_n}^{(1)})^{\tau} - (M_{\kappa_n}^{(2)}/2)^{\tau/2}}{(M_{\kappa_n}^{(2)}/2)^{\tau/2} - (M_{\kappa_n}^{(3)}/6)^{\tau/3}} & \text{if } \tau > 0, \\ \frac{\log(M_{\kappa_n}^{(1)}) - \frac{1}{2}\log(M_{\kappa_n}^{(2)}/2)}{\frac{1}{2}\log(M_{\kappa_n}^{(2)}/2) - \frac{1}{3}\log(M_{\kappa_n}^{(3)}/6)} & \text{if } \tau = 0. \end{cases}$$

The estimator of ρ that we consider is a simple function of $T_{\kappa_n}^{(\tau)}$:

$$\hat{\rho}_{\kappa_n}^{(\tau)} = - \left| \frac{3(T_{\kappa_n}^{(\tau)} - 1)}{T_{\kappa_n}^{(\tau)} - 3} \right|. \qquad (A.1)$$

This estimator is implemented in the R function mop, available in the Expectrem package<fnref n="4" />. In this package, with $\kappa_n = \lfloor n^{0.999} \rfloor$, and a choice of $\tau \in \{0, 1\}$ is made based on a stability criterion for $\kappa \mapsto \hat{\rho}_{\kappa}^{(\tau)}$ for large $\kappa$ (see Section 3.2 in Gomes et al. (2016) for more details). According to Proposition 2.1 in Caeiro et al. (2005), these choices ensure, if ρ is large enough (a calculation analogous to that in Remark 2.2 in Caeiro et al. (2005) provides roughly $\rho > -249.75$, which will cover all practical applications), that $(\hat{\rho}_{\kappa_n}^{(\tau)} - \rho) \log(n) = o_P(1)$ as required in our asymptotic results. An estimator of b is then

$$\hat{b}_{\kappa_n} = \left(\frac{\kappa_n}{n}\right)^{\bar{\rho}} \frac{\left(\frac{1}{\kappa_n} \sum_{i=1}^{\kappa_n} \left(\frac{i}{\kappa_n}\right)^{-\bar{\rho}}\right) \left(\frac{1}{\kappa_n} \sum_{i=1}^{\kappa_n} U_i\right) - \left(\frac{1}{\kappa_n} \sum_{i=1}^{\kappa_n} \left(\frac{i}{\kappa_n}\right)^{-\bar{\rho}} U_i\right)}{\left(\frac{1}{\kappa_n} \sum_{i=1}^{\kappa_n} \left(\frac{i}{\kappa_n}\right)^{-\bar{\rho}}\right) \left(\frac{1}{\kappa_n} \sum_{i=1}^{\kappa_n} \left(\frac{i}{\kappa_n}\right)^{-\bar{\rho}} U_i\right) - \left(\frac{1}{\kappa_n} \sum_{i=1}^{\kappa_n} \left(\frac{i}{\kappa_n}\right)^{-2\bar{\rho}} U_i\right)}, \quad (A.2)$$

<footnote n="4">⁴Taken from the package evt0, unfortunately not available anymore from CRAN at the time this paper was written.</footnote>
where $\bar{\rho} = \hat{\rho}_{\kappa_n}^{(\tau)}$ and the $U_i = i \log(Y_{n-i+1,n}/Y_{n-i,n})$ are the weighted log-spacings. This estimator is also available from the R function mop. The aforementioned choice of $\kappa_n$ ensures that $\bar{b} = \hat{b}_{\kappa_n}$ is consistent, see Proposition 2.2 in Caeiro et al. (2005). Our results will therefore always feature the conditions that $(\bar{\rho} - \rho) \log(n) = o_P(1)$ and $\bar{b} - b = o_P(1)$.

# B Detailed calculations for the construction of LAWS and quantile-based extreme expectile confidence intervals

## B.1 LAWS confidence interval construction

Our argument is motivated by an inspection of the existing proof of the joint asymptotic normality of an intermediate empirical expectile and the corresponding expectile-based tail index estimator $\hat{\gamma}_{\tau_n}^E$ (see Theorem 2 in Stupfler and Usseglio-Carleve (2023) and Theorem 2 in Girard et al. (2022b)). The idea of the proof is not new; however, carefully identifying the crucial steps of this proof and the key approximations that are made will be instrumental in our subsequent construction of refined confidence intervals in which the novelty of the present work resides. We assume throughout this section that $\mathbb{E}|\min(Y, 0)|^2 < \infty$ and $\gamma < 1/2$.

### B.1.1 Preliminary steps

The basic idea of the proof behind Theorem 2 in Stupfler and Usseglio-Carleve (2023) is that for any $z_1, z_2 \in \mathbb{R}$,

$$ \left\{ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1 - \tau_n} - (\gamma^{-1} - 1) \right) \le z_1, \quad \sqrt{n(1 - \tau_n)} \left( \frac{\hat{\xi}_{\tau_n}}{\hat{\xi}_{\tau_n}} - 1 \right) \le z_2 \right\} \\ = \left\{ \hat{\xi}_{\tau_n} \ge \hat{q}_{1-(1-\tau_n)(\gamma^{-1}-1+z_1/\sqrt{n(1-\tau_n)}),} \quad \sqrt{n(1 - \tau_n)} \left( \frac{\hat{\xi}_{\tau_n}}{\hat{\xi}_{\tau_n}} - 1 \right) \le z_2 \right\}. \quad (B.1) $$

Evaluating the joint uncertainty in $(\hat{\gamma}_{\tau_n}^E, \hat{\xi}_{\tau_n})$ is therefore equivalent to evaluating the joint uncertainty in $(\hat{\xi}_{\tau_n}, \hat{q}_{\alpha_n})$, where $\alpha_n$ is a sequence satisfying $1 - \alpha_n = (\gamma^{-1} - 1)(1 - \tau_n)(1 + o(1))$. Pick then $u, v \in \mathbb{R}$ and set

$$ A_n(u, v) = \left\{ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{\xi}_{\tau_n}}{\hat{\xi}_{\tau_n}} - 1 \right) \le u \right\} \cap \left\{ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{q}_{\alpha_n}}{q_{\alpha_n}} - 1 \right) \le v \right\}. $$

Setting $y_n = y_n(u) = \xi_{\tau_n}(1 + u/\sqrt{n(1 - \tau_n)})$ and $z_n = z_n(v) = q_{\alpha_n}(1 + v/\sqrt{n(1 - \tau_n)})$, it follows that

$$ A_n(u, v) = \left\{ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{E}_n(y_n)}{\bar{E}(y_n)} - 1 \right) \le \sqrt{n(1 - \tau_n)} \left( \frac{\bar{E}(\xi_{\tau_n})}{\bar{E}(y_n)} - 1 \right) \right\} \\ \cap \left\{ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{F}_n(z_n)}{\bar{F}(z_n)} - 1 \right) \le \sqrt{n(1 - \tau_n)} \left( \frac{\bar{F}(q_{\alpha_n})}{\bar{F}(z_n)} - 1 \right) \right\}. \quad (B.2) $$

There are two crucial points of note here. On the one hand, it is straightforward to obtain, via a Lyapunov central limit theorem, that

$$ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{E}_n(y_n)}{\bar{E}(y_n)} - 1, \frac{\hat{F}_n(z_n)}{\bar{F}(z_n)} - 1 \right) \xrightarrow{d} \mathcal{N} \left( (0, 0), \frac{\gamma}{1 - \gamma} \Sigma \right), \quad (B.3) $$
where the 2 × 2 symmetric matrix Σ has entries $\Sigma_{11} = 2(1 - \gamma)/(1 - 2\gamma)$, $\Sigma_{12} = \Sigma_{22} = 1$. On the other hand,

$$ \sqrt{n(1 - \tau_n)} \left( \frac{\overline{E}(\xi_{\tau_n})}{\overline{E}(y_n)} - 1 \right) = \frac{u}{\gamma}(1 + o(1)) \quad (B.4) $$

$$ \text{and} \quad \sqrt{n(1 - \tau_n)} \left( \frac{\overline{F}(q_{\alpha_n})}{\overline{F}(z_n)} - 1 \right) = \frac{v}{\gamma}(1 + o(1)). $$

This is shown using the second-order regular variation property and Lemma A.3(iv) in Stupfler and Usseglio-Carleve (2023). One now readily concludes that

$$ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1, \frac{\hat{q}_{\alpha_n}}{q_{\alpha_n}} - 1 \right) \xrightarrow{d} \mathcal{N} \left( (0,0), \frac{\gamma^3}{1 - \gamma} \Sigma \right). \quad (B.5) $$

The final step is to recall Equation (B.1) and to use the fact that

$$ \{\hat{\xi}_{\tau_n} \ge \hat{q}_{\alpha_n}\} = \left\{ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1 \right) \ge \sqrt{n(1 - \tau_n)} \left( \frac{\hat{q}_{\alpha_n}}{q_{\alpha_n}} - 1 \right) \frac{q_{\alpha_n}}{\xi_{\tau_n}} + \sqrt{n(1 - \tau_n)} \left( \frac{q_{\alpha_n}}{\xi_{\tau_n}} - 1 \right) \right\}. $$

Now recall that, from Proposition 1 in Daouia et al. (2020) and the second-order regular variation assumption,

$$ \frac{q_{\alpha_n}}{\xi_{\tau_n}} = \frac{q_{\alpha_n}}{q_{\tau_n}} \times \frac{q_{\tau_n}}{\xi_{\tau_n}} = \left( 1 - \frac{\gamma^2}{1 - \gamma} \times \frac{z_1}{\sqrt{n(1 - \tau_n)}} (1 + o(1)) \right) \times \left( 1 - \frac{\gamma(\gamma^{-1} - 1)^\gamma}{q_{\tau_n}} (\mathbb{E}(Y) + o(1)) - \left( \frac{(\gamma^{-1} - 1)^{-\rho}}{1 - \gamma - \rho} + o(1) \right) A((1 - \tau_n)^{-1}) \right) $$

when $\alpha_n = \alpha_n(z_1) = 1 - (1 - \tau_n)(\gamma^{-1} - 1 + z_1/\sqrt{n(1 - \tau_n)})$. Neglecting bias terms and isolating $z_1$, one now readily finds, via (B.1), that

$$ \left\{ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1 - \tau_n} - (\gamma^{-1} - 1) \right) \le z_1, \sqrt{n(1 - \tau_n)} \left( \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1 \right) \le z_2 \right\} \\ = \left\{ -\frac{1 - \gamma}{\gamma^2} \times \sqrt{n(1 - \tau_n)} \left( \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1 \right) + \frac{1 - \gamma}{\gamma^2} \times \sqrt{n(1 - \tau_n)} \left( \frac{\hat{q}_{\alpha_n}}{q_{\alpha_n}} - 1 \right) + o_P(1) \le z_1, \right. \\ \left. \sqrt{n(1 - \tau_n)} \left( \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1 \right) \le z_2 \right\}. $$

This means that the asymptotic covariance matrix of

$$ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1 - \tau_n} - (\gamma^{-1} - 1), \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1 \right) $$

is essentially that of

$$ \left( \begin{array}{cc} -\frac{1 - \gamma}{\gamma^2} & \frac{1 - \gamma}{\gamma^2} \\ 1 & 0 \end{array} \right) \times \sqrt{n(1 - \tau_n)} \left( \frac{\frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1}{\frac{\hat{q}_{1-(1-\tau_n)(\gamma^{-1}-1)}}{q_{1-(1-\tau_n)(\gamma^{-1}-1)}} - 1} \right). \quad (B.6) $$
The asymptotic covariance structure of $\sqrt{n(1-\tau_n)}(\hat{\gamma}_{\tau_n}^E - \gamma, \hat{\xi}_{\tau_n}/\xi_{\tau_n} - 1)$ is then obtained by applying the delta-method with the function $x \mapsto (h(x), y) = (1/(1+x), y)$: one finds

$$ \sqrt{n(1-\tau_n)} \left( \hat{\gamma}_{\tau_n}^E - \gamma, \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1 \right) \xrightarrow{d} \mathcal{N} \left( \left( \frac{\gamma(\gamma^{-1}-1)^{1-\rho}}{1-\gamma-\rho} \lambda_1 + \gamma^2(\gamma^{-1}-1)^{\gamma+1} \mathbb{E}(Y)\lambda_2, 0 \right), \frac{\gamma^3}{1-2\gamma} \begin{pmatrix} 1-\gamma & 1 \\ 1 & 2 \end{pmatrix} \right) \quad (B.7) $$

where $\lambda_1 = \lim_{n \to \infty} \sqrt{n(1-\tau_n)} A((1-\tau_n)^{-1})$ and $\lambda_2 = \lim_{n \to \infty} \sqrt{n(1-\tau_n)}/q_{\tau_n}$, assumed henceforth to be finite. Recall finally that, if $(1-\tau'_n)/(1-\tau_n) \to 0$,

$$ \left( \frac{1 - \tau'_n}{1 - \tau_n} \right)^{\gamma} \frac{\xi_{\tau'_n}}{\xi_{\tau_n}} = 1 + O(1/\sqrt{n(1 - \tau_n)}), $$

see Proposition 1 in Daouia et al. (2020) and the proof of Theorem 4.3.8 p.139 in de Haan and Ferreira (2006). Using the identity

$$ \log \left( \frac{\hat{\xi}_{\tau_n}^*}{\xi_{\tau_n}^*} \right) = (\hat{\gamma}_{\tau_n}^E - \gamma) \log \left( \frac{1 - \tau_n}{1 - \tau_n^*} \right) + \log \left( \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} \right) - \log \left( \left( \frac{1 - \tau_n^*}{1 - \tau_n} \right)^{\gamma} \frac{\xi_{\tau_n^*}}{\xi_{\tau_n}} \right) $$

provides (2.3). The fact that $\bar{b} - b = o_P(1)$ and $(\bar{\rho} - \rho) \log(n) = o_P(1)$ then yields the same asymptotic normality result for $\hat{\xi}_{\tau_n}^{*,BR}$, only with the asymptotic mean being 0.

## B.1.2 Calculations

The idea is to construct versions of steps (B.3), (B.4), (B.5), (B.6) and (B.7) that are more accurate in finite samples. We do so based on the sequence $\alpha_n = 1 - \bar{F}(\xi_{\tau_n})$, which is asymptotically equivalent (but not equal) to the sequence $1 - (1 - \tau_n)(\gamma^{-1} - 1)$ corresponding to the choice $z_1 = 0$ in (B.1). This choice is motivated by the fact that for the construction of asymptotic Gaussian confidence regions, one is typically interested in an accurate calculation of the probability of events such as (B.1) for small values of ($z_1, z_2$), because this is the region where most of the mass of the Gaussian limiting distribution concentrates. Taking $\alpha_n = 1 - \bar{F}(\xi_{\tau_n})$, rather than $\alpha_n = 1 - (1 - \tau_n)(\gamma^{-1} - 1)$, makes it possible to more accurately quantify the uncertainty in $\hat{F}_n(\hat{\xi}_{\tau_n})$ by comparing it directly to its population counterpart $\bar{F}(\xi_{\tau_n})$.

We first examine what room for improvement there is when using the asymptotic approximation (B.3). We do so with $y_n$ replaced by $\xi_{\tau_n}$ and $z_n$ replaced by $q_{\alpha_n}$, i.e. we take $u = v = 0$ in $y_n = y_n(u)$ and $z_n = z_n(v)$, which is once again appropriate for our purpose of correctly evaluating the probability of $A_n(u, v)$ for small $u$ and $v$. It is straightforward to prove that the covariance matrix of

$$ \sqrt{n(1-\tau_n)} \left( \frac{\hat{\varphi}_n^{(1)}(\xi_{\tau_n})}{\varphi^{(1)}(\xi_{\tau_n})} - 1, \frac{\hat{F}_n(q_{\alpha_n})}{\bar{F}(q_{\alpha_n})} - 1 \right) $$

is exactly the $2 \times 2$ symmetric matrix $M^{\varphi}$, where $M^{\varphi} = M_n^{\varphi}$ has components:

$$ M_{n,11}^{\varphi} = (1 - \tau_n) \left( \frac{\varphi^{(2)}(\xi_{\tau_n})}{[\varphi^{(1)}(\xi_{\tau_n})]^2} - 1 \right), \\ M_{n,12}^{\varphi} = (1 - \tau_n) \left( \frac{\varphi^{(1)}(\xi_{\tau_n} \vee q_{\alpha_n}) + (\xi_{\tau_n} \vee q_{\alpha_n} - \xi_{\tau_n}) \bar{F}(\xi_{\tau_n} \vee q_{\alpha_n})}{\varphi^{(1)}(\xi_{\tau_n})(1 - \alpha_n)} - 1 \right) \\ = (1 - \tau_n) \frac{\alpha_n}{1 - \alpha_n} \quad \text{and} \quad M_{n,22}^{\varphi} = (1 - \tau_n) \frac{\alpha_n}{1 - \alpha_n} = M_{12}^{\varphi}. $$
In particular

$$ M_n^\varphi \to \begin{pmatrix} \frac{2\gamma}{1-2\gamma} & \frac{\gamma}{1-\gamma} \\ \frac{\gamma}{1-\gamma} & \frac{\gamma}{1-\gamma} \end{pmatrix} \quad \text{as } n \to \infty. $$

Recall that $\hat{E}_n(y) = \hat{\varphi}_n^{(1)}(y)/(2\hat{\varphi}_n^{(1)}(y) + y - \bar{Y}_n)$ and therefore

$$ \begin{aligned} \frac{\hat{E}_n(\xi_{\tau_n})}{\bar{E}(\xi_{\tau_n})} - 1 &= \frac{\hat{\varphi}_n^{(1)}(\xi_{\tau_n})(\xi_{\tau_n} - \mathbb{E}(Y)) - \varphi^{(1)}(\xi_{\tau_n})(\xi_{\tau_n} - \bar{Y}_n)}{\varphi^{(1)}(\xi_{\tau_n})(2\hat{\varphi}_n^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \bar{Y}_n)} \\ &= \frac{(\hat{\varphi}_n^{(1)}(\xi_{\tau_n}) - \varphi^{(1)}(\xi_{\tau_n}))(\xi_{\tau_n} - \mathbb{E}(Y)) + \varphi^{(1)}(\xi_{\tau_n})(\bar{Y}_n - \mathbb{E}(Y))}{\varphi^{(1)}(\xi_{\tau_n})(2\hat{\varphi}_n^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \bar{Y}_n)} \\ &\approx \frac{\xi_{\tau_n} - \mathbb{E}(Y)}{2\varphi^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \mathbb{E}(Y)} \left( \frac{\hat{\varphi}_n^{(1)}(\xi_{\tau_n})}{\varphi^{(1)}(\xi_{\tau_n})} - 1 \right) + \frac{1}{2\varphi^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \mathbb{E}(Y)} (\bar{Y}_n - \mathbb{E}(Y)). \end{aligned} $$

Since $\bar{Y}_n - \mathbb{E}(Y)$ converges to 0 at the rate $1/\sqrt{n}$ by the central limit theorem, this identity suggests that it is in fact reasonable to write

$$ \frac{\hat{E}_n(\xi_{\tau_n})}{\bar{E}(\xi_{\tau_n})} - 1 \approx \frac{\xi_{\tau_n} - \mathbb{E}(Y)}{2\varphi^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \mathbb{E}(Y)} \left( \frac{\hat{\varphi}_n^{(1)}(\xi_{\tau_n})}{\varphi^{(1)}(\xi_{\tau_n})} - 1 \right). $$

As a consequence, the joint distribution of

$$ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{E}_n(\xi_{\tau_n})}{\bar{E}(\xi_{\tau_n})} - 1, \frac{\hat{F}_n(q_{\alpha_n})}{\bar{F}(q_{\alpha_n})} - 1 \right) $$

has a covariance matrix that can be more accurately approximated by $M^E = M_n^E$ having components:

$$ M_{n,11}^E = \frac{(\xi_{\tau_n} - \mathbb{E}(Y))^2 M_{11}^\varphi}{(2\varphi^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \mathbb{E}(Y))^2}, \quad M_{n,12}^E = \frac{(\xi_{\tau_n} - \mathbb{E}(Y))M_{12}^\varphi}{2\varphi^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \mathbb{E}(Y)}, \quad M_{n,22}^E = M_{22}^\varphi. $$

[Clearly $M_n^E$ is elementwise asymptotically equivalent to $M_n^\varphi$.] We now focus on obtaining a more accurate version of (B.5). Recall from (B.2) that

$$ \begin{aligned} \mathbb{P}(A_n(u, v)) &= \mathbb{P} \left( \left\{ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{E}_n(y_n)}{\bar{E}(y_n)} - 1 \right) \le \sqrt{n(1 - \tau_n)} \left( \frac{\bar{E}(\xi_{\tau_n})}{\bar{E}(y_n)} - 1 \right) \right\} \right. \\ & \qquad \left. \cap \left\{ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{F}_n(z_n)}{\bar{F}(z_n)} - 1 \right) \le \sqrt{n(1 - \tau_n)} \left( \frac{\bar{F}(q_{\alpha_n})}{\bar{F}(z_n)} - 1 \right) \right\} \right) \\ & \approx \mathbb{P} \left( \left\{ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{E}_n(\xi_{\tau_n})}{\bar{E}(\xi_{\tau_n})} - 1 \right) \le \sqrt{n(1 - \tau_n)} \left( \frac{\bar{E}(\xi_{\tau_n})}{\bar{E}(y_n)} - 1 \right) \right\} \right. \\ & \qquad \left. \cap \left\{ \sqrt{n(1 - \tau_n)} \left( \frac{\hat{F}_n(q_{\alpha_n})}{\bar{F}(q_{\alpha_n})} - 1 \right) \le \sqrt{n(1 - \tau_n)} \left( \frac{\bar{F}(q_{\alpha_n})}{\bar{F}(z_n)} - 1 \right) \right\} \right) \end{aligned} $$
where $y_n = y_n(u) = \xi_{\tau_n}(1 + u/\sqrt{n(1 - \tau_n)})$ and $z_n = z_n(v) = q_{\alpha_n}(1 + v/\sqrt{n(1 - \tau_n)})$; the approximation is motivated by the fact that $\overline{E}_n(y_n)/\overline{E}(y_n) - 1$ and $\overline{E}_n(\xi_{\tau_n})/\overline{E}(\xi_{\tau_n}) - 1$ have the same asymptotic behavior, and similarly for $\overline{F}_n(z_n)/\overline{F}(z_n) - 1$ and $\overline{F}_n(q_{\alpha_n})/\overline{F}(q_{\alpha_n}) - 1$. Note also that the function $\overline{E}$ is absolutely continuous, because $y \mapsto \varphi^{(1)}(y) = \int_y^\infty \overline{F}(t) dt$ is Lipschitz continuous and the denominator $y \mapsto \mathbb{E}(|Y - y|) = 2\varphi^{(1)}(y) + y - \mathbb{E}(Y)$ defines a Lipschitz continuous function that is bounded away from zero. The function $\overline{E}$ has Lebesgue derivative

$$ \overline{E}'(y) = - \frac{\varphi^{(1)}(y) + \overline{F}(y)(y - \mathbb{E}(Y))}{(2\varphi^{(1)}(y) + y - \mathbb{E}(Y))^2}. $$

This suggests the use of a Taylor expansion to obtain the following more precise version of (B.4):

$$ \frac{\overline{E}(\xi_{\tau_n})}{\overline{E}(y_n)} - 1 = u\xi_{\tau_n} \frac{1}{\sqrt{n(1-\tau_n)}} \frac{\varphi^{(1)}(\xi_{\tau_n}) + \overline{F}(\xi_{\tau_n})(\xi_{\tau_n} - \mathbb{E}(Y))}{\varphi^{(1)}(\xi_{\tau_n})(2\varphi^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \mathbb{E}(Y))} (1 + o(1)). $$

Therefore, using (B.2), the distribution of $(\xi_{\tau_n}/\xi_{\tau_n} - 1, \hat{q}_{\alpha_n}/q_{\alpha_n} - 1)$ is approximated by that of

$$ \left( \frac{1}{\xi_{\tau_n}} \frac{\varphi^{(1)}(\xi_{\tau_n})(2\varphi^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \mathbb{E}(Y))}{\varphi^{(1)}(\xi_{\tau_n}) + \overline{F}(\xi_{\tau_n})(\xi_{\tau_n} - \mathbb{E}(Y))} \left( \frac{\overline{E}_n(\xi_{\tau_n})}{\overline{E}(\xi_{\tau_n})} - 1 \right), \gamma \left( \frac{\overline{F}_n(q_{\alpha_n})}{\overline{F}(q_{\alpha_n})} - 1 \right) \right). $$

An approximation of the covariance matrix of $\sqrt{n(1 - \tau_n)}(\xi_{\tau_n}/\xi_{\tau_n} - 1, \hat{q}_{\alpha_n}/q_{\alpha_n} - 1)$ is then given by $\mathbf{M}^\xi = \mathbf{M}_n^\xi$, whose elements are

$$ M_{n,11}^{\xi} = \left( \frac{\varphi^{(1)}(\xi_{\tau_n})(2\varphi^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \mathbb{E}(Y))}{\xi_{\tau_n}(\varphi^{(1)}(\xi_{\tau_n}) + \overline{F}(\xi_{\tau_n})(\xi_{\tau_n} - \mathbb{E}(Y)))} \right)^2 M_{n,11}^E, $$

$$ M_{n,12}^{\xi} = \gamma \left( \frac{\varphi^{(1)}(\xi_{\tau_n})(2\varphi^{(1)}(\xi_{\tau_n}) + \xi_{\tau_n} - \mathbb{E}(Y))}{\xi_{\tau_n}(\varphi^{(1)}(\xi_{\tau_n}) + \overline{F}(\xi_{\tau_n})(\xi_{\tau_n} - \mathbb{E}(Y)))} \right) M_{n,12}^E, \text{ and } M_{n,22}^{\xi} = \gamma^2 M_{n,22}^E. $$

Note that indeed

$$ \mathbf{M}_n^\xi \to \frac{\gamma^3}{1-\gamma} \mathbf{\Sigma} = \begin{pmatrix} \frac{2\gamma^3}{1-2\gamma} & \frac{\gamma^3}{1-\gamma} \\ \frac{\gamma^3}{1-\gamma} & \frac{\gamma^3}{1-\gamma} \end{pmatrix} \quad \text{as } n \to \infty, $$

as a straightforward consequence of Lemma 3 in Stupfler and Usseglio-Carleve (2023) involving burdensome calculations, which use the asymptotic equivalents

$$ \varphi^{(1)}(y) = \frac{\gamma}{1-\gamma} y \overline{F}(y) (1+o(1)) \quad \text{and} \quad \varphi^{(2)}(y) = \frac{2\gamma^2}{(1-\gamma)(1-2\gamma)} y^2 \overline{F}(y) (1+o(1)), $$

combined with the convergence $\overline{F}(\xi_\tau)/(1-\tau) \to \gamma^{-1}-1$ as $\tau \uparrow 1$. This means that $\mathbf{M}^\xi$ provides a corrected version of the asymptotic covariance matrix of $\sqrt{n(1-\tau_n)}(\xi_{\tau_n}/\xi_{\tau_n} - 1, \hat{q}_{\alpha_n}/q_{\alpha_n} - 1)$ that is consistent with the asymptotics in (B.5).

We turn to finding an improvement of (B.6). Note that

$$ \left\{ \sqrt{n(1-\tau_n)} \left( \frac{\overline{F}_n(\hat{\xi}_{\tau_n})}{1-\tau_n} - \frac{\overline{F}(\xi_{\tau_n})}{1-\tau_n} \right) \le z \right\} = \left\{ \hat{\xi}_{\tau_n} \ge \hat{q}_{1-(1-\tau_n)(\overline{F}(\xi_{\tau_n})/(1-\tau_n)+z/\sqrt{n(1-\tau_n)})} \right\}. $$
The intermediate quantile level $\beta_n = \beta_n(z) = 1 - (1 - \tau_n)(\overline{F}(\xi_{\tau_n})/(1 - \tau_n) + z/\sqrt{n(1 - \tau_n)})$ in the right-hand side is such that $1 - \beta_n \sim \overline{F}(\xi_{\tau_n}) \sim (\gamma^{-1} - 1)(1 - \tau_n)$. Then

$$
\begin{aligned}
& \mathbb{P}\left(\left\{\sqrt{n(1-\tau_n)}\left(\frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1-\tau_n} - \frac{\overline{F}(\xi_{\tau_n})}{1-\tau_n}\right) \le z\right\} \cap \left\{\sqrt{n(1-\tau_n)}\left(\frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1\right) \le z'\right\}\right) \\
& = \mathbb{P}\left(\left\{\sqrt{n(1-\tau_n)}\left(\frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1\right) \ge \sqrt{n(1-\tau_n)}\left(\frac{\hat{q}_{\beta_n}}{q_{\beta_n}} - 1\right) \frac{q_{\beta_n}}{\xi_{\tau_n}} + \sqrt{n(1-\tau_n)}\left(\frac{q_{\beta_n}}{\xi_{\tau_n}} - 1\right)\right\} \cap \left\{\sqrt{n(1-\tau_n)}\left(\frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1\right) \le z'\right\}\right) \\
& \approx \mathbb{P}\left(\left\{\sqrt{n(1-\tau_n)}\left(\frac{\hat{q}_{\alpha_n}}{q_{\alpha_n}} - 1\right) \frac{q_{\alpha_n}}{\xi_{\tau_n}} - \sqrt{n(1-\tau_n)}\left(\frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1\right) \le -\sqrt{n(1-\tau_n)}\left(\frac{q_{\beta_n}}{\xi_{\tau_n}} - 1\right)\right\} \cap \left\{\sqrt{n(1-\tau_n)}\left(\frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1\right) \le z'\right\}\right)
\end{aligned}
$$

where again $\alpha_n = \beta_n(0) = F(\xi_{\tau_n})$, because $(\hat{q}_{\beta_n}/q_{\beta_n} - 1)q_{\beta_n}/\xi_{\tau_n}$ and $(\hat{q}_{\alpha_n}/q_{\alpha_n} - 1)q_{\alpha_n}/\xi_{\tau_n}$ have the same asymptotic behavior. From Proposition 1 in Daouia et al. (2020) and the second-order regular variation assumption,

$$ \sqrt{n(1-\tau_n)} \left( \frac{q_{\beta_n}}{\xi_{\tau_n}} - 1 \right) = \sqrt{n(1-\tau_n)} \left( \frac{q_{\beta_n}}{q_{\tau_n}} \times \frac{q_{\tau_n}}{\xi_{\tau_n}} - 1 \right) \approx -\gamma \frac{1-\tau_n}{\overline{F}(\xi_{\tau_n})} z $$

up to bias terms, which will be taken care of in the bias reduction procedure (and therefore can be neglected at this stage). Conclude that the asymptotic distribution of

$$ \sqrt{n(1-\tau_n)} \left( \frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1-\tau_n} - \frac{\overline{F}(\xi_{\tau_n})}{1-\tau_n}, \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1 \right) $$

can be more accurately approximated by that of

$$ \left( \begin{matrix} -\frac{1}{\gamma} \frac{\overline{F}(\xi_{\tau_n})}{1-\tau_n} & \frac{1}{\gamma} \frac{\overline{F}(\xi_{\tau_n})}{1-\tau_n} \\ 1 & 0 \end{matrix} \right) \sqrt{n(1-\tau_n)} \left( \begin{matrix} \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1 \\ \frac{\hat{q}_{\alpha_n}}{q_{\alpha_n}} - 1 \end{matrix} \right) $$

As a consequence, this asymptotic distribution can be approximated by a centered Gaussian distribution having covariance matrix $\mathbf{M} = \mathbf{M}_n$ defined by its elements

$$
\begin{aligned}
M_{n,11} &= \frac{1}{\gamma^2} \left( \frac{\overline{F}(\xi_{\tau_n})}{1-\tau_n} \right)^2 [M_{n,11}^{\xi} - 2M_{n,12}^{\xi} + M_{n,22}^{\xi}], \\
M_{n,12} &= \frac{1}{\gamma} \frac{\overline{F}(\xi_{\tau_n})}{1-\tau_n} [M_{n,12}^{\xi} - M_{n,11}^{\xi}] \text{ and } M_{n,22} = M_{n,11}^{\xi}.
\end{aligned}
$$

This is the required refinement of (B.6), satisfying

$$ \mathbf{M}_n = \begin{pmatrix} M_{n,11} & M_{n,12} \\ M_{n,12} & M_{n,22} \end{pmatrix} \rightarrow \begin{pmatrix} \frac{1-\gamma}{\gamma(1-2\gamma)} & -\frac{\gamma}{1-2\gamma} \\ -\frac{\gamma}{1-2\gamma} & \frac{2\gamma^3}{1-2\gamma} \end{pmatrix} \quad \text{as } n \to \infty. $$
It remains to find a more precise version of (B.7) tailored to $\hat{\gamma}_{\tau_n}^{\text{E,BR}}$. This convergence relies on applying the delta-method to $\hat{F}_n(\hat{\xi}_{\tau_n})/(1 - \tau_n)$ with the transformation $\psi: x \mapsto 1/(1+x)$. An even more accurate approximation is possible by, first, noting that $\hat{\gamma}_{\tau_n}^{\text{E}} = \psi(\hat{F}_n(\hat{\xi}_{\tau_n})/(1-\tau_n))$ is an (asymptotically) unbiased estimator of $\psi(F(\xi_{\tau_n})/(1-\tau_n))$, and then by writing a formal power series expansion of $\psi$ in order to assess the behavior of $\hat{\gamma}_{\tau_n}^{\text{E}} - \psi(F(\xi_{\tau_n})/(1-\tau_n))$. This is done as follows:

$$ 
\begin{aligned}
& \sqrt{n(1-\tau_n)}(\hat{\gamma}_{\tau_n}^{\text{E}} - \psi(\overline{F}(\xi_{\tau_n}))/(1-\tau_n))) \\
& = \sqrt{n(1-\tau_n)}(\psi(\overline{F}_n(\hat{\xi}_{\tau_n}))/(1-\tau_n)) - \psi(\overline{F}(\xi_{\tau_n}))/(1-\tau_n))) \\
& = \sum_{j=1}^{\infty} \frac{\psi^{(j)}(\overline{F}(\xi_{\tau_n}))/(1-\tau_n))}{j!} [n(1-\tau_n)]^{(1-j)/2} \left\{ \sqrt{n(1-\tau_n)} \left( \frac{\overline{F}_n(\hat{\xi}_{\tau_n})}{1-\tau_n} - \frac{\overline{F}(\xi_{\tau_n})}{1-\tau_n} \right) \right\}^j.
\end{aligned}
$$

Evaluating the variance of the right-hand side and its covariance with $\sqrt{n(1-\tau_n)}(\hat{\xi}_{\tau_n}/\xi_{\tau_n} - 1)$ and truncating the resulting (formal) power series at an order $1/[n(1-\tau_n)]^J$, for a suitably chosen $J \ge 1$, will lead to a matrix $\mathfrak{M}(J) = \mathfrak{M}_n(J)$ whose first two rows and columns will represent a correction for the asymptotic variance matrix in (B.7), with higher values of $J$ intuitively linked to a more accurate approximation. Of course, first of all, $\mathfrak{M}_{n,22}(J) = M_{n,22}$ since this asymptotic variance of $\sqrt{n(1-\tau_n)}(\hat{\xi}_{\tau_n}/\xi_{\tau_n} - 1)$ does not involve the above power series. Then, recalling that

$$ \sqrt{n(1-\tau_n)} \left( \frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1-\tau_n} - \frac{\overline{F}(\xi_{\tau_n})}{1-\tau_n} \right) \approx \sqrt{M_{n,11}} Z $$

where $Z$ is standard Gaussian, $\mathfrak{M}_{n,11}(J)$ is obtained by calculating first

$$ \mathfrak{M}_{n,11}(\infty) = \text{Var} \left( \sum_{j=1}^{\infty} \frac{\psi^{(j)}(\overline{F}(\xi_{\tau_n}))/(1-\tau_n))}{j!} [n(1-\tau_n)]^{(1-j)/2} M_{n,11}^{j/2} Z^j \right). $$

Recall that $\mathbb{E}(Z^j) = 0$ if $j$ is odd and $(j-1)!!$ if $j$ is even, where the double factorial $N!!$ denotes the product of all integers from 1 to $N$ having the same parity as $N$, i.e., when $N$ is even, $(N-1)!! = (N-1)(N-3)\cdots 3 \cdot 1 = N!/(2^{N/2}(N/2)!)$. It follows that

$$ \text{Cov}(Z^i, Z^j) = \begin{cases} 0 & \text{when } i+j \text{ is odd,} \\ (i+j-1)!! & \text{when } i \text{ and } j \text{ are odd,} \\ (i+j-1)!! - (i-1)!!(j-1)!! & \text{when } i \text{ and } j \text{ are even.} \end{cases} $$

Paired with the fact that $\psi^{(j)}(x) = (-1)^j j!(1+x)^{-j-1}$, we obtain, after straightforward calculations and gathering together negative powers of $n(1-\tau_n)$,

$$ 
\begin{aligned}
\mathfrak{M}_{n,11}(\infty) &= (1 + \overline{F}(\xi_{\tau_n})/(1 - \tau_n))^{-4} M_{n,11} \\
&\quad \times \left( 1 + \sum_{j=1}^{\infty} \frac{(1 + \overline{F}(\xi_{\tau_n})/(1 - \tau_n))^{-2j}}{[n(1 - \tau_n)]^j} M_{n,11}^j \left[ (2j + 1)!!(2j + 1) - \sum_{i=1}^{j} (2i - 1)!!(2j + 1 - 2i)!! \right] \right) \\
&= \sum_{j=0}^{\infty} \frac{\mathfrak{m}_{j,11}}{[n(1 - \tau_n)]^j}.
\end{aligned}
$$
We next seek $\mathfrak{M}_{n,12}(J)$, which requires evaluating

$$ \mathrm{Cov} \left( \left\{ \sqrt{n(1-\tau_n)} \left( \frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1-\tau_n} - \frac{\bar{F}(\xi_{\tau_n})}{1-\tau_n} \right) \right\}^j, \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1 \right), \quad \text{for } j \ge 1. $$

Recall that if $(X_1, X_2)$ is a Gaussian centered random pair then $\mathbb{E}(X_2|X_1) = (\mathrm{Cov}(X_1, X_2)/\mathrm{Var}(X_1))X_1$. It follows that, for any $j \ge 1$,

$$ \mathrm{Cov}(X_1^j, X_2) = \mathbb{E}(X_1^j X_2) = \frac{\mathrm{Cov}(X_1, X_2)}{\mathrm{Var}(X_1)} \mathbb{E}(X_1^{j+1}). $$

Therefore

$$ \begin{aligned} \mathfrak{M}_{n,12}(\infty) &= \sum_{j=1}^{\infty} \frac{\psi^{(j)}(\overline{F}(\xi_{\tau_n})/(1-\tau_n))}{j!} [n(1-\tau_n)]^{(1-j)/2} M_{n,12} M_{n,11}^{(j-1)/2} \mathbb{E}(Z^{j+1}) \\ &= -(1+\overline{F}(\xi_{\tau_n})/(1-\tau_n))^{-2} M_{n,12} \left( 1 + \sum_{j=1}^{\infty} \frac{(1+\overline{F}(\xi_{\tau_n})/(1-\tau_n))^{-2j}}{[n(1-\tau_n)]^j} M_{n,11}^j (2j+1)!! \right) \\ &=: \sum_{j=0}^{\infty} \frac{\mathfrak{m}_{j,12}}{[n(1-\tau_n)]^j}. \end{aligned} $$

The matrix $\mathfrak{M}(J) = \mathfrak{M}_n(J)$, approximating the covariance matrix of $\sqrt{n(1-\tau_n)}(\hat{\gamma}_{\tau_n}^E - \psi(\overline{F}(\xi_{\tau_n})/(1-\tau_n)), \hat{\xi}_{\tau_n}/\xi_{\tau_n} - 1)$, is finally obtained by truncating each of the series defining $\mathfrak{M}_{n,11}(\infty)$ and $\mathfrak{M}_{n,12}(\infty)$ at order $1/[n(1-\tau_n)]^J$. We now use this in order to approximate the uncertainty about $\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}$. One has

$$ \begin{aligned} \hat{\gamma}_{\tau_n}^{\mathrm{E,BR}} &= \left( 1 + (2\tau_n - 1)(1/\hat{\gamma}_{\tau_n}^E - 1) \left( 1 - \frac{\bar{Y}_n}{\hat{\xi}_{\tau_n}} \right)^{-1} (1 + \mathbb{O}_P(A((1-\tau_n)^{-1}))) \right)^{-1} \\ &= u_n(\hat{\gamma}_{\tau_n}^E, \hat{\xi}_{\tau_n}/\xi_{\tau_n}) + \mathbb{O}_P(A((1-\tau_n)^{-1})), \end{aligned} $$

where

$$ u_n(x, y) = \left( 1 + (2\tau_n - 1) \left( \frac{1}{x} - 1 \right) \left( 1 - \frac{\mathbb{E}(Y)}{y\xi_{\tau_n}} \right)^{-1} \right)^{-1}. $$

Here and in the rest of this section, we have neglected the finite-sample variability in $\bar{Y}_n$, which typically converges faster to $\mathbb{E}(Y)$ than the other (extreme value) terms do to their respective limits. We have also neglected any term proportional to $A((1-\tau_n)^{-1})$ or $A((1-\tau'_n)^{-1})$. The rationale behind this choice is that, since $A(t) = \gamma bt^\rho$, keeping these terms in this kind of calculation would then entail approximating the correlation of estimators of the second-order parameters $\rho$ and $b$ with estimators of other extreme value parameters, here $\hat{\gamma}_{\tau_n}^E$ and $\hat{\xi}_{\tau_n}$. This is a joint convergence problem which, to the best of our knowledge, remains open, and whose solution deserves a separate in-depth study. We now write

$$ \begin{aligned} & \sqrt{n(1-\tau_n)} \left( \hat{\gamma}_{\tau_n}^{\mathrm{E,BR}} - \gamma, \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1 \right) \\ &= \sqrt{n(1-\tau_n)} \left( u_n(\hat{\gamma}_{\tau_n}^E, \hat{\xi}_{\tau_n}/\xi_{\tau_n}) - u_n((1+\overline{F}(\xi_{\tau_n})/(1-\tau_n))^{-1}, 1), \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1 \right) + \mathbb{O}_P(1). \end{aligned} $$
The random vector

$$ \sqrt{n(1 - \tau_n)} \left( \hat{\gamma}_{\tau_n}^{\text{E,BR}} - \gamma, \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - 1 \right) $$

can thus be considered asymptotically Gaussian centered with a covariance matrix approximated by $\mathfrak{M}^{\text{BR}}(J) = \mathfrak{M}_n^{\text{BR}}(J)$ whose elements are

$$ \mathfrak{M}_{n,11}^{\text{BR}}(J) = (\partial_1 u_n((1 + \bar{F}(\xi_{\tau_n})/(1 - \tau_n))^{-1}, 1))^2 \mathfrak{M}_{n,11}(J) \\ + (\partial_2 u_n((1 + \bar{F}(\xi_{\tau_n})/(1 - \tau_n))^{-1}, 1))^2 \mathfrak{M}_{n,22}(J) \\ + 2\partial_1 u_n((1 + \bar{F}(\xi_{\tau_n})/(1 - \tau_n))^{-1}, 1)\partial_2 u_n((1 + \bar{F}(\xi_{\tau_n})/(1 - \tau_n))^{-1}, 1)\mathfrak{M}_{n,12}(J), $$

$$ \mathfrak{M}_{n,12}^{\text{BR}}(J) = \partial_1 u_n((1 + \bar{F}(\xi_{\tau_n})/(1 - \tau_n))^{-1}, 1)\mathfrak{M}_{n,12}(J) \\ + \partial_2 u_n((1 + \bar{F}(\xi_{\tau_n})/(1 - \tau_n))^{-1}, 1)\mathfrak{M}_{n,22}(J), $$

$$ \mathfrak{M}_{n,22}^{\text{BR}}(J) = \mathfrak{M}_{n,22}(J). $$

Here

$$ \partial_1 u_n(x, y) = \frac{\xi_{\tau_n}(2\tau_n - 1)(\xi_{\tau_n}y - \mathbb{E}(Y)y)}{[(\mathbb{E}(Y) - 2\xi_{\tau_n}y(1 - \tau_n))x - \xi_{\tau_n}y(2\tau_n - 1)]^2} \quad (B.8) $$

$$ \text{and } \partial_2 u_n(x, y) = \frac{\xi_{\tau_n}(2\tau_n - 1)\mathbb{E}(Y)x(1-x)}{[(\mathbb{E}(Y) - 2\xi_{\tau_n}y(1 - \tau_n))x - \xi_{\tau_n}y(2\tau_n - 1)]^2}. \quad (B.9) $$

Our final step is to combine all these elements in order to accurately quantify the uncertainty in $\log(\hat{\xi}_{\tau_n}^{\text{BR}}/\xi_{\tau_n}')$. We rewrite this quantity as

$$ \log \frac{\hat{\xi}_{\tau_n}^{\text{BR}}}{\xi_{\tau_n}'} = (\hat{\gamma}_{\tau_n}^{\text{E,BR}} - \gamma) \log \left( \frac{1 - \tau_n}{1 - \tau_n'} \right) + \log \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} - \log(1 + B_{1,n}) \\ + \hat{\gamma}_{\tau_n}^{\text{E,BR}} \log \left( \frac{1 - \bar{Y}_n}{2\tau_n - 1} \right) - \log(1 + B_{2,n}) \\ - \hat{\gamma}_{\tau_n}^{\text{E,BR}} \log \left( \frac{1 - \bar{Y}_n}{2\tau_n' - 1} \right) - \log(1 + B_{3,n}) \\ + \mathbb{O}_{\mathbb{P}}(A((1 - \tau_n)^{-1})). $$

Since we neglect any term proportional to $A((1 - \tau_n)^{-1})$ or $A((1 - \tau_n')^{-1})$, we now:

* Neglect $B_{1,n}$,

* Approximate $1 + B_{2,n}$ by $(1 + r(\tau_n))^{\gamma}$, where $1 + r(\tau_n)$ is itself approximated by $(1 - \mathbb{E}(Y)/\xi_{\tau_n})/(2\tau_n - 1)$,

* Approximate $1 + B_{3,n}$ by $(1 + r(\tau_n'))^{-\gamma}$, where $1 + r(\tau_n')$ is itself approximated by $(1 - \mathbb{E}(Y)/\xi_{\tau_n'})/(2\tau_n' - 1)$.
Hence the approximation

$$
\begin{aligned}
\log \frac{\hat{\xi}_{\tau_n}^{*,\mathrm{BR}}}{\xi_{\tau_n}} &\approx \left[ \log \left( \frac{1-\tau_n}{1-\tau_n'} \right) + \log \left( \frac{2\tau_n'-1}{2\tau_n-1} \right) \right] (\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}} - \gamma) + \log \frac{\hat{\xi}_{\tau_n}}{\xi_{\tau_n}} \\
&\quad + \hat{\gamma}_{\tau_n}^{\mathrm{E,BR}} \log \left( 1 - \frac{\mathbb{E}(Y)}{\hat{\xi}_{\tau_n}} \right) - \gamma \log \left( 1 - \frac{\mathbb{E}(Y)}{\xi_{\tau_n}} \right) \\
&\quad - \left[ \hat{\gamma}_{\tau_n}^{\mathrm{E,BR}} \log \left( 1 - \frac{\mathbb{E}(Y)}{\left(\frac{1-\tau_n'}{1-\tau_n}\right)^{-\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}} \hat{\xi}_{\tau_n}} \right) - \gamma \log \left( 1 - \frac{\mathbb{E}(Y)}{\xi_{\tau_n}} \right) \right].
\end{aligned}
$$

We finally write

$$
\begin{aligned}
& \frac{1}{\xi_{\tau_n'}} - \frac{1}{\left(\frac{1-\tau_n'}{1-\tau_n}\right)^{-\gamma} \xi_{\tau_n}} \\
&= \frac{1}{\xi_{\tau_n'}} \left( 1 - \frac{\xi_{\tau_n'}}{q_{\tau_n'}} \times \frac{q_{\tau_n'}}{\xi_{\tau_n}} \times \frac{q_{\tau_n'}}{\left(\frac{1-\tau_n'}{1-\tau_n}\right)^{-\gamma} q_{\tau_n}} \right) \\
&= O\left( \frac{1}{\xi_{\tau_n'}} \left( O(|A((1-\tau_n)^{-1})|) + O(|A((1-\tau_n')^{-1})|) + \frac{1}{\xi_{\tau_n}} + \frac{1}{\xi_{\tau_n'}} \right) \right)
\end{aligned}
$$

which is typically very small due to the presence of the factor $1/\xi_{\tau_n'}$ (the inverse of an expectile at the properly extreme level $\tau_n'$). This motivates approximating $\log(1 - \mathbb{E}(Y)/\xi_{\tau_n'})$ by

$$
\log \left( 1 - \frac{\mathbb{E}(Y)}{\left(\frac{1-\tau_n'}{1-\tau_n}\right)^{-\gamma} \xi_{\tau_n}} \right).
$$

This eventually leads to the approximation

$$
\frac{\sqrt{n(1-\tau_n)}}{\log\left(\frac{(1-\tau_n)}{(1-\tau_n')\tau_n'\tau_n}\right)} \log \frac{\hat{\xi}_{\tau_n}^{*,\mathrm{BR}}}{\xi_{\tau_n}} \approx \sqrt{n(1-\tau_n)} \left\{ g_n(\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}, \hat{\xi}_{\tau_n}/\xi_{\tau_n}) - g_n(\gamma, 1) \right\}
$$

where

$$
\begin{aligned}
g_n(x, y) = & x \left( 1 + \frac{\log\left(\frac{2\tau_n' - 1}{2\tau_n - 1}\right)}{\log\left(\frac{1 - \tau_n}{1 - \tau_n'}\right)} \right) + \frac{\log(y)}{\log\left(\frac{1 - \tau_n}{1 - \tau_n'}\right)} \\
& + \frac{x}{\log\left(\frac{1 - \tau_n}{1 - \tau_n'}\right)} \log\left(1 - \frac{\mathbb{E}(Y)}{\xi_{\tau_n} y}\right) \\
& - \frac{x}{\log\left(\frac{1 - \tau_n}{1 - \tau_n'}\right)} \log\left(1 - \frac{\mathbb{E}(Y)}{\left(\frac{1 - \tau_n'}{1 - \tau_n}\right)^{-x} \xi_{\tau_n} y}\right).
\end{aligned}
$$

The variance of $(\sqrt{n(1-\tau_n)}/\log((1-\tau_n)/(1-\tau_n'))\log(\hat{\xi}_{\tau_n}^{*,\mathrm{BR}}/\xi_{\tau_n}))$ is then well approximated by

$$
s_n^2(J) = (\partial_1 g_n(\gamma, 1))^2 \mathfrak{M}_{n,11}^{\mathrm{BR}}(J) + 2\partial_1 g_n(\gamma, 1)\partial_2 g_n(\gamma, 1) \mathfrak{M}_{n,12}^{\mathrm{BR}}(J) + (\partial_2 g_n(\gamma, 1))^2 \mathfrak{M}_{n,22}^{\mathrm{BR}}(J)
$$
where

$$
\begin{align}
& \partial_1 g_n(\gamma, 1) \nonumber \\
& = 1 + \frac{\log((2\tau'_n - 1)/(2\tau_n - 1))}{\log((1 - \tau_n)/(1 - \tau'_n))} + \frac{\log(1 - \mathbb{E}(Y)/\xi_{\tau_n})}{\log((1 - \tau_n)/(1 - \tau'_n))} \nonumber \\
& \quad - \frac{\log(1 - \mathbb{E}(Y)/((1 - \tau'_n)/(1 - \tau_n))^{-\gamma} \xi_{\tau_n})}{\log((1 - \tau_n)/(1 - \tau'_n))} - \frac{\gamma \mathbb{E}(Y)}{((1 - \tau'_n)/(1 - \tau_n))^{-\gamma} \xi_{\tau_n} - \mathbb{E}(Y)} \tag{B.10}
\end{align}
$$

and

$$
\begin{equation}
\begin{split}
\partial_2 g_n(\gamma, 1) ={}& \frac{1}{\log((1-\tau_n)/(1-\tau'_n))} \left( 1 - \frac{\gamma \mathbb{E}(Y)}{((1-\tau'_n)/(1-\tau_n))^{-\gamma} \xi_{\tau_n} - \mathbb{E}(Y)} + \frac{\gamma \mathbb{E}(Y)}{\xi_{\tau_n} - \mathbb{E}(Y)} \right) .
\end{split}
\tag{B.11}
\end{equation}
$$

In order to produce confidence intervals, we estimate $\gamma$ by $\hat{\gamma}_{\tau_n}^{\text{E,BR}}$, $\xi_{\tau_n}$ by $\hat{\xi}_{\tau_n}$, $\mathbb{E}(Y)$ by $\bar{Y}_n$, $\bar{F}(\xi_{\tau_n})$ by $\bar{F}_n(\hat{\xi}_{\tau_n})$, $\varphi^{(1)}(\xi_{\tau_n})$ by

$$
\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n}) = \frac{1}{n} \sum_{i=1}^{n} (Y_i - \hat{\xi}_{\tau_n}) \mathbb{1}\{Y_i > \hat{\xi}_{\tau_n}\}.
$$

The estimation of $\varphi^{(2)}(\xi_{\tau_n})$ is more complex, because its naive empirical counterpart
$\hat{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n})$ is unbiased but highly skewed, and therefore tends to vastly underestimate
$\varphi^{(2)}(\xi_{\tau_n})$. A second-order approximation of the underlying distribution function $\bar{F}$ in
a neighborhood of infinity suggests the expansion

$$
\varphi^{(2)}(\xi_{\tau_n}) \approx 2\bar{F}(\xi_{\tau_n})\xi_{\tau_n}^2 \gamma^2 \left( \frac{1}{(1-\gamma)(1-2\gamma)} + \frac{b\bar{F}(\xi_{\tau_n})^{-\rho}}{\rho} \left\{ \frac{1}{(1-\gamma-\rho)(1-2\gamma-\rho)} - \frac{1}{(1-\gamma)(1-2\gamma)} \right\} \right).
$$

A simple plug-in then suggests the estimator

$$
\tilde{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n}) = 2\hat{F}_n(\hat{\xi}_{\tau_n}) \hat{\xi}_{\tau_n}^2 (\hat{\gamma}_{\tau_n}^{\text{E,BR}})^2 \left( \frac{1}{(1-\hat{\gamma}_{\tau_n}^{\text{E,BR}})(1-2\hat{\gamma}_{\tau_n}^{\text{E,BR}})} + \frac{\bar{b}(\hat{F}_n(\hat{\xi}_{\tau_n}))^{-\bar{\rho}}}{\bar{\rho}} \left\{ \frac{1}{(1-\hat{\gamma}_{\tau_n}^{\text{E,BR}}-\bar{\rho})(1-2\hat{\gamma}_{\tau_n}^{\text{E,BR}}-\bar{\rho})} - \frac{1}{(1-\hat{\gamma}_{\tau_n}^{\text{E,BR}})(1-2\hat{\gamma}_{\tau_n}^{\text{E,BR}})} \right\} \right)
$$

which tends to be much more stable in finite samples than $\hat{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n})$. We deduce from these
calculations a corrected asymptotic Gaussian confidence interval for $\xi_{\tau_n'}$ at level $1 - \alpha$ as

$$
\hat{I}_{\tau_n'}^{(2)}(\alpha) = \hat{I}_{\tau_n'}^{(2)}(\alpha; J) = \left[ \hat{\xi}_{\tau_n'}^{\star, \text{BR}} \exp \left( \pm \frac{\log((1-\tau_n)/(1-\tau_n'))}{\sqrt{n(1-\tau_n)}} \sqrt{\hat{s}_n^2(J)} \times z_{1-\alpha/2} \right) \right]
$$

where $\hat{s}_n^2(J)$ is obtained from $s_n^2(J)$ by replacing $\gamma$ by $\hat{\gamma}_{\tau_n}^{\text{E,BR}}$, $\xi_{\tau_n}$ by $\hat{\xi}_{\tau_n}$, $\mathbb{E}(Y)$ by $\bar{Y}_n$, $\bar{F}(\xi_{\tau_n})$ by $\bar{F}_n(\hat{\xi}_{\tau_n})$, $\varphi^{(1)}(\xi_{\tau_n})$ by $\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n})$ and $\varphi^{(2)}(\xi_{\tau_n})$ by $\tilde{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n})$.
## B.2 Quantile-based confidence interval construction

### B.2.1 Preliminary steps

The quantile-based intermediate expectile estimator is $\tilde{\xi}_{\tau_n} = (1/\hat{\gamma}_{\tau_n}^H - 1)^{-\hat{\gamma}_{\tau_n}^H} \hat{q}_{\tau_n}$. Here we take $\tau_n = 1 - k_n/n$, where $k_n$ is a sequence of integers, $\hat{q}_{\tau_n} = Y_{n-\lfloor n(1-\tau_n) \rfloor, n} = Y_{n-k_n, n}$ is the corresponding intermediate order statistic, and

$$ \hat{\gamma}_{\tau_n}^H = \hat{\gamma}_{1-k_n/n}^H = \frac{1}{k_n} \sum_{i=1}^{k_n} \log Y_{n-i+1, n} - \log Y_{n-k_n, n} $$

is the usual Hill estimator of $\gamma$ calculated upon the top $k_n$ log-spacings. Then

$$ \log \left( \frac{\tilde{\xi}_{1-k_n/n}}{(\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n}} \right) = \phi(\hat{\gamma}_{1-k_n/n}^H) - \phi(\gamma) + \log \left( \frac{\hat{q}_{1-k_n/n}}{q_{1-k_n/n}} \right). \quad (B.12) $$

where $\phi(x) = -x \log(x^{-1} - 1)$. It is known that, if $\lambda_1 = \lim_{n \to \infty} \sqrt{k_n} A(n/k_n) \in \mathbb{R}$,

$$ \sqrt{k_n} \left( \hat{\gamma}_{1-k_n/n}^H - \gamma, \log \left( \frac{\hat{q}_{1-k_n/n}}{q_{1-k_n/n}} \right) \right) \xrightarrow{d} \mathcal{N} \left( \left( \frac{\lambda_1}{1-\rho}, 0 \right), \gamma^2 \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix} \right) $$

(combine e.g. Theorem 2.4.8 p.52, Lemma 3.2.3 p.71 and Theorem 3.2.5 p.74 in de Haan and Ferreira (2006)). The asymptotic correlation of the second term in (B.12) with the first one is then exactly 0, and the third term is a bias term that is evaluated by employing Proposition 1 in Daouia et al. (2020). A straightforward application of the delta-method with the function $g(x, y) = (x, y - x \log(x^{-1} - 1))$ leads to:

$$ \sqrt{k_n} \left( \hat{\gamma}_{1-k_n/n}^H - \gamma, \log \left( \frac{\tilde{\xi}_{1-k_n/n}}{(\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n}} \right) \right) \xrightarrow{d} \mathcal{N} \left( \frac{\lambda_1}{1-\rho} (1, m(\gamma)), \mathbf{V} \right), \quad (B.13) $$

where $\lambda_1 = \lim_{n \to \infty} \sqrt{k_n} A(n/k_n)$, $\lambda_2 = \lim_{n \to \infty} \sqrt{k_n}/q_{1-k_n/n}$,

$$ B_1 = \frac{(\gamma^{-1}-1)^{-\rho}}{1-\gamma-\rho} + \frac{(\gamma^{-1}-1)^{-\rho}-1}{\rho} \quad \text{and} \quad B_2 = \gamma(\gamma^{-1}-1)^{\gamma} \mathbb{E}(Y), $$

with $m(\gamma) = \phi'(\gamma) = (1 - \gamma)^{-1} - \log(\gamma^{-1} - 1)$, and $\mathbf{V}$ is the symmetric matrix having entries $V_{11} = \gamma^2$, $V_{12} = \gamma^2 m(\gamma)$ and $V_{22} = \gamma^2(1 + (m(\gamma))^2)$. Write finally

$$ \log \left( \frac{\tilde{\xi}_{\tau_n}^*}{\xi_{\tau_n}^*} \right) = (\hat{\gamma}_{1-k_n/n}^H - \gamma) \log \left( \frac{k_n}{n(1-\tau_n')} \right) + \log \left( \frac{\tilde{\xi}_{1-k_n/n}}{(\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n}} \right) \\ + \log \left( (\gamma^{-1}-1)^{-\gamma} \frac{q_{\tau_n'}}{\xi_{\tau_n}^*} \right) + \log \left( \left[ \frac{n(1-\tau_n')}{k_n} \right]^{-\gamma} \frac{q_{1-k_n/n}}{q_{\tau_n'}} \right). \quad (B.14) $$

Convergence (2.4) follows using Proposition 1 in Daouia et al. (2020) (for the first non-random term in the second line of (B.14)) and the equation at the top of p.139 in the proof of Theorem 4.3.8 in de Haan and Ferreira (2006) (for the second nonrandom term), and since $\bar{b} - b = o_p(1)$ and $(\bar{\rho} - \rho) \log(n) = o_p(1)$, $\tilde{\xi}_{\tau_n}^{*,BR}$ satisfies the same asymptotic normality result but with asymptotic mean 0.
## B.2.2 Detailed calculation

Convergence (2.4) is obtained by neglecting the finite-sample uncertainty in the estimator $\tilde{\xi}_{1-k_n/n}$. Combining (B.13) and (B.14) suggests that the asymptotic variance of

$$ \frac{\sqrt{k_n}}{\log(k_n/(n(1-\tau'_n)))} \log\left(\frac{\tilde{\xi}_{\tau'_n}^*}{\xi_{\tau'_n}}\right) $$

may be better approximated by

$$ \begin{aligned} & V_{11} + 2 \frac{V_{12}}{\log(k_n/(n(1-\tau'_n)))} + \frac{V_{22}}{\log^2(k_n/(n(1-\tau'_n)))} \\ &= \gamma^2 \left[ 1 + 2 \frac{m(\gamma)}{\log(k_n/(n(1-\tau'_n)))} + \frac{1+(m(\gamma))^2}{\log^2(k_n/(n(1-\tau'_n)))} \right]. \end{aligned} \quad (B.15) $$

This approximation substantially improves upon the expression of the variance obtained through first-order asymptotics, but behaves somewhat poorly when $\gamma$ is in a neighborhood of 1. An even more accurate approximation is possible by using the same trick as in Section B.1 and writing a formal power series expansion of $\phi$ instead. Use (B.12) in order to write

$$ \begin{aligned} & \sqrt{k_n} \log \left( \frac{\tilde{\xi}_{1-k_n/n}}{(\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n}} \right) \\ &= \sqrt{k_n} \log \left( \frac{\hat{q}_{1-k_n/n}}{q_{1-k_n/n}} \right) + \sum_{j=1}^{\infty} \frac{\phi^{(j)}(\gamma)}{j!} k_n^{\frac{1-j}{2}} \left\{ \sqrt{k_n} (\hat{\gamma}_{1-k_n/n}^{\mathrm{H}} - \gamma) \right\}^j. \end{aligned} $$

As in Section B.1, we approximate the variance of $\sqrt{k_n} \log(\tilde{\xi}_{1-k_n/n}/((\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n}))$ and its covariance with $\sqrt{k_n}(\hat{\gamma}_{1-k_n/n}^{\mathrm{H}} - \gamma)$ by working with the above right-hand side in a formal way and truncating the resulting (formal) variance and covariance power series at an order $1/k_n^J$, for a suitably chosen $J \ge 1$. This provides a finite-sample correction $V(J) = V_n(J)$ for the matrix $V = V(0)$ in (B.13), and hence a further correction for the asymptotic variance of the extrapolated quantile-based extreme expectile estimator. The resulting matrix $V(J)$ has the following coefficients: the first variance term $V_{11}(J) = V_{11}(0) = \gamma^2$, obtained from the asymptotic variance of $\sqrt{k_n}(\hat{\gamma}_{1-k_n/n}^{\mathrm{H}} - \gamma)$, does not depend on $J$. To calculate the second variance term $V_{22}(J)$, obtained by evaluating the variance of the right-hand side in the above expansion and Equation (B.13), recall first that $\sqrt{k_n}(\hat{\gamma}_{1-k_n/n}^{\mathrm{H}} - \gamma, \hat{q}_{1-k_n/n}/q_{1-k_n/n} - 1)$ converges to a Gaussian random pair with independent components having variance $\gamma^2$. We shall then obtain $V_{22}(J)$ after the calculation of

$$ V_{22}(\infty) = \gamma^2 + \mathrm{Var} \left( \sum_{j=1}^{\infty} \frac{\phi^{(j)}(\gamma)}{j!} k_n^{\frac{1-j}{2}} \gamma^j Z^j \right) $$
where Z is standard Gaussian. One finds

$$
\begin{align*}
V_{22}(\infty) &= \gamma^2 + \sum_{\substack{i,j \ge 1 \\ i+j \text{ even}}} \frac{\phi^{(i)}(\gamma)\phi^{(j)}(\gamma)(i+j-1)!!}{i!j!k_n^{(i+j)/2-1}} \gamma^{i+j} \\
&\quad - \sum_{\substack{i,j \ge 1 \\ i \text{ and } j \text{ even}}} \frac{\phi^{(i)}(\gamma)\phi^{(j)}(\gamma)}{i!j!k_n^{(i+j)/2-1}} (i-1)!!(j-1)!!\gamma^{i+j} \\
&= \gamma^2(1+(m(\gamma))^2) + \sum_{\substack{i,j \ge 1 \\ i+j \ge 2 \\ i+j \text{ even}}} \frac{\phi^{(i)}(\gamma)\phi^{(j)}(\gamma)(i+j-1)!!}{i!j!k_n^{(i+j)/2-1}} \gamma^{i+j} \\
&\quad - \sum_{i,j \ge 1} \frac{\phi^{(2i)}(\gamma)\phi^{(2j)}(\gamma)}{i!j!2^{i+j}k_n^{i+j-1}} \gamma^{2(i+j)}.
\end{align*}
$$

Gathering powers of $k_n$ together, we find

$$
\begin{align*}
V_{22}(\infty) &= \gamma^2(1 + (m(\gamma))^2) + \sum_{j=1}^{\infty} \frac{(2j+1)!!\gamma^{2j+2}}{k_n^j} \sum_{i=1}^{2j+1} \frac{\phi^{(i)}(\gamma)\phi^{(2j+2-i)}(\gamma)}{i!(2j+2-i)!} \\
&\quad - \sum_{j=1}^{\infty} \frac{\gamma^{2j+2}}{2^{j+1}k_n^j} \sum_{i=1}^{j} \frac{\phi^{(2i)}(\gamma)\phi^{(2j+2-2i)}(\gamma)}{i!(j+1-i)!} \\
&= \gamma^2 \left( 1 + (m(\gamma))^2 + \sum_{j=1}^{\infty} \frac{\gamma^{2j}}{k_n^j} \left( (2j+1)!! \sum_{i=1}^{2j+1} \frac{\phi^{(i)}(\gamma)\phi^{(2j+2-i)}(\gamma)}{i!(2j+2-i)!} \right. \right. \\
&\qquad \left. \left. - \frac{1}{2^{j+1}} \sum_{i=1}^{j} \frac{\phi^{(2i)}(\gamma)\phi^{(2j+2-2i)}(\gamma)}{i!(j+1-i)!} \right) \right) =: \sum_{j=0}^{\infty} \frac{v_{j,22}}{k_n^j}.
\end{align*}
$$

An expression of the formal covariance term $V_{12}(\infty)$ is similarly easily derived: if Z is
standard Gaussian,

$$
\begin{align*}
V_{12}(\infty) &= \sum_{j=0}^{\infty} \frac{\phi^{(j)}(\gamma)}{j!} k_n^{\frac{1-j}{2}} \gamma^{j+1} \mathbb{E}(Z^{j+1}) \\
&= \sum_{\substack{j \ge 0 \\ j \text{ odd}}} \frac{\phi^{(j)}(\gamma)}{j!} k_n^{\frac{1-j}{2}} \gamma^{j+1} j!! \\
&= \sum_{j=0}^{\infty} \frac{\phi^{(2j+1)}(\gamma)}{2^j j! k_n^j} \gamma^{2j+2} =: \sum_{j=0}^{\infty} \frac{v_{j,12}}{k_n^j}.
\end{align*}
$$

Closed-form expressions of the coefficients $v_{j,12}$ and $v_{j,22}$, for $j \ge 1$, can then be obtained
by remarking that

$$
\phi''(\gamma) = m'(\gamma) = \frac{1}{(1-\gamma)^2} + \frac{1}{\gamma(1-\gamma)} = \frac{1}{\gamma} + \frac{1}{1-\gamma} + \frac{1}{(1-\gamma)^2}
$$

and therefore

$$
\begin{align*}
\forall j \ge 2, \phi^{(j)}(\gamma) &= \frac{(-1)^{j-2}(j-2)!}{\gamma^{j-1}} + \frac{(j-2)!}{(1-\gamma)^{j-1}} + \frac{(j-1)!}{(1-\gamma)^j} \\
&= \frac{(-1)^{j-2}(j-2)! \sum_{i=0}^{j-2} \binom{j}{i} (-1)^i \gamma^i}{\gamma^{j-1}(1-\gamma)^j}.
\end{align*}
$$
In practice we truncate the power series at order $1/k_n^J$, for a given $J \ge 1$, so that a refined approximation to the covariance matrix of $\sqrt{k_n}(\hat{\gamma}_{1-k_n/n}^H - \gamma, \log(\tilde{\xi}_{1-k_n/n}/((\gamma^{-1} - 1)^{-\gamma}q_{1-k_n/n})))$ is the symmetric matrix $V(J)$, where the elements of $V(J)$ are

$$V_{11}(J) = \gamma^2, V_{12}(J) = \sum_{j=0}^{J} \frac{v_{j,12}}{k_n^j}, V_{22}(J) = \sum_{j=0}^{J} \frac{v_{j,22}}{k_n^j}.$$

In our implementation we stop the approximation at order $1/k_n^2$, that is, $J = 2$. A use of this refinement will be shown to dramatically improve the coverage probabilities of the asymptotic Gaussian confidence interval, especially in the challenging case when $\gamma$ is close to 1 and the right tail of $Y$ is very heavy. The explicit expressions of the coefficients $v_{j,12}$ and $v_{j,22}$, for $j = 0, 1, 2$, are provided in Table 3.

<table><caption>Table 3: First values of the coefficients $v_{j,12}/\gamma^2$ and $v_{j,22}/\gamma^2$, for $j = 0, 1, 2$.</caption><thead><tr><th>j</th><th>v<sub>j,12</sub>/γ<sup>2</sup></th><th>v<sub>j,22</sub>/γ<sup>2</sup></th></tr></thead><tbody><tr><td>0</td><td>m(γ)</td><td>(1 + (m(γ))<sup>2</sup>)</td></tr><tr><td>1</td><td>(3γ - 1) / (2(1 - γ)<sup>3</sup>)</td><td>(3γ - 1)m(γ) / (1 - γ)<sup>3</sup> + 1 / (2(1 - γ)<sup>4</sup>)</td></tr><tr><td>2</td><td>(3(10γ<sup>3</sup> - 10γ<sup>2</sup> + 5γ - 1) / (4(1 - γ)<sup>5</sup>))</td><td>(3(10γ<sup>3</sup> - 10γ<sup>2</sup> + 5γ - 1)m(γ) / (2(1 - γ)<sup>5</sup>) + (6γ<sup>2</sup> - 4γ + 1) / (1 - γ)<sup>6</sup> + (5(3γ - 1)<sup>2</sup> / (12(1 - γ)<sup>6</sup>))</td></tr></tbody></table>

We may now proceed with the construction of a confidence interval about $\xi_{\tau'_n}$. Recall Equation (2.6), which can be rewritten as

$$\log \frac{\tilde{\xi}_{\tau'_n}^{\star, \mathrm{BR}}}{\xi_{\tau'_n}} = \log \left( \frac{k_n}{n(1 - \tau'_n)} \right) (\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}} - \gamma) + \log \left( \frac{(1/\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}} - 1)^{-\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}} \hat{q}_{1-k_n/n}}{(1/\gamma - 1)^{-\gamma} q_{1-k_n/n}} \right) \\ + \log \left( (\gamma^{-1} - 1)^{-\gamma} \frac{q_{\tau'_n}}{\xi_{\tau'_n}} \right) + \log \left( \left[ \frac{n(1 - \tau'_n)}{k_n} \right]^{-\gamma} \frac{q_{1-k_n/n}}{q_{\tau'_n}} (1 + \tilde{B}_{1,n})(1 + \tilde{B}_{3,n}) \right).$$

Carrying out accurate inference about $\xi_{\tau'_n}$ based on $\tilde{\xi}_{\tau'_n}^{\star, \mathrm{BR}}$ requires providing a good approximation of the variance of the right-hand side. For that purpose, as in Section B.1, we neglect any term proportional to $A(n/k_n)$ or $A((1 - \tau'_n)^{-1})$. We therefore:

* Approximate $\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}$ by $\hat{\gamma}_{1-k_n/n}^{\mathrm{H}}$,

* Write

$$\log \left( (\gamma^{-1} - 1)^{-\gamma} \frac{q_{\tau'_n}}{\xi_{\tau'_n}} \right) \approx \gamma \log \left( 1 - \frac{\mathbb{E}(Y)}{\xi_{\tau'_n}} \right) - \gamma \log(2\tau'_n - 1),$$

see Equation (12) in Girard et al. (2022b),

* Neglect $\log([n(1 - \tau'_n)/k_n]^{-\gamma}q_{1-k_n/n}/q_{\tau'_n})$, in virtue of the equation at the top of p.139 in de Haan and Ferreira (2006),

* Neglect $\tilde{B}_{1,n}$,

* Approximate $1 + \tilde{B}_{3,n}$ by $(1 + \tilde{r}^*(\tau'_n))^{-\hat{\gamma}_{1-k_n/n}^H}$, where $1 + \tilde{r}^*(\tau'_n)$ is itself approximated by $(1 - \mathbb{E}(Y)/\tilde{\xi}_{\tau'_n}^*)/(2\tau'_n - 1)$.
This suggests the simpler asymptotic approximation

$$ \log \frac{\tilde{\xi}_{\tau_n'}^{\star, \text{BR}}}{\xi_{\tau_n'}} \approx \left( \log \left( \frac{k_n}{n(1-\tau_n')} \right) + \log(2\tau_n' - 1) \right) \left( \hat{\gamma}_{1-k_n/n}^{\text{H}} - \gamma \right) + \log \frac{\tilde{\xi}_{1-k_n/n}}{(\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n}} \\ - \left( \hat{\gamma}_{1-k_n/n}^{\text{H}} \log \left( 1 - \frac{\mathbb{E}(Y)}{\tilde{\xi}_{\tau_n'}^{\star}} \right) - \gamma \log \left( 1 - \frac{\mathbb{E}(Y)}{\xi_{\tau_n'}} \right) \right). $$

We finally write

$$ \frac{1}{\xi_{\tau_n'}} - \frac{1}{\left(\frac{n(1-\tau_n')}{k_n}\right)^{-\gamma} (\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n}} \\ = \frac{1}{\xi_{\tau_n'}} \left( 1 - (\gamma^{-1}-1)^{\gamma} \frac{\xi_{\tau_n'}}{q_{\tau_n'}} \times \frac{q_{\tau_n'}}{\left(\frac{n(1-\tau_n')}{k_n}\right)^{-\gamma} q_{1-k_n/n}} \right) \\ = O\left( \frac{1}{\xi_{\tau_n'}} \left( O(|A(n/k_n)|) + O(|A((1-\tau_n')^{-1})|) + \frac{1}{\xi_{\tau_n'}} \right) \right) $$

which is typically very small due to the presence of the factor $1/\xi_{\tau_n'}$ (the inverse of an expectile at the properly extreme level $\tau_n'$). This motivates approximating $\log(1 - \mathbb{E}(Y)/\xi_{\tau_n'})$ by

$$ \log \left( 1 - \frac{\mathbb{E}(Y)}{\left(\frac{n(1-\tau_n')}{k_n}\right)^{-\gamma} (\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n}} \right). $$

This eventually leads to the approximation

$$ \frac{\sqrt{k_n}}{\log(k_n/(n(1-\tau_n')))} \log \frac{\tilde{\xi}_{\tau_n'}^{\star, \text{BR}}}{\xi_{\tau_n'}} \approx \sqrt{k_n} \left\{ h_n \left( \hat{\gamma}_{1-k_n/n}^{\text{H}}, \frac{\tilde{\xi}_{1-k_n/n}}{(\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n}} \right) - h_n(\gamma, 1) \right\} $$

where

$$ h_n(x, y) = \left( 1 + \frac{\log(2\tau_n' - 1)}{\log(k_n/(n(1-\tau_n')))} \right) x + \frac{\log(y)}{\log(k_n/(n(1-\tau_n')))} \\ - \frac{x}{\log(k_n/(n(1-\tau_n')))} \log \left( 1 - \frac{\mathbb{E}(Y)}{\left(\frac{n(1-\tau_n')}{k_n}\right)^{-x} (\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n} y} \right). $$

This suggests that the variance of $(\sqrt{k_n}/\log(k_n/(n(1-\tau_n')))) \log(\tilde{\xi}_{\tau_n'}^{\star, \text{BR}}/\xi_{\tau_n'})$ is well approximated by

$$ \sigma_n^2(J) = (\partial_1 h_n(\gamma, 1))^2 V_{11}(J) + 2\partial_1 h_n(\gamma, 1)\partial_2 h_n(\gamma, 1)V_{12}(J) + (\partial_2 h_n(\gamma, 1))^2 V_{22}(J) $$

where

$$ \partial_1 h_n(\gamma, 1) = 1 - \frac{\gamma \mathbb{E}(Y)}{\left(\frac{n(1-\tau_n')}{k_n}\right)^{-\gamma} (\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n} - \mathbb{E}(Y)} \\ + \frac{1}{\log(k_n/(n(1-\tau_n')))} \left( \log(2\tau_n' - 1) - \log \left( 1 - \frac{\mathbb{E}(Y)}{\left(\frac{n(1-\tau_n')}{k_n}\right)^{-\gamma} (\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n}} \right) \right) \quad (B.16) $$
and

$$ \partial_2 h_n(\gamma, 1) = \frac{1}{\log(k_n/(n(1-\tau'_n)))} \left( 1 - \frac{\gamma \mathbb{E}(Y)}{\left(\frac{n(1-\tau'_n)}{k_n}\right)^{-\gamma} (\gamma^{-1}-1)^{-\gamma} q_{1-k_n/n} - \mathbb{E}(Y)} \right). \quad (B.17) $$

Note that this is indeed a further refinement of (B.15), since $\partial_1 h_n(\gamma, 1)$ converges to 1 and $\partial_2 h_n(\gamma, 1)$ is asymptotically equivalent to $1/\log(k_n/(n(1-\tau'_n))) \to 0$. We deduce from these calculations a corrected asymptotic Gaussian confidence interval for $\xi_{\tau'_n}$ at level $1-\alpha$ as

$$ \tilde{I}_{\tau'_n}^{(2)}(\alpha) = \tilde{I}_{\tau'_n}^{(2)}(\alpha; J) = \left[ \tilde{\xi}_{\tau'_n}^{\star, \mathrm{BR}} \exp \left( \pm \frac{\log(k_n/(n(1-\tau'_n)))}{\sqrt{k_n}} \sqrt{\tilde{\sigma}_n^2(J)} \times z_{1-\alpha/2} \right) \right] $$

where

$$ \tilde{\sigma}_n^2(J) = (\partial_1 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}, 1))^2 (\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}})^2 + (\partial_2 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}, 1))^2 \hat{V}_{22}(J) \\ + 2 \partial_1 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}, 1) \partial_2 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}, 1) \hat{V}_{12}(J) $$

is obtained by replacing $\gamma$ by $\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}$, $q_{1-k_n/n}$ by $\hat{q}_{1-k_n/n}$ and $\mathbb{E}(Y)$ by $\bar{Y}_n$ in the quantities $V_{11}(J) \equiv \gamma^2$, $V_{12}(J)$, $V_{22}(J)$, $\partial_1 h_n(\gamma, 1)$ and $\partial_2 h_n(\gamma, 1)$, thus producing their respective estimators $\hat{V}_{11}(J) \equiv (\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}})^2$, $\hat{V}_{12}(J)$, $\hat{V}_{22}(J)$, $\partial_1 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}, 1)$ and $\partial_2 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}, 1)$.

# C Theoretical results and proofs

## C.1 Auxiliary results

All distributions in Table 4 have natural location parameter 0. Some of them, such as the Burr and Dagum distributions, have scale parameters which never appear in the values of $\rho$ and $b$ because of the invariance of the ratio of two quantiles by changes of scale. Location parameters can be handled using the following lemma.

**Lemma C.1.** Assume that a given distribution has a survival function $\bar{F}$ and a probability density function $f$ (with respect to the Lebesgue measure) that satisfies

$$ f(t) = ct^{-\alpha-1}(1 + dt^{-\beta}(1 + o(1))) \quad \text{as } t \to \infty $$

where $\alpha, \beta, c > 0$ and $d \neq 0$. Set $\gamma = 1/\alpha$.

(i) Then $\bar{F}$ satisfies condition $C_2(\gamma, -\beta/\alpha, A)$, where

$$ A(t) = \left( -\frac{d\beta}{\alpha + \beta} (c/\alpha)^{-\beta/\alpha} \right) \times \frac{1}{\alpha} \times t^{-\beta/\alpha} = b\gamma t^{-\beta/\alpha} $$

$$ \text{with } b = -\frac{d\beta}{\alpha + \beta} (c/\alpha)^{-\beta/\alpha}. $$

Assume moreover that

$$ f(t) = ct^{-\alpha-1}(1 + dt^{-\beta} + d't^{-\beta-\beta'}(1 + o(1))) \quad \text{as } t \to \infty $$

where $\alpha, \beta, \beta', c > 0$ and $d, d' \neq 0$. Let $m \neq 0$ and define a survival function $\bar{G}$ having probability density function $g: t \mapsto f(t-m)$. Then:
(ii) If $\beta < 1$, then $\bar{G}$ satisfies condition $C_2(\gamma, -\beta/\alpha, A)$, where A is as in (i),

(iii) If $\beta > 1$, then $\bar{G}$ satisfies condition $C_2(\gamma, -1/\alpha, B_1)$, where

$$B_1(t) = -m(c/\alpha)^{-1/\alpha} \times \frac{1}{\alpha} \times t^{-1/\alpha} = b_1\gamma t^{-1/\alpha}$$
$$\text{with } b_1 = -m(c/\alpha)^{-1/\alpha}.$$

(iv) If $\beta = 1$ and $(\alpha + 1)m + d \neq 0$, then $\bar{G}$ satisfies condition $C_2(\gamma, -1/\alpha, B_2)$, where

$$B_2(t) = \left( -\frac{(\alpha + 1)m + d}{\alpha + 1} (c/\alpha)^{-1/\alpha} \right) \times \frac{1}{\alpha} \times t^{-1/\alpha} = b_2 \gamma t^{-1/\alpha}$$
$$\text{with } b_2 = -\frac{(\alpha + 1)m + d}{\alpha + 1} (c/\alpha)^{-1/\alpha}.$$

(v) If $\beta = 1$ and $(\alpha + 1)m + d = 0$, but $\beta' < 1$, then $\bar{G}$ satisfies condition $C_2(\gamma, -(1 + \beta')/\alpha, B_3)$, where

$$B_3(t) = \left( -\frac{d'(1 + \beta')}{\alpha + \beta' + 1} (c/\alpha)^{-(1+\beta')/\alpha} \right) \times \frac{1}{\alpha} \times t^{-(1+\beta')/\alpha} = b_3 \gamma t^{-(1+\beta')/\alpha}$$
$$\text{with } b_3 = -\frac{d'(1 + \beta')}{\alpha + \beta' + 1} (c/\alpha)^{-(1+\beta')/\alpha}.$$

(vi) If $\beta = 1$ and $(\alpha + 1)m + d = 0$, but $\beta' > 1$, then $\bar{G}$ satisfies condition $C_2(\gamma, -2/\alpha, B_4)$, where

$$B_4(t) = \frac{d^2}{\alpha + 1} (c/\alpha)^{-2/\alpha} \times \frac{1}{\alpha} \times t^{-2/\alpha} = b_4 \gamma t^{-2/\alpha}$$
$$\text{with } b_4 = \frac{d^2}{\alpha + 1} (c/\alpha)^{-2/\alpha}.$$

(vii) If $\beta = 1$, $(\alpha + 1)m + d = 0$, $\beta' = 1$ and $2d'(\alpha + 1) \neq d^2(\alpha + 2)$, then $\bar{G}$ satisfies condition $C_2(\gamma, -2/\alpha, B_5)$, where

$$B_5(t) = \left( -\frac{2}{\alpha + 2} \left[ d' - d^2 \frac{\alpha + 2}{2(\alpha + 1)} \right] (c/\alpha)^{-2/\alpha} \right) \times \frac{1}{\alpha} \times t^{-2/\alpha} = b_5 \gamma t^{-2/\alpha}$$
$$\text{with } b_5 = -\frac{2}{\alpha + 2} \left[ d' - d^2 \frac{\alpha + 2}{2(\alpha + 1)} \right] (c/\alpha)^{-2/\alpha}.$$

Heavy-tailed distributions with negative second-order parameter typically satisfy the third-order expansion

$$f(t) = ct^{-\alpha-1}(1 + dt^{-\beta} + d't^{-\beta-\beta'}(1 + o(1))) \text{ as } t \to \infty$$

where $\alpha, \beta, \beta', c > 0$ and $d, d' \neq 0$, as well as (depending on $m$) one of the conditions in (ii)-(vii). The only exception among the usual families of heavy-tailed distributions seems to be the Generalized Pareto distribution, for which it can be the case that $f(t-m) = ct^{-\alpha-1}$ for suitable $\alpha, c > 0$, and thus no remainder term is present. This, however, represents the ideal case when the shifted distribution is exactly Pareto and no bias is incurred by the use of the heavy tail assumption.
Proof. The assumption

$$f(t) = ct^{-\alpha-1}(1 + dt^{-\beta}(1 + o(1))) \text{ as } t \to \infty$$

yields

$$\bar{F}(t) = \frac{c}{\alpha} t^{-\alpha} \left( 1 + \frac{d\alpha}{\alpha + \beta} t^{-\beta} (1 + o(1)) \right).$$

In particular, for any $x > 0$,

$$\frac{\bar{F}(tx)}{\bar{F}(t)} - x^{-\alpha} = \left(-\frac{d\beta}{\alpha + \beta}\right) \frac{1}{\alpha} t^{-\beta} \times x^{-\alpha} \frac{x^{-\beta} - 1}{(-\beta/\alpha) \times (1/\alpha)} + o(t^{-\beta}).$$

It is then straightforward to check statement (i). The key to showing statements (ii)-(vii) is to note that when

$$f(t) = ct^{-\alpha-1}(1 + dt^{-\beta} + d't^{-\beta-\beta'}(1 + o(1))) \text{ as } t \to \infty,$$

it similarly holds that

$$\begin{aligned} f(t-m) &= c(t-m)^{-\alpha-1}(1+d(t-m)^{-\beta}+d't^{-\beta-\beta'}(1+o(1))) \\ &= ct^{-\alpha-1}(1-mt^{-1})^{-\alpha-1}(1+dt^{-\beta}(1-mt^{-1})^{-\beta}+d't^{-\beta-\beta'}(1+o(1))) \\ &= ct^{-\alpha-1}\left(1+(\alpha+1)mt^{-1}+\frac{(\alpha+1)(\alpha+2)m^2}{2}t^{-2}+o(t^{-2})\right) \\ &\quad \times (1+dt^{-\beta}+\beta dmt^{-\beta-1}+d't^{-\beta-\beta'}+o(t^{-\beta-1})+o(t^{-\beta-\beta'})) \\ &= ct^{-\alpha-1}\left(1+(\alpha+1)mt^{-1}+dt^{-\beta}\right. \\ &\qquad \left. +\frac{(\alpha+1)(\alpha+2)m^2}{2}t^{-2}+(\alpha+\beta+1)dmt^{-\beta-1}+d't^{-\beta-\beta'}+o(t^{-2})+o(t^{-\beta-1})+o(t^{-\beta-\beta'})\right). \end{aligned}$$

We then find that, as $t \to \infty$:

* In the setup of (ii),

$$f(t-m) = ct^{-\alpha-1}(1 + dt^{-\beta}(1 + o(1))).$$

* In the setup of (iii),

$$f(t-m) = ct^{-\alpha-1}(1 + (\alpha + 1)mt^{-1}(1 + o(1))).$$

* In the setup of (iv),

$$f(t-m) = ct^{-\alpha-1}(1 + ((\alpha+1)m + d)t^{-1}(1 + o(1))).$$

* In the setup of (v),

$$f(t-m) = ct^{-\alpha-1}(1 + d't^{-1-\beta'}(1 + o(1))).$$

* In the setup of (vi),

$$f(t-m) = ct^{-\alpha-1}\left(1 - \frac{d^2(\alpha+2)}{2(\alpha+1)}t^{-2}(1 + o(1))\right).$$

* In the setup of (vii),

$$f(t-m) = ct^{-\alpha-1}\left(1 + \left(d' - \frac{d^2(\alpha+2)}{2(\alpha+1)}\right) t^{-2}(1 + o(1))\right).$$

The result then immediately follows from (i). □
<table><caption>Table 4: A list of standard continuous heavy-tailed distributions satisfying $C_2(γ, ρ, A)$ with $A(t) = bγt^ρ$, with the associated values of $γ$, $ρ$ and $b$. The Burr distribution in this parametrization is also called the Singh-Maddala distribution, and it coincides with the log-logistic distribution when $β = 1$. The Dagum distribution is also called Burr type III distribution on p.59 of Beirlant et al. (2004). $Γ$ stands for Euler's Gamma function, and in the density function of the Davis distribution, $ζ$ denotes Riemann's Zeta function. In the Slash density function, $χ_{1+α}^2$ stands for the cumulative distribution function of the $χ^2$ distribution with $1+α$ degrees of freedom; the definition is derived from the univariate version of the slash distribution in Wang and Genton (2006).</caption><thead><tr><th>Distribution (parameters)</th><th>Density function</th><th>γ</th><th>ρ</th><th>b</th></tr></thead><tbody><tr><td>Pareto (α > 0)</td><td>αt<sup>-α-1</sup> (t > 1)</td><td>1/α</td><td>-∞</td><td>0 (by convention)</td></tr><tr><td>Hall-Weiss (α, β > 0)</td><td>(αt<sup>-α-1</sup> + (α + β)t<sup>-α-β-1</sup>)/2 (t > 1)</td><td>1/α</td><td>-β/α</td><td>-β 2<sup>β/α</sup>/α</td></tr><tr><td>Burr (α, β, λ > 0)</td><td>αβ(1/λ)<sup>α-1</sup>(1 + (t/λ)<sup>α</sup>)<sup>-β-1</sup> (t > 0)</td><td>1/(αβ)</td><td>-1/β</td><td>1</td></tr><tr><td>Dagum (α, β, λ > 0)</td><td>αβ(1/λ)<sup>-α-1</sup>(1 + (t/λ)<sup>-α</sup>)<sup>-β-1</sup> (t > 0)</td><td>1/α</td><td>-1</td><td>β + 1/2β</td></tr><tr><td>Fréchet (α > 0)</td><td>αt<sup>-α-1</sup> exp(-t<sup>-α</sup>) (t > 0)</td><td>1/α</td><td>-1</td><td>1/2</td></tr><tr><td>Generalized Pareto (σ, ξ > 0)</td><td>σ<sup>-1</sup>(1 + ξt/σ)<sup>-1-1/ξ</sup> (t > 0)</td><td>ξ</td><td>-ξ</td><td>1</td></tr><tr><td>Fisher (ν₁, ν₂ > 0)</td><td>(ν₁/ν₂)<sup>ν₁/2</sup> t<sup>ν₁/2-1</sup>(1 + ν₁t/ν₂)<sup>-ν₁+ν₂/2</sup> (t > 0)</td><td>2/ν₂</td><td>-2/ν₂</td><td>ν₁ + ν₂(B(ν₁/2, ν₂/2))<sup>2/ν₂</sup>(ν₂/2)<sup>2/ν₂</sup></td></tr><tr><td>Lévy (λ > 0)</td><td>√(λ/2π)t<sup>-3/2</sup> exp(-λ/(2t)) (t > 0)</td><td>2</td><td>-2</td><td>π/6</td></tr><tr><td>Inverse-G (α, λ > 0)</td><td>λα(1/α)t<sup>-α-1</sup> exp(-λ/t) (t > 0)</td><td>1/α</td><td>-1/α</td><td>(Γ(α+1))<sup>1/α</sup>/(α+1)</td></tr><tr><td>Davis (α > 1, λ > 0)</td><td>λαt<sup>-α-1</sup>(1/α)<sup>α/(α-1)</sup>(e<sup>λt</sup> - 1) (t > 0)</td><td>1/(α-1)</td><td>-1/(α-1)</td><td>1/(α-1)Γ(α)<sup>α/(α-1)</sup></td></tr><tr><td>Slash (α > 0)</td><td>2<sup>α/2-1</sup>αΓ(α+1/2)χ<sup>2</sup>(t<sup>2</sup>)t<sup>-1-α</sup> (t ∈ ℝ \ {0})</td><td>1/α</td><td>-∞</td><td>0 (by convention)</td></tr><tr><td>Cauchy (c > 0)</td><td>c/(π(c² + t<sup>2</sup>)) (t ∈ ℝ)</td><td>1</td><td>-2</td><td>2π²/3</td></tr><tr><td>Student (ν > 0)</td><td>1/√(νπ)Γ(ν/2) (1 + t<sup>2</sup>/ν)<sup>-(ν+1)/2</sup> (t ∈ ℝ)</td><td>1/ν</td><td>-2/ν</td><td>ν(ν+1)/(ν+2) (Γ(ν+1/2)ν<sup>(ν-1)/2</sup>)/(νπΓ(ν/2))<sup>-2/ν</sup></td></tr></tbody></table>
The following lemma is a law of large numbers, with rate of convergence, for $\hat{F}_n$ and $\hat{\varphi}_n^{(1)}$ at random thresholds. It is the crucial element for our proof of Theorem 1(i), with the quantification of the rate being a prerequisite for the proof of Lemma C.3 below, which is itself key to the proof of Theorem 1(ii). Before that we recall the following fact: if condition $C_2(\gamma, \rho, A)$ is satisfied and $\gamma < 1/k$, where $k$ is a positive integer, then

$$ 
\begin{aligned}
\varphi^{(k)}(x) &= kx^k \overline{F}(x) \left( \int_1^\infty (v-1)^{k-1} v^{-1/\gamma} dv + \int_1^\infty (v-1)^{k-1} \left\{ \frac{\overline{F}(vx)}{\overline{F}(x)} - v^{-1/\gamma} \right\} dv \right) \\
&= x^k \overline{F}(x) \left( \frac{1}{\gamma} B(k+1, \gamma^{-1}-k) \right. \\
&\qquad \left. + \frac{A(1/\overline{F}(x))}{\gamma \rho} \left\{ \frac{1-\rho}{\gamma} B(k+1, \gamma^{-1}(1-\rho)-k) - \frac{1}{\gamma} B(k+1, \gamma^{-1}-k) + o(1) \right\} \right) \quad (C.1)
\end{aligned}
$$

as $x \to \infty$, using uniform second-order regular variation inequalities: see Theorem B.2.18 in de Haan and Ferreira (2006), which applies since second-order regular variation can be rewritten as extended regular variation (this is explicitly stated as Lemma 5 in Daouia et al. (2020)). In this equation, $B(\cdot, \cdot)$ denotes the Beta function. In particular,

$$ \varphi^{(k)}(x) = \frac{B(k+1, \gamma^{-1}-k)}{\gamma} x^k \overline{F}(x) (1 + O(A(1/\overline{F}(x)))) \quad \text{as } x \to \infty. \quad (C.2) $$

**Lemma C.2.** Assume that condition $C_2(\gamma, \rho, A)$ is satisfied. Let $Y_1, \dots, Y_n$ be independent random variables with distribution function $F$, and assume that $\tau_n \uparrow 1$ satisfies $n(1-\tau_n) \to \infty$ and $\sqrt{n(1-\tau_n)}A((1-\tau_n)^{-1}) = O(1)$ as $n \to \infty$. Let finally $\hat{c}_n$ be a random sequence converging in probability to a fixed constant $c > 0$ at the rate $1/\sqrt{n(1-\tau_n)}$, namely, $\hat{c}_n - c = O_P(1/\sqrt{n(1-\tau_n)})$.

(i) Then

$$ \frac{\hat{F}_n(\hat{c}_n q_{\tau_n})}{\overline{F}(cq_{\tau_n})} = \frac{1}{n} \sum_{i=1}^{n} \frac{1\{Y_i > \hat{c}_n q_{\tau_n}\}}{\mathbb{P}(Y > cq_{\tau_n})} = 1 + O_P(1/\sqrt{n(1-\tau_n)}). $$

(ii) If $0 < \gamma < 1/2$ then

$$ \frac{\hat{\varphi}_n^{(1)}(\hat{c}_n q_{\tau_n})}{\varphi^{(1)}(cq_{\tau_n})} = \frac{1}{n} \sum_{i=1}^{n} \frac{(Y_i - \hat{c}_n q_{\tau_n}) 1\{Y_i > \hat{c}_n q_{\tau_n}\}}{\mathbb{E}((Y - cq_{\tau_n}) 1\{Y > cq_{\tau_n}\})} = 1 + O_P(1/\sqrt{n(1-\tau_n)}). $$

*Proof of Lemma C.2.* (i) Obviously

$$ \frac{\hat{F}_n(cq_{\tau_n})}{\overline{F}(cq_{\tau_n})} = 1 + O_P(1/\sqrt{n(1-\tau_n)}). \quad (C.3) $$

This is immediate because the left-hand side is a mean of independent and identically distributed random variables for each $n$, having mean 1, and

$$ \mathrm{Var} \left( \frac{\hat{F}_n(cq_{\tau_n})}{\overline{F}(cq_{\tau_n})} \right) = \frac{1}{n} \times \frac{\mathrm{Var}(1\{Y > cq_{\tau_n}\})}{(\overline{F}(cq_{\tau_n}))^2} = \frac{1}{n} \left( \frac{1}{\overline{F}(cq_{\tau_n})} - 1 \right) $$

is asymptotically equivalent to a multiple of $1/(n(1-\tau_n))$, by the regular variation assumption on $\overline{F}$. It then suffices to check that

$$ \frac{\hat{F}_n(\hat{c}_n q_{\tau_n})}{\hat{F}_n(cq_{\tau_n})} = 1 + O_P(1/\sqrt{n(1-\tau_n)}). $$
Write

$$ \left| \frac{\hat{F}_n(\hat{c}_n q_{\tau_n})}{\hat{F}_n(cq_{\tau_n})} - 1 \right| \le \frac{1 + o_P(1)}{\bar{F}(cq_{\tau_n})} \times \frac{1}{n} \sum_{i=1}^{n} |\mathbb{1}\{Y_i > \hat{c}_n q_{\tau_n}\} - \mathbb{1}\{Y_i > cq_{\tau_n}\}|. $$

It follows from $\hat{c}_n - c = O_P(1/\sqrt{n(1-\tau_n)})$ and (C.3) that, for any $\varepsilon > 0$, there is $M > 0$ such that, for $n$ large enough, one has

$$ \begin{aligned} & \left| \frac{\hat{F}_n(\hat{c}_n q_{\tau_n})}{\hat{F}_n(cq_{\tau_n})} - 1 \right| \\ & \le \frac{2}{\bar{F}(cq_{\tau_n})} \times \frac{1}{n} \sum_{i=1}^{n} (\mathbb{1}\{Y_i > c(1 - M/\sqrt{n(1-\tau_n)})q_{\tau_n}\} - \mathbb{1}\{Y_i > c(1 + M/\sqrt{n(1-\tau_n)})q_{\tau_n}\}) \\ & = \frac{2(\hat{F}_n(c(1 - M/\sqrt{n(1-\tau_n)})q_{\tau_n}) - \hat{F}_n(c(1 + M/\sqrt{n(1-\tau_n)})q_{\tau_n}))}{\bar{F}(cq_{\tau_n})} \\ & = \frac{2(\bar{F}(c(1 - M/\sqrt{n(1-\tau_n)})q_{\tau_n}) - \bar{F}(c(1 + M/\sqrt{n(1-\tau_n)})q_{\tau_n}))}{\bar{F}(cq_{\tau_n})} + O_P(1/\sqrt{n(1-\tau_n)}) \end{aligned} $$

with probability larger than $1 - \varepsilon/2$. Now

$$ \frac{2(\bar{F}(c(1 - M/\sqrt{n(1 - \tau_n)})q_{\tau_n}) - \bar{F}(c(1 + M/\sqrt{n(1 - \tau_n)})q_{\tau_n}))}{\bar{F}(cq_{\tau_n})} = O(1/\sqrt{n(1 - \tau_n)}) $$

using the assumption $\sqrt{n(1-\tau_n)}A((1-\tau_n)^{-1}) = O(1)$ and the locally uniform character of second-order regular variation (see again Theorem B.2.18 in de Haan and Ferreira (2006)). Conclude that there is $M' > 0$ such that, for $n$ large enough, one has

$$ \sqrt{n(1 - \tau_n)} \left| \frac{\hat{F}_n(\hat{c}_n q_{\tau_n})}{\hat{F}_n(cq_{\tau_n})} - 1 \right| \le M' $$

with probability larger than $1 - \varepsilon$. The result follows.

(ii) We mimic the proof of (i) with a couple of adaptations. First of all

$$ \frac{\hat{\varphi}_n^{(1)}(cq_{\tau_n})}{\varphi^{(1)}(cq_{\tau_n})} = 1 + O_P(1/\sqrt{n(1-\tau_n)}). \qquad (C.4) $$

Indeed, the left-hand side is a mean of independent and identically distributed random variables for each $n$, having mean 1, and

$$ \mathrm{Var} \left( \frac{\hat{\varphi}_n^{(1)}(cq_{\tau_n})}{\varphi^{(1)}(cq_{\tau_n})} \right) = \frac{1}{n} \times \frac{\mathrm{Var}((Y - cq_{\tau_n})\mathbb{1}\{Y > cq_{\tau_n}\})}{(\varphi^{(1)}(cq_{\tau_n}))^2} = \frac{1}{n} \left( \frac{\varphi^{(2)}(cq_{\tau_n})}{(\varphi^{(1)}(cq_{\tau_n}))^2} - 1 \right) $$

which, by Lemma 3(i) in Stupfler and Usseglio-Carleve (2023), is asymptotically equivalent to a multiple of $1/(n(1-\tau_n))$. It is then enough to prove that

$$ \frac{\hat{\varphi}_n^{(1)}(\hat{c}_n q_{\tau_n})}{\hat{\varphi}_n^{(1)}(cq_{\tau_n})} = 1 + O_P(1/\sqrt{n(1-\tau_n)}). $$
Note then that

$$
\begin{align*}
\frac{\hat{\varphi}_n^{(1)}(\hat{c}_n q_{\tau_n})}{\hat{\varphi}_n^{(1)}(cq_{\tau_n})} - 1
&= \frac{1 + o_{\mathbb{P}}(1)}{\varphi^{(1)}(cq_{\tau_n})} \times \frac{1}{n} \sum_{i=1}^{n} (Y_i - cq_{\tau_n}) (\mathbb{1}\{Y_i > \hat{c}_n q_{\tau_n}\} - \mathbb{1}\{Y_i > cq_{\tau_n}\}) \\
&\quad - (\hat{c}_n - c) \frac{q_{\tau_n} \bar{F}(cq_{\tau_n})}{\varphi^{(1)}(cq_{\tau_n})} \times \frac{\hat{F}_n(\hat{c}_n q_{\tau_n})}{\bar{F}(cq_{\tau_n})} (1 + o_{\mathbb{P}}(1)) \\
&= \frac{1 + o_{\mathbb{P}}(1)}{\varphi^{(1)}(cq_{\tau_n})} \times \frac{1}{n} \sum_{i=1}^{n} (Y_i - cq_{\tau_n}) (\mathbb{1}\{Y_i > \hat{c}_n q_{\tau_n}\} - \mathbb{1}\{Y_i > cq_{\tau_n}\}) + O_{\mathbb{P}}(1/\sqrt{n(1-\tau_n)})
\end{align*}
$$

by (i) of the present lemma and Lemma 3(i) in Stupfler and Usseglio-Carleve (2023). Then,
for any $\varepsilon > 0$, there is $M > 0$ such that, for $n$ large enough, one has

$$
\begin{align*}
\left| \frac{1}{\varphi^{(1)}(cq_{\tau_n})} \times \frac{1}{n} \sum_{i=1}^{n} (Y_i - cq_{\tau_n}) (\mathbb{1}\{Y_i > \hat{c}_n q_{\tau_n}\} - \mathbb{1}\{Y_i > cq_{\tau_n}\}) \right|
&\le \frac{1}{\varphi^{(1)}(cq_{\tau_n})} \times \frac{1}{n} \sum_{i=1}^{n} |Y_i - cq_{\tau_n}| (\mathbb{1}\{Y_i > c(1 - M/\sqrt{n(1-\tau_n)})q_{\tau_n}\} - \mathbb{1}\{Y_i > c(1 + M/\sqrt{n(1-\tau_n)})q_{\tau_n}\}) \\
&\le \frac{1}{\varphi^{(1)}(cq_{\tau_n})} \left( \begin{aligned}[t]
& \frac{1}{n} \sum_{i=1}^{n} (Y_i - c(1 - M/\sqrt{n(1-\tau_n)})q_{\tau_n}) \mathbb{1}\{Y_i > c(1 - M/\sqrt{n(1-\tau_n)})q_{\tau_n}\} \\
& - \frac{1}{n} \sum_{i=1}^{n} (Y_i - c(1 + M/\sqrt{n(1-\tau_n)})q_{\tau_n}) \mathbb{1}\{Y_i > c(1 + M/\sqrt{n(1-\tau_n)})q_{\tau_n}\} \\
& + \frac{M}{\sqrt{n(1-\tau_n)}} cq_{\tau_n} \times \frac{1}{n} \sum_{i=1}^{n} (\mathbb{1}\{Y_i > c(1 - M/\sqrt{n(1-\tau_n)})q_{\tau_n}\} - \mathbb{1}\{Y_i > c(1 + M/\sqrt{n(1-\tau_n)})q_{\tau_n}\})
\end{aligned} \right) \\
&= \frac{\varphi^{(1)}(c(1 - M/\sqrt{n(1-\tau_n)})q_{\tau_n})}{\varphi^{(1)}(cq_{\tau_n})} - \frac{\varphi^{(1)}(c(1 + M/\sqrt{n(1-\tau_n)})q_{\tau_n})}{\varphi^{(1)}(cq_{\tau_n})} \\
&\quad + \frac{M}{\sqrt{n(1-\tau_n)}} cq_{\tau_n} \left\{ \frac{\bar{F}(c(1 - M/\sqrt{n(1-\tau_n)})q_{\tau_n})}{\varphi^{(1)}(cq_{\tau_n})} - \frac{\bar{F}(c(1 + M/\sqrt{n(1-\tau_n)})q_{\tau_n})}{\varphi^{(1)}(cq_{\tau_n})} \right\} \\
&\quad + O_{\mathbb{P}}(1/\sqrt{n(1-\tau_n)})
\end{align*}
$$

with probability larger than $1 - \varepsilon/2$ (where (C.3) and (C.4) were used again). Using (C.2) and again the locally uniform character of second-order regular variation, we find that the upper bound is a $O_{\mathbb{P}}(1/\sqrt{n(1-\tau_n)})$, thus completing the proof. $\square$

Our final lemma quantifies the error made in the estimation of the various matrices intro-
duced in the construction of our refined confidence interval for the LAWS estimator.

**Lemma C.3.** Assume that $\mathbb{E}(|\min(Y, 0)|^{2+\delta}) < \infty$ for some $\delta > 0$, and condition $C_2(\gamma, \rho, A)$ is satisfied with $0 < \gamma < 1/2$, $\rho < 0$ and $A(t) = b\gamma t^\rho$. Let $\tau_n \uparrow 1$ satisfy $n(1-\tau_n) \to \infty$ and $\sqrt{n(1-\tau_n)}(1/q_{\tau_n} + A((1-\tau_n)^{-1})) = O(1)$ as $n \to \infty$. Throughout this lemma, let $\hat{M}_n^\varphi, \hat{M}_n^E, \hat{M}_n^\xi, \hat{M}_n$ and $\hat{M}_n(J)$ be the random matrices deduced from $M_n^\varphi, M_n^E, M_n^\xi, M_n$ and $M_n(J)$, respectively, by replacing $\gamma$ by $\hat{\gamma}_{\tau_n}^{\text{E,BR}}$, $\xi_{\tau_n}$ by $\hat{\xi}_{\tau_n}$, $\mathbb{E}(Y)$ by $\bar{Y}_n$, $\bar{F}(\xi_{\tau_n})$ by $\hat{F}_n(\hat{\xi}_{\tau_n})$, $\varphi^{(1)}(\xi_{\tau_n})$ by $\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n})$ and $\varphi^{(2)}(\xi_{\tau_n})$ by $\hat{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n})$ (where in the latter, $\bar{\rho}$ and $\bar{b}$
are consistent estimators of $\rho$ and $b$ such that $(\bar{\rho} - \rho) \log(n) = o_P(1)$), where throughout $\alpha_n = F(\xi_{\tau_n})$. Then, elementwise, for any $J \ge 1$,

$$ \begin{aligned} \widehat{M}_n^E &= \widehat{M}_n^\varphi + O_P(1/\sqrt{n(1-\tau_n)}) \\ &= \begin{pmatrix} \frac{2\gamma}{1-2\gamma} & \frac{\gamma}{1-\gamma} \\ \frac{\gamma}{1-\gamma} & \frac{\gamma}{1-\gamma} \end{pmatrix} + O_P(1/\sqrt{n(1-\tau_n)}), \end{aligned} \quad (C.5) $$

$$ \widehat{M}_n^\xi = \begin{pmatrix} \frac{2\gamma^3}{1-2\gamma} & \frac{\gamma^3}{1-\gamma} \\ \frac{\gamma^3}{1-\gamma} & \frac{\gamma^3}{1-\gamma} \end{pmatrix} + O_P(1/\sqrt{n(1-\tau_n)}), \quad (C.6) $$

$$ \widehat{M}_n = \begin{pmatrix} \frac{1-\gamma}{\gamma(1-2\gamma)} & \frac{\gamma}{1-2\gamma} \\ -\frac{\gamma}{1-2\gamma} & \frac{2\gamma^3}{1-2\gamma} \end{pmatrix} + O_P(1/\sqrt{n(1-\tau_n)}) \quad (C.7) $$

$$ \text{and } \widehat{\mathfrak{M}}_n(J) = \begin{pmatrix} \frac{\gamma^3(1-\gamma)}{(1-2\gamma)} & \frac{\gamma^3}{1-2\gamma} \\ \frac{\gamma^3}{1-2\gamma} & \frac{2\gamma^3}{1-2\gamma} \end{pmatrix} + O_P(1/\sqrt{n(1-\tau_n)}). \quad (C.8) $$

Proof of Lemma C.3. Under the stated assumptions:

* $\hat{\gamma}_{\tau_n}^{\text{E,BR}}$ is a $\sqrt{n(1-\tau_n)}$-consistent estimator of $\gamma$, see Theorem 1 in Girard et al. (2022b).

* $\hat{\xi}_{\tau_n}$ is a $\sqrt{n(1-\tau_n)}$-relatively consistent estimator of $\xi_{\tau_n}$ in the sense that $\hat{\xi}_{\tau_n}/\xi_{\tau_n} - 1 = O_P(1/\sqrt{n(1-\tau_n)})$, see Theorem 2 in Daouia et al. (2018) or Theorem 1 in Daouia et al. (2020).

* $\bar{Y}_n$ is a $\sqrt{n}$-consistent estimator of $\mathbb{E}(Y)$ by the standard central limit theorem.

* $\hat{F}_n(\hat{\xi}_{\tau_n})$ and $\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n})$ are $\sqrt{n(1-\tau_n)}$-relatively consistent estimators of $\bar{F}(\xi_{\tau_n})$ and $\varphi^{(1)}(\xi_{\tau_n})$, see Lemma C.2.

* $\hat{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n})$ is a $\sqrt{n(1-\tau_n)}$-relatively consistent estimator of $\varphi^{(2)}(\xi_{\tau_n})$, because

$$ \begin{aligned} \hat{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n}) &= \frac{2\gamma^2}{(1-\gamma)(1-2\gamma)} \xi_{\tau_n}^2 \bar{F}(\xi_{\tau_n}) (1 + O_P(1/\sqrt{n(1-\tau_n)})) \\ &= \varphi^{(2)}(\xi_{\tau_n}) (1 + O_P(1/\sqrt{n(1-\tau_n)})), \end{aligned} $$

see Equation (C.2).

These convergences will be used without further mention throughout the proof. We start by noting that

$$ \frac{\hat{\xi}_{\tau_n} - \bar{Y}_n}{2\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n}) + \hat{\xi}_{\tau_n} - \bar{Y}_n} - 1 = -\frac{2\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n})}{2\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n}) + \hat{\xi}_{\tau_n} - \bar{Y}_n} = -2\frac{\varphi^{(1)}(\xi_{\tau_n})}{\xi_{\tau_n}}(1 + o_P(1)) $$
by Lemma C.2. By Equation (C.2) then,

$$ \frac{\hat{\xi}_{\tau_n} - \bar{Y}_n}{2\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n}) + \hat{\xi}_{\tau_n} - \bar{Y}_n} - 1 = O_P(\bar{F}(\hat{\xi}_{\tau_n})) = O_P(1 - \tau_n) = o_P(1/q_{\tau_n}) = o_P(1/\sqrt{n(1 - \tau_n)}) $$

because $q_{\tau_n}$ is asymptotically equivalent to a multiple of $(1 - \tau_n)^{-\gamma}$ (see the remark below Example 2.3.11 p.49 in de Haan and Ferreira (2006)), and $\gamma < 1$. Consequently $\hat{M}_{n,11}^E = \hat{M}_{n,11}^{\varphi}(1 + o_P(1/\sqrt{n(1 - \tau_n)}))$ and likewise $\hat{M}_{n,12}^E = \hat{M}_{n,12}^{\varphi}(1 + o_P(1/\sqrt{n(1 - \tau_n)}))$. Now, since $1 - \alpha_n = \bar{F}(\hat{\xi}_{\tau_n})$ is estimated by $\hat{F}_n(\hat{\xi}_{\tau_n})$,

$$ \hat{M}_{n,11}^{\varphi} = (1 - \tau_n) \left( \frac{\tilde{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n})}{[\tilde{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n})]^2} - 1 \right) \quad \text{and} \quad \hat{M}_{n,12}^{\varphi} = \hat{M}_{n,22}^{\varphi} = (1 - \tau_n) \frac{1 - \hat{F}_n(\hat{\xi}_{\tau_n})}{\hat{F}_n(\hat{\xi}_{\tau_n})} $$

Using Equation (C.2) and assumption $\sqrt{n(1 - \tau_n)}(1/q_{\tau_n} + A((1 - \tau_n)^{-1})) = O(1)$, this entails

$$ \hat{M}_{n,11}^{\varphi} = \frac{2\gamma}{1 - 2\gamma} + O_P(1/\sqrt{n(1 - \tau_n)}) \quad \text{and} \quad \hat{M}_{n,12}^{\varphi} = \hat{M}_{n,22}^{\varphi} = \frac{\gamma}{1 - \gamma} + O_P(1/\sqrt{n(1 - \tau_n)}) $$

This proves Equation (C.5). Then, similarly,

$$ \frac{\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n})(2\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n}) + \hat{\xi}_{\tau_n} - \bar{Y}_n)}{\hat{\xi}_{\tau_n}(\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n}) + \hat{F}_n(\hat{\xi}_{\tau_n})(\hat{\xi}_{\tau_n} - \bar{Y}_n))} = \gamma + O_P(1/\sqrt{n(1 - \tau_n)}) $$

because of Equation (C.2), and therefore

$$ \hat{M}_{n,11}^{\xi} = \frac{2\gamma^3}{1 - 2\gamma} + O_P(1/\sqrt{n(1 - \tau_n)}), \quad \hat{M}_{n,12}^{\xi} = \frac{\gamma^3}{1 - \gamma} + O_P(1/\sqrt{n(1 - \tau_n)}) \\ \text{and } \hat{M}_{n,22}^{\xi} = \frac{\gamma^3}{1 - \gamma} + O_P(1/\sqrt{n(1 - \tau_n)}) $$

This proves Equation (C.6). As a consequence

$$ \begin{aligned} \hat{M}_{n,11} &= \frac{1}{(\hat{\gamma}_{\tau_n}^{\text{E,BR}})^2} \left( \frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1 - \tau_n} \right)^2 [\hat{M}_{n,11}^{\xi} - 2\hat{M}_{n,12}^{\xi} + \hat{M}_{n,22}^{\xi}] \\ &= \frac{1 - \gamma}{\gamma(1 - 2\gamma)} + O_P(1/\sqrt{n(1 - \tau_n)}), \\ \hat{M}_{n,12} &= \frac{1}{\hat{\gamma}_{\tau_n}^{\text{E,BR}}} \frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1 - \tau_n} [\hat{M}_{n,12}^{\xi} - \hat{M}_{n,11}^{\xi}] = -\frac{\gamma}{1 - 2\gamma} + O_P(1/\sqrt{n(1 - \tau_n)}) \\ \text{and } \hat{M}_{n,22} &= \hat{M}_{n,11}^{\xi} = \frac{2\gamma^3}{1 - 2\gamma} + O_P(1/\sqrt{n(1 - \tau_n)}) \end{aligned} $$
which is Equation (C.7). Finally, $\mathfrak{M}_{n,22}(J) = \hat{M}_{n,22}$, and besides

$$ \begin{aligned} \mathfrak{M}_{n,11}(J) &= \left( 1 + \frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1 - \tau_n} \right)^{-4} \hat{M}_{n,11} \left( 1 + \mathbb{O}_{\mathbb{P}}\left( \frac{1}{n(1 - \tau_n)} \right) \right) \\ &= \frac{\gamma^3(1 - \gamma)}{1 - 2\gamma} + \mathbb{O}_{\mathbb{P}}\left( \frac{1}{\sqrt{n(1 - \tau_n)}} \right) \end{aligned} $$

$$ \begin{aligned} \text{and } \mathfrak{M}_{n,12}(J) &= - \left( 1 + \frac{\hat{F}_n(\hat{\xi}_{\tau_n})}{1 - \tau_n} \right)^{-2} \hat{M}_{n,12} \left( 1 + \mathbb{O}_{\mathbb{P}}\left( \frac{1}{n(1 - \tau_n)} \right) \right) \\ &= \frac{\gamma^3}{1 - 2\gamma} + \mathbb{O}_{\mathbb{P}}\left( \frac{1}{\sqrt{n(1 - \tau_n)}} \right). \end{aligned} $$

This proves Equation (C.8) and completes the proof. □

## C.2 Proofs of the main results

*Proof of Theorem 1.* Recall from the proof of Lemma C.3 that one has, under the stated assumptions:

* $\hat{\gamma}_{\tau_n}^{\text{E,BR}}$ is a $\sqrt{n(1 - \tau_n)}$-consistent estimator of $\gamma$.

* $\hat{\xi}_{\tau_n}$ is a $\sqrt{n(1 - \tau_n)}$-relatively consistent estimator of $\xi_{\tau_n}$.

* $\bar{Y}_n$ is a $\sqrt{n}$-consistent estimator of $\mathbb{E}(Y)$.

* $\hat{F}_n(\hat{\xi}_{\tau_n})$ and $\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n})$ are $\sqrt{n(1 - \tau_n)}$-relatively consistent estimators of $\bar{F}(\xi_{\tau_n})$ and $\varphi^{(1)}(\xi_{\tau_n})$.

* $\hat{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n})$ is a $\sqrt{n(1 - \tau_n)}$-relatively consistent estimator of $\varphi^{(2)}(\xi_{\tau_n})$.

We turn to the proof of the desired results.

(i) First of all

$$ \frac{\sqrt{n(1 - \tau_n)}}{\log\left(\frac{1 - \tau_n}{1 - \tau'_n}\right)} \log \left( \frac{\hat{\xi}_{\tau_n}^{\star, \text{BR}}}{\xi_{\tau'_n}} \right) \xrightarrow{d} \mathcal{N} \left( 0, \frac{\gamma^3(1 - \gamma)}{1 - 2\gamma} \right) $$

by Theorem 2(i) in Girard et al. (2022b). By Slutsky's lemma, it follows that it is sufficient to prove that $\hat{s}_n^2(J) \to \gamma^3(1 - \gamma)/(1 - 2\gamma)$ in probability. Now, from straightforward calculations,

$$ \partial_1 \hat{u}_n \left( \left( 1 + \hat{F}_n(\hat{\xi}_{\tau_n})/(1 - \tau_n) \right)^{-1}, 1 \right) \to 1 \text{ and } \partial_2 \hat{u}_n \left( \left( 1 + \hat{F}_n(\hat{\xi}_{\tau_n})/(1 - \tau_n) \right)^{-1}, 1 \right) \to 0 $$

in probability as $n \to \infty$, where $\hat{u}_n$ is obtained from the function $u_n$ by replacing $\gamma$ by $\hat{\gamma}_{\tau_n}^{\text{E,BR}}$, $\xi_{\tau_n}$ by $\hat{\xi}_{\tau_n}$, $\mathbb{E}(Y)$ by $\bar{Y}_n$, $\bar{F}(\xi_{\tau_n})$ by $\hat{F}_n(\hat{\xi}_{\tau_n})$, $\varphi^{(1)}(\xi_{\tau_n})$ by $\hat{\varphi}_n^{(1)}(\hat{\xi}_{\tau_n})$ and $\varphi^{(2)}(\xi_{\tau_n})$ by
$\tilde{\varphi}_n^{(2)}(\hat{\xi}_{\tau_n})$. Consequently, by Lemma C.3,

$$
\begin{aligned}
\mathfrak{M}_{n,11}^{\text{BR}}(J) &= (\partial_1 \hat{u}_n((1 + \hat{F}_n(\hat{\xi}_{\tau_n})/(1-\tau_n))^{-1}, 1))^2 \mathfrak{M}_{n,11}(J) \\
&\quad + (\partial_2 \hat{u}_n((1 + \hat{F}_n(\hat{\xi}_{\tau_n})/(1-\tau_n))^{-1}, 1))^2 \mathfrak{M}_{n,22}(J) \\
&\quad + 2\partial_1 \hat{u}_n((1 + \hat{F}_n(\hat{\xi}_{\tau_n})/(1-\tau_n))^{-1}, 1)\partial_2 \hat{u}_n((1 + \hat{F}_n(\hat{\xi}_{\tau_n})/(1-\tau_n))^{-1}, 1) \mathfrak{M}_{n,12}(J), \\
&\quad \to \frac{\gamma^3(1-\gamma)}{1-2\gamma},
\end{aligned}
$$

$$
\begin{aligned}
\mathfrak{M}_{n,12}^{\text{BR}}(J) &= \partial_1 \hat{u}_n((1 + \hat{F}_n(\hat{\xi}_{\tau_n})/(1 - \tau_n))^{-1}, 1) \mathfrak{M}_{n,12}(J) \\
&\quad + \partial_2 \hat{u}_n((1 + \hat{F}_n(\hat{\xi}_{\tau_n})/(1 - \tau_n))^{-1}, 1) \mathfrak{M}_{n,22}(J), \\
&\quad \to \frac{\gamma^3}{1 - 2\gamma}, \quad \text{and}
\end{aligned}
$$

$$
\mathfrak{M}_{n,22}^{\text{BR}}(J) \to \frac{2\gamma^3}{1 - 2\gamma}
$$

in probability as $n \to \infty$, with the same notational convention. With that convention
again,

$$
\partial_1 \hat{g}_n(\hat{\gamma}_{\tau_n}^{\text{E,BR}}, 1) \to 1 \text{ and } \partial_2 \hat{g}_n(\hat{\gamma}_{\tau_n}^{\text{E,BR}}, 1) \to 0
$$

in probability as $n \to \infty$, so that

$$
\begin{aligned}
\hat{s}_n^2(J) ={}& (\partial_1 \hat{g}_n(\hat{\gamma}_{\tau_n}^{\text{E,BR}}, 1))^2 \mathfrak{M}_{n,11}^{\text{BR}}(J) + 2\partial_1 \hat{g}_n(\hat{\gamma}_{\tau_n}^{\text{E,BR}}, 1) \partial_2 \hat{g}_n(\hat{\gamma}_{\tau_n}^{\text{E,BR}}, 1) \mathfrak{M}_{n,12}^{\text{BR}}(J) \\
& + (\partial_2 \hat{g}_n(\hat{\gamma}_{\tau_n}^{\text{E,BR}}, 1))^2 \mathfrak{M}_{n,22}^{\text{BR}}(J) \to \frac{\gamma^3(1-\gamma)}{1-2\gamma}
\end{aligned}
$$

in probability as $n \to \infty$, as required.

(ii) Observe that

$$
\exp \left( \pm \frac{\log\left(\frac{(1-\tau_n)}{(1-\tau'_n)}\right)}{\sqrt{n(1-\tau_n)}} \sqrt{\hat{s}_n^2(J)} \times z_{1-\alpha/2} \right) \\
- \exp \left( \pm \frac{\log\left(\frac{(1-\tau_n)}{(1-\tau'_n)}\right)}{\sqrt{n(1-\tau_n)}} \sqrt{\frac{(\hat{\gamma}_{\tau_n}^{\text{E,BR}})^3 (1-\hat{\gamma}_{\tau_n}^{\text{E,BR}})}{1-2\hat{\gamma}_{\tau_n}^{\text{E,BR}}}} \times z_{1-\alpha/2} \right)
$$

is, in probability, asymptotically equivalent to

$$
\pm \frac{\log((1 - \tau_n)/(1 - \tau'_n))}{\sqrt{n(1 - \tau_n)}} \left( \sqrt{\hat{s}_n^2(J)} - \sqrt{\frac{(\hat{\gamma}_{\tau_n}^{\text{E,BR}})^3 (1 - \hat{\gamma}_{\tau_n}^{\text{E,BR}})}{1 - 2\hat{\gamma}_{\tau_n}^{\text{E,BR}}}} \right) \times z_{1-\alpha/2},
$$

and therefore also to

$$
\pm \frac{1}{2} \sqrt{\frac{1-2\gamma}{\gamma^3(1-\gamma)}} \frac{\log\left(\frac{(1-\tau_n)}{(1-\tau'_n)}\right)}{\sqrt{n(1-\tau_n)}} \left( \hat{s}_n^2(J) - \frac{(\hat{\gamma}_{\tau_n}^{\text{E,BR}})^3 (1-\hat{\gamma}_{\tau_n}^{\text{E,BR}})}{1-2\hat{\gamma}_{\tau_n}^{\text{E,BR}}} \right) \times z_{1-\alpha/2}.
$$

This means that

$$
\begin{aligned}
& \text{length}(\hat{I}_{\tau_n'}^{(2)}(\alpha; J)) - \text{length}(\hat{I}_{\tau_n'}^{(1)}(\alpha)) \\
& = \hat{\xi}_{\tau_n'}^{\text{BR}} \sqrt{\frac{1-2\gamma}{\gamma^3(1-\gamma)}} \frac{\log\left(\frac{(1-\tau_n)}{(1-\tau_n')}\right)}{\sqrt{n(1-\tau_n)}} \left( \hat{s}_n^2(J) - \frac{(\hat{\gamma}_{\tau_n}^{\text{E,BR}})^3 (1-\hat{\gamma}_{\tau_n}^{\text{E,BR}})}{1-2\hat{\gamma}_{\tau_n}^{\text{E,BR}}} \right) z_{1-\alpha/2}(1+o_{\mathbb{P}}(1)).
\end{aligned}
$$
Now

$$
\begin{align*}
\hat{s}_n^2(J) - \frac{(\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}})^3 (1 - \hat{\gamma}_{\tau_n}^{\mathrm{E,BR}})}{1 - 2\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}} &= ((\partial_1 \hat{g}_n(\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}, 1))^2 - 1) \hat{\mathfrak{M}}_{n,11}^{\mathrm{BR}}(J) \\
&\quad + \hat{\mathfrak{M}}_{n,11}^{\mathrm{BR}}(J) - \frac{(\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}})^3 (1 - \hat{\gamma}_{\tau_n}^{\mathrm{E,BR}})}{1 - 2\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}} + 2\partial_1 \hat{g}_n(\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}, 1) \partial_2 \hat{g}_n(\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}, 1) \hat{\mathfrak{M}}_{n,12}^{\mathrm{BR}}(J) \\
&\quad + (\partial_2 \hat{g}_n(\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}, 1))^2 \hat{\mathfrak{M}}_{n,22}^{\mathrm{BR}}(J).
\end{align*}
$$

Equation (B.10) shows that $1 - \partial_1 \hat{g}_n(\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}, 1) = o_P(1/\log((1-\tau_n)/(1-\tau'_n)))$, and Equation (B.11) shows that $\partial_2 \hat{g}_n(\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}, 1)$ is asymptotically equivalent to $1/\log((1-\tau_n)/(1-\tau'_n))$ in probability. Therefore

$$
\begin{equation}
\begin{split}
\hat{s}_n^2(J) - \frac{(\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}})^3 (1 - \hat{\gamma}_{\tau_n}^{\mathrm{E,BR}})}{1 - 2\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}} ={}& \hat{\mathfrak{M}}_{n,11}^{\mathrm{BR}}(J) - \frac{\gamma^3(1-\gamma)}{1-2\gamma} \\
& + \frac{1}{\log((1-\tau_n)/(1-\tau'_n))} \left( \frac{2\gamma^3}{1-2\gamma} + o_P(1) \right).
\end{split}
\tag{C.9}
\end{equation}
$$

Besides

$$
\begin{align*}
& \hat{\mathfrak{M}}_{n,11}^{\mathrm{BR}}(J) - \frac{\gamma^3(1-\gamma)}{1-2\gamma} \\
&= ((\partial_1 \hat{u}_n((1+\hat{F}_n(\hat{\xi}_{\tau_n}))/(1-\tau_n))^{-1}, 1))^2 - 1) \hat{\mathfrak{M}}_{n,11}(J) + \hat{\mathfrak{M}}_{n,11}(J) \\
&\quad - \frac{\gamma^3(1-\gamma)}{1-2\gamma} \\
&\quad + (\partial_2 \hat{u}_n((1+\hat{F}_n(\hat{\xi}_{\tau_n}))/(1-\tau_n))^{-1}, 1))^2 \hat{\mathfrak{M}}_{n,22}(J) \\
&\quad + 2\partial_1 \hat{u}_n((1+\hat{F}_n(\hat{\xi}_{\tau_n}))/(1-\tau_n))^{-1}, 1) \partial_2 \hat{u}_n((1+\hat{F}_n(\hat{\xi}_{\tau_n}))/(1-\tau_n))^{-1}, 1) \hat{\mathfrak{M}}_{n,12}(J).
\end{align*}
$$

By Equations (B.8) and (B.9), $1 - \partial_1 \hat{u}_n((1+\hat{F}_n(\hat{\xi}_{\tau_n}))/(1-\tau_n))^{-1}, 1) = O_P(1/\xi_{\tau_n}) + O_P(1 - \tau_n) = O_P(1/q_{\tau_n}) = O_P(1/\sqrt{n(1-\tau_n)}) = o_P(1/\log((1-\tau_n)/(1-\tau'_n)))$, and similarly $\partial_2 \hat{u}_n((1+\hat{F}_n(\hat{\xi}_{\tau_n}))/(1-\tau_n))^{-1}, 1) = O_P(1/\xi_{\tau_n}) = o_P(1/\log((1-\tau_n)/(1-\tau'_n)))$. Applying Lemma C.3 thus yields

$$
\hat{\mathfrak{M}}_{n,11}^{\mathrm{BR}}(J) - \frac{\gamma^3(1-\gamma)}{1-2\gamma} = o_P(1/\log((1-\tau_n)/(1-\tau'_n))). \quad (C.10)
$$

Combining Equations (C.9) and (C.10) results in

$$
\hat{s}_n^2(J) - \frac{(\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}})^3 (1 - \hat{\gamma}_{\tau_n}^{\mathrm{E,BR}})}{1 - 2\hat{\gamma}_{\tau_n}^{\mathrm{E,BR}}} = \frac{1}{\log((1-\tau_n)/(1-\tau'_n))} \left( \frac{2\gamma^3}{1-2\gamma} + o_P(1) \right).
$$

From this one can conclude that

$$
\begin{align*}
\text{length}(\hat{I}_{\tau_n'}^{(2)}(\alpha; J)) - \text{length}(\hat{I}_{\tau_n'}^{(1)}(\alpha)) &= 2 \times \frac{\hat{\xi}_{\tau_n'}^{\star, \mathrm{BR}}}{\sqrt{n(1-\tau_n)}} \left( \sqrt{\frac{\gamma^3}{(1-\gamma)(1-2\gamma)}} z_{1-\alpha/2} + o_P(1) \right) \\
&= 2 \times \frac{\xi_{\tau_n'}}{\sqrt{n(1-\tau_n)}} \left( \sqrt{\frac{\gamma^3}{(1-\gamma)(1-2\gamma)}} z_{1-\alpha/2} + o_P(1) \right).
\end{align*}
$$

This completes the proof. $\square$
Proof of Theorem 2. The proof generally follows the ideas of the proof of Theorem 1 with the necessary adaptations. We discuss the proof in detail below.

(i) First of all

$$ \frac{\sqrt{k_n}}{\log(k_n/(n(1-\tau'_n)))} \log\left(\frac{\tilde{\xi}_{\tau'_n}^{\star,\mathrm{BR}}}{\xi_{\tau'_n}}\right) \xrightarrow{d} \mathcal{N}(0, \gamma^2) $$

by Theorem 2(ii) in Girard et al. (2022b). By Slutsky's lemma, it follows that it is sufficient to prove that $\tilde{\sigma}_n^2(J) \to \gamma^2$ in probability. To this end note that each term $v_{j,12}$ and $v_{j,22}$ appearing in the expressions of $V_{12}(J)$ and $V_{22}(J)$ is a continuous function of $\gamma \in (0,1)$, because $\phi: x \mapsto -x \log(x^{-1}-1)$ is infinitely differentiable on $(0,1)$. As a consequence, the estimators $\hat{V}_{12}(J)$ and $\hat{V}_{22}(J)$ appearing in $\tilde{\sigma}_n^2(J)$ converge in probability to $V_{12}(0) = V_{12} = \gamma^2 m(\gamma)$ and $V_{22}(0) = V_{22} = \gamma^2(1+(m(\gamma))^2)$, respectively. Moreover, $\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}$ is a $\sqrt{k_n}$-consistent estimator of $\gamma$ under the stated assumptions (as recalled in e.g. Girard et al. (2022b)), $\bar{Y}_n$ is a consistent estimator of $\mathbb{E}(Y)$ by the weak law of large numbers, and

$$ \log \left[ \frac{1}{\xi_{\tau'_n}} \left( \frac{n(1-\tau'_n)}{k_n} \right)^{-\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}} \left( 1/\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}} - 1 \right)^{-\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}} \hat{q}_{1-k_n/n} \right] \to 0 $$

in probability, using the asymptotic proportionality relationship $\xi_{\tau'_n} \sim (\gamma^{-1}-1)^{-\gamma} q_{\tau'_n}$ and the consistency of the Weissman extreme quantile estimator. The estimators $\partial_1 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}, 1)$ and $\partial_2 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}, 1)$ in $\tilde{\sigma}_n^2(J)$ thus converge in probability to 1 and 0, respectively, by (B.16) and (B.17). Conclude that $\tilde{\sigma}_n^2(J) \to \gamma^2$ in probability, as required.

(ii) The quantity

$$ \exp \left( \pm \frac{\log(k_n/(n(1-\tau'_n)))}{\sqrt{k_n}} \sqrt{\tilde{\sigma}_n^2(J)} \times z_{1-\alpha/2} \right) - \exp \left( \pm \frac{\log(k_n/(n(1-\tau'_n)))}{\sqrt{k_n}} \hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}} \times z_{1-\alpha/2} \right) $$

is, in probability, asymptotically equivalent to

$$ \pm \frac{1}{2\gamma} \frac{\log(k_n/(n(1-\tau'_n)))}{\sqrt{k_n}} (\tilde{\sigma}_n^2(J) - (\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}})^2) \times z_{1-\alpha/2}. $$

Then clearly

$$ \begin{aligned} \tilde{\sigma}_n^2(J) - (\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}})^2 & = (\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}})^2 ((\partial_1 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}, 1))^2 - 1) + (\partial_2 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}, 1))^2 \hat{V}_{22}(J) \\ & \quad + 2\partial_1 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}, 1) \partial_2 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}, 1) \hat{V}_{12}(J) \end{aligned} $$

and it suffices to control this expression to obtain the asymptotic behavior of the difference in length of the two asymptotic confidence intervals. As in the proof of Theorem 1(ii), Equations (B.16) and (B.17) show that $1 - \partial_1 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}, 1) = o_P(1/\log(k_n/(n(1-\tau'_n))))$ and $\partial_2 \hat{h}_n(\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}}, 1)$ is asymptotically equivalent (in probability) to $1/\log(k_n/(n(1-\tau'_n)))$. As a consequence

$$ \tilde{\sigma}_n^2(J) - (\hat{\gamma}_{1-k_n/n}^{\mathrm{H,BR}})^2 = 2\gamma^2 m(\gamma)/\log(k_n/(n(1-\tau'_n))) + o_P(1/\log(k_n/(n(1-\tau'_n)))) $$
Conclude that

$$
\begin{align*}
\text{length}(\tilde{I}_{\tau_n'}^{(2)}(\alpha; J)) - \text{length}(\tilde{I}_{\tau_n'}^{(1)}(\alpha)) &= 2 \times \frac{\tilde{\xi}_{\tau_n'}^{\star, \text{BR}}}{\sqrt{k_n}} (\gamma m(\gamma) z_{1-\alpha/2} + o_{\mathbb{P}}(1)) \\
&= 2 \times \frac{\tilde{\xi}_{\tau_n'}^{\star}}{\sqrt{k_n}} (\gamma m(\gamma) z_{1-\alpha/2} + o_{\mathbb{P}}(1)).
\end{align*}
$$

The proof is complete.

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

D Further finite-sample illustrations and results

D.1 Simulated data

This section contains a complete set of results in our 18 models, for both of our inference
methods (LAWS and quantile-based) and our three sample sizes n = 200, 1,000, 5,000.

D.2 Real data

This section contains extra results about our cyber insurance real data analysis. We
provide estimates of the tail index (Figure D.29), extreme quantile (Figure D.30) and
extreme expectile (Figure D.31) in each cluster formed by stratifying upon the type of
cyber breach.