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

HAL
open science

# Approximate Bayesian Computation of reduced-bias extreme risk measures from heavy-tailed distributions

Jonathan El Methni, Stéphane Girard

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

## To cite this version:

Jonathan El Methni, Stéphane Girard. Approximate Bayesian Computation of reduced-bias extreme risk measures from heavy-tailed distributions. 2025. hal-04965629v2

HAL Id: hal-04965629

https://inria.hal.science/hal-04965629v2

Preprint submitted on 6 Oct 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Approximate Bayesian Computation of reduced-bias
extreme risk measures from heavy-tailed distributions

Jonathan El Methni¹† and Stéphane Girard¹*†

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

*Corresponding author(s). E-mail(s): stephane.girard@inria.fr;

Contributing authors: jonathan.el-methni@univ-grenoble-alpes.fr;

†These authors contributed equally to this work.

Abstract

Most of extrapolation methods dedicated to the estimation of extreme risk measures rely on
the approximation of the excesses distribution above a high threshold by a Generalized Pareto
Distribution (GPD). We propose an alternative to the GPD, called the Refined Pareto Distri-
bution (RPD), which allows for a second-order approximation of the excesses distribution. The
parameters of the RPD are estimated using an Approximate Bayesian Computation (ABC)
method, and reduced-bias estimators of extreme risk measures are then derived together with
the associated credible intervals. The ABC estimator demonstrates good performance over a
wide range of heavy-tailed distributions. Its usefulness is also illustrated on two data sets of
insurance claims.

**Keywords:** Extreme-value statistics, Generalized Pareto Distribution, Heavy-tailed distributions, Risk measures.

# 1 Introduction

Assessing the extreme behavior of a random phenomenon is a major issue in environmental science, quantitative finance, and insurance. As an example, extreme weather events usually have strong negative consequences such as loss of life, infrastructure damage, as well as decrease of agricultural production, see Rootzén and Tajvidi (1997) for the impacts of wind storms in terms of insurance. When the phenomenon of interest is modeled by a quantitative random variable, the risk may be assessed by the associated quantile of level 1 − α ∈ (0, 1), referred to as the Value at Risk (VaR) in finance or insurance and to as the return level in environmental sciences. Focusing on extreme risks, when α is small, the quantile of interest is likely to be larger than the greatest observation, extrapolation methods are then necessary to estimate such an extreme quantile.

Peaks Over Threshold (POT) and block maxima methods are the two main statistical approaches for conducting inference in distributions tails. The first one is based on the approximation of the excesses distribution above a high threshold by a Generalized Pareto Distribution (GPD) while the second one relies on the approximation of the maxima distribution by an extreme-value distribution. Here, we focus on the POT method (see for instance (Beirlant et al., 2006, Section 5.3) for an overview) whose principle may be described as follows.

Let $X$ be the random variable of interest with cumulative distribution function $F$. The associated tail quantile function is defined as $t \ge 1 \mapsto U(t) = F^{-1}(1 - 1/t)$ with $F^{-1}(\cdot) = \inf\{x \in$
$\mathbb{R}$, $F(x) > \cdot$} the generalized inverse of $F$. Let $u \ge 1$ be a threshold, and consider $Y = X/u$ given $X \ge u$ the excess of $X$ above $u$. Then, clearly, the tail quantile function of $Y$ is $t \ge 1 \mapsto U(t/\alpha)/u$ with $\alpha = \mathbb{P}(X \ge u)$ and, following the Pickands-Balkema-de Haan theorem (Pickands, 1975, Theorem 7), it converges, as $u$ tends to the endpoint $U(\infty)$, to the tail quantile function of the GPD given for all $t \ge 1$ by

$$U_{\text{GPD}}(t | \gamma, \sigma, \mu) = \mu + \frac{\sigma}{\gamma} \int_{1}^{t} x^{\gamma-1} dx =: \mu + \sigma K_{\gamma}(t), \qquad (1)$$

where $\mu \in \mathbb{R}$, $\sigma > 0$ and $\gamma \in \mathbb{R}$ are respectively the location, scale and shape parameters of the GPD. Note that $\gamma$ coincides with the extreme-value index i.e. the shape parameter of the extreme-value distribution. Consequently, the tail quantile function of $X$ can be approximated by

$$U(t) \approx u U_{\text{GPD}}(\alpha t | \gamma, \sigma, \mu = 1), \qquad (2)$$

when $t$ is large, which provides a theoretical basis for extrapolating in the distribution tails. Most estimators of extreme risk measures take advantage of this property: See (de Haan and Ferreira, 2006, Section 4.3) for quantiles, Daouia et al. (2019) for extremiles, Daouia et al. (2018) for expec- tiles and El Methni and Stupfler (2017) for Wang risk measures. This approach however suffers from some limitations due to the fact that, in some cases, the GPD approximation may be too crude. First, the asymptotic approximation (2) of $U(t)$ only holds for $t$ large enough, and as such, it can- not be used to model the bulk of the distribution. Second, Equation (2) can be seen as a first-order approximation of the tail behavior and may thus be oversimplified for some distributions.

To overcome these vexing effects, a first line of works proposed extensions of the GPD including additional parameters to get more flexibility. For instance, Beirlant et al. (2009); Papastathopoulos and Tawn (2013) introduce an additional shape parameter in the GPD while keeping the tail behavior unaffected. Naveau et al. (2016) and Tencalic et al. (2020) build on these works to model the full range of $X$ (assumed to be a positive random variable). Both the lower and the upper tails are distributed with extended GPDs while allowing a smooth transition in-between for the bulk, see also Stein (2021). Similarly, the parametric GPD extensions proposed in Gamet and Jalbert (2022) provide better estimate for low thresholds $u$ than the usual GPD model. Semiparametric models are introduced in Carreau and Bengio (2009); Carreau et al. (2009); Carreau and Vrac (2011) by stitching a GPD with a mixture of Gaussian densities. Another working direction consists in refining (2) by considering second-order expansions of the excess tail quantile function $U(t/\alpha)/u$ as $t \to \infty$ based on the so-called second-order condition, see (de Haan and Ferreira, 2006, Section 2.3). We refer to Allouche et al. (2023); Gomes and Pestana (2007) for an application to the reduced bias estimation of quantiles, see also Allouche et al. (2024) and Daouia et al. (2024a) for conditional tail expectation and expectile frameworks.

Our approach lies at the intersection of these two lines of work. We exhibit a distribution for which the second-order expansion of the excess tail quantile function (induced by the second-order condition) is exact. This new distribution is referred to as the Refined Pareto Distribution (RPD) and is shown to be the natural extension of the GPD to model excesses. It enables the construction of extrapolation procedures with reduced bias compared to the extreme-value techniques based on the first-order GPD approximation (2). The RPD is defined through its tail quantile function $U_{\text{RPD}}$ which cannot be inverted in an explicit way. Its cumulative and density distribution functions are thus not closed form, making an inference based on the likelihood not tractable. In contrast, simu- lating data from the RPD is straightforward via $U_{\text{RPD}}$ and the inverse transform sampling method. Inference can then be performed in an efficient way using Approximate Bayesian Computation (ABC), see Karabatsos and Leisen (2018); Marin et al. (2012); Robert (2021) for more information on this topic. Numerical experiments conducted in a wide variety of heavy-tailed situations show that the proposed ABC estimator outperforms several competitors including classical ones such as GPD and Weissman estimators (Weissman, 1978) but also recent reduced bias estimators (Gomes and Pestana, 2007). Moreover, a nice side product of the method is the availability of empirical
credible intervals without extra computation, which makes it possible to quantify the estimation uncertainty.

The remaining of the paper is organized as follows. First, we introduce in Section 2 the Refined Pareto Distribution. Its main properties are highlighted in terms of approximation of the excesses distribution (proofs are available in the Supplementary material). The application to the estimation of extreme risk measures (specifically: quantiles and conditional tail moments) is described in Section 3 with a special focus on the ABC procedure. Next, the finite-sample properties of the estimators are illustrated in Section 4 on simulated data. We further investigate in Section 5 an application to two data sets of insurance claims. The respective code in R is publicly available on the github page of the first author.

# 2 Statistical framework

## 2.1 Extreme-value theory

Let {$X_1, \dots, X_n$} be an independent sample from a cumulative distribution function $F$ and let $X_{1,n} \le \dots \le X_{n,n}$ denote the associated order statistics. Let $U$ be the tail quantile function defined as $U(t) = F^{-1}(1 - 1/t)$ for all $t \ge 1$. We assume that $F$ belongs to the Fréchet maximum domain of attraction, which is equivalent to assuming that $U$ is regularly-varying with index $\gamma > 0$:

$$ \lim_{t \to \infty} \frac{U(ty)}{U(t)} = y^{\gamma}, \qquad (3) $$

for all $y > 0$. See Bingham et al. (1989) for a detailed account on regular variation. From (de Haan and Ferreira, 2006, Theorem 1.2.1), the first-order condition (3) is also equivalent to supposing that the survival function $(1-F)$ is regularly-varying with index $-1/\gamma$. This latter property shows that the associated distribution is heavy-tailed and emphasizes the crucial role of the parameter $\gamma$ referred to as the tail-index: The larger $\gamma$, the heavier the tail. Heavy-tailed distributions have been revealed useful to describe the tail structure of actuarial data (Resnick, 2007, Page 1) as well as extreme events in climatic risk studies (Kottegoda and Rosso, 2008, Section 7). Examples include the GPD, Burr, Fisher, Inverse gamma and Student $t$ distributions, see Section 4.

The following second-order condition refines the initial heavy-tail assumption by providing the rate of convergence in (3). The tail quantile function $U$ is assumed to be second-order regularly-varying with index $\gamma > 0$, second-order parameter $\rho < 0$ and auxiliary function $A$ having ultimately constant sign and converging to zero at infinity, i.e.

$$ \lim_{t \to \infty} \frac{1}{A(t)} \left( y^{-\gamma} \frac{U(ty)}{U(t)} - 1 \right) = K_{\rho}(y) = \frac{y^{\rho} - 1}{\rho}, \qquad (4) $$

for all $y > 0$. The auxiliary function $A$ drives the bias of most extreme-value estimators. In particular, the asymptotic sign of $A$ determines the sign of the asymptotic bias. Besides, necessarily, $|A|$ is regularly-varying with index $\rho$, see for instance (de Haan and Ferreira, 2006, Theorem 2.3.9). Since $A(t) \to 0$ as $t \to \infty$, the limit (4) can be equivalently rewritten as

$$ \lim_{t \to \infty} \frac{\log U(ty) - \log U(t) - \gamma \log(y)}{A(t)} = K_{\rho}(y), \qquad (5) $$

for all $y > 0$. We shall suppose that the auxiliary function in (4) is given by

$$ A(t) = \beta \gamma t^{\rho}, \qquad (6) $$

with $\beta \neq 0$. This situation (6) arises for instance in the Hall-Welsh class of heavy-tailed distributions (Hall, 1982) defined for all $t \ge 1$ by

$$ U(t) = C t^{\gamma} \left(1 + \beta \gamma t^{\rho} / \rho + o(t^{\rho})\right), \quad \text{with } C > 0, \qquad (7) $$
see (Daouia et al., 2024a, Table C2) for examples of classical distributions belonging to this class. Clearly, the Pareto distribution defined by $U(t) = U(1)t^γ$ verifies the first-order condition (3) in a non-asymptotic way, i.e. for all $t ≥ 1$ and $y > 0$. The goal of the next paragraph is to derive the analogous distribution for the second-order condition (5).

## 2.2 Refined Pareto Distribution

Let $R^* = \{x \in R, x \neq 0\}$. Our first result provides the set of tail functions $U$ such that the limit as $t \to \infty$ in (5) holds in fact for all $t \ge 1$.

**Proposition 1** Let $U$ be a differentiable tail quantile function such that

$$ \frac{1}{A(t)} (\log U(ty) - \log U(t) - \gamma \log y) = K_{\rho}(y), \quad (8) $$

for all $t \ge 1$ and $y > 0$, where $A(t) = \beta\gamma t^\rho$ with $(\gamma, \rho, \beta) \in (0, +\infty) \times (-\infty, 0) \times R^*$. Then, necessarily, for all $t \ge 1$, $U(t) = U(t | \gamma, \rho, \beta) = U(1) t^\gamma \exp(\beta\gamma K_\rho(t))$, which is a proper tail quantile function if $\beta \ge -1$.

The distribution associated with the above tail quantile function is referred to as the Refined Pareto Distribution.

**Definition 1** The Refined Pareto Distribution RPD($\gamma, \rho, \beta$) is defined by the tail quantile function

$$ U_{\text{RPD}}(t | \gamma, \rho, \beta) = U_{\text{RPD}}(1) t^{\gamma} \exp(\beta\gamma K_{\rho}(t)), \quad t \ge 1, $$

with $(\gamma, \rho, \beta) \in (0, +\infty) \times (-\infty, 0) \times [-1, \infty)$.

**Remark 1** Asymptotic expansions.

(i) $U_{\text{RPD}}(t | \gamma, \rho, \beta) = U_{\text{RPD}}(1) \exp(-\beta\gamma/\rho) t^\gamma (1 + \beta\gamma t^\rho/\rho + o(t^\rho))$ as $t \to \infty$. From (3), $U_{\text{RPD}}$ is thus regularly-varying with index $\gamma > 0$, which in turn proves that the RPD belongs to the Fréchet maximum domain of attraction. Moreover, comparing with (7), it appears that the RPD belongs to Hall-Welsh class with second-order parameter $\rho$. In contrast, in the GPD case, necessarily $\rho = -\gamma$ (Daouia et al., 2024a, Table C2).

(ii) $U_{\text{RPD}}(t | \gamma, \rho, \beta) \to U_{\text{RPD}}(1) t^\gamma$ as $\rho \to -\infty$, i.e. the RPD can be approximated by a standard Pareto distribution when the second-order parameter $\rho$ is small.

(iii) $U_{\text{RPD}}(t | \gamma, -\gamma, \beta) = U_{\text{RPD}}(1) ((1 + \beta)t^\gamma - \beta) + O(\beta^2)$ as $\beta \to 0$, i.e. the RPD can be approximated by a GPD with location parameter $\mu = U_{\text{RPD}}(1)$, scale parameter $\sigma = U_{\text{RPD}}(1)\gamma(1 + \beta)$ and shape parameter $\gamma$ when $\beta$ is close to zero, see (1).

In view of the above remarks, the RPD can thus be considered as the natural extension to the Pareto distribution taking account the second-order terms in the maximum domain attraction of Fréchet. The next Proposition provides an additional argument.

**Proposition 2** Let $X$ be a random variable with tail quantile function $U$ such that (5) and (6) hold. Define the excess above the threshold $u = U(1/\alpha)$ as $Y = X/u$ when $X \ge u$ with $\alpha \in (0, 1)$, and denote by $U_\alpha$ the associated tail quantile function.

(i) If $U_{\text{RPD}}(1) = 1$, then, for all $t \ge 1$,

$$ \lim_{\alpha \to 0} \frac{U_{\alpha}(t)}{U_{\text{RPD}}(t | \gamma, \rho, \beta\alpha^{-\rho})} = 1 + o(\alpha^{-\rho}). $$

(ii) If, moreover, $X \sim \text{RPD}(\gamma, \rho, \beta)$, then $Y \sim \text{RPD}(\gamma, \rho, \beta\alpha^{-\rho})$ for all $\alpha \in (0, 1)$.

Note that the condition $U_{\text{RPD}}(1) = 1$ ensures that the starting point of the RPD coincides with the one of the excess distribution. It thus appears that, similarly to the GPD, the RPD can be used to approximate the distribution of the excess above a high threshold: First, just like the
GPD, the RPD is stable with respect to thresholding. Up to our knowledge, these are the only two distributions sharing this property. Second, the RPD can approximate any excess distribution with a controlled error. The latter approximation rate is upper bounded by $\alpha^{-\rho}$ which naturally arises from the second-order condition. Third, in the RPD case, the parameters $\gamma$ and $\rho$ are independent (while in the GPD case $\rho = -\gamma$) which allows greater flexibility in modeling excesses: One can deal with the difficult situation where the excesses have both large $\gamma$ and $\rho$. Let us now provide an explicit derivation of two families of moments associated with the RPD that will reveal useful in the sequel.

## 2.3 Moments

First, a slightly different version from Hosking and Wallis (1987) of the Probability weighted moment (PWM) is adopted. Considering the moments of $\log Y$ instead of $Y$ allows an explicit calculation in the case of the RPD:

**Proposition 3** Let $Y \sim \text{RPD}(\gamma, \rho, \beta)$ with $(\gamma, \rho, \beta) \in (0, +\infty) \times (-\infty, 0) \times [-1, \infty)$ such that $U_{\text{RPD}}(1) = 1$. Consider the PWM defined by $m_Y(a) = (a+1)^2 \mathbb{E}[(\log Y)(\bar{F}_{\text{RPD}}(Y))^a]$, with $a \in \mathbb{R}$ and where $\bar{F}_{\text{RPD}}$ denotes the survival function associated with $U_{\text{RPD}}$. Then, $m_Y(a)$ is finite for all $a > -1$ and is given by

$$m_Y(a) = \gamma \left( 1 + \frac{\beta(a+1)}{a+1-\rho} \right).$$

Let us then introduce for all $x \neq y \in (-1, \infty)$ and $\rho \in (-\infty, 0)$,

$$r_{\rho}(x, y) = \frac{m_Y(y) - m_Y(0)}{m_Y(x) - m_Y(0)} = \frac{y x + 1 - \rho}{x y + 1 - \rho},$$

from Proposition 3, so that all three parameters can be expressed as functions of the PWMs:

$$\gamma = m_Y(0) + \frac{(y-x)(m_Y(y) - m_Y(0))}{y(x+1) - r_{\rho}(x, y)x(y+1)}, \qquad (9)$$

$$\rho = y + 1 + \frac{y(y - x)}{r_{\rho}(x, y)x - y}, \qquad (10)$$

$$\beta = \frac{xy(1 - r_{\rho}(x, y))}{r_{\rho}(x, y)x - y} \times \frac{(x - y)(m_Y(y) - m_Y(0))}{m_Y(y)(y - x) + x(y + 1)(1 - r_{\rho}(x, y))m_Y(0)}. \qquad (11)$$

This remark will allow us to build PWM estimators of $(\gamma, \rho, \beta)$ in Paragraph 3.4.

Second, the Conditional tail moment (CTM) introduced in El Methni et al. (2014) permits to recover various risk measures such as the Expected Shortfall or the conditional tail variance (Valdez, 2005), see (El Methni and Stupfler, 2017, Table 2) for additional examples.

**Proposition 4** Let $Y \sim \text{RPD}(\gamma, \rho, \beta)$ with $(\gamma, \rho, \beta) \in (0, +\infty) \times (-\infty, 0) \times [-1, \infty)$ such that $U_{\text{RPD}}(1) = 1$. Consider for $a \in \mathbb{R}$ and $\alpha \in (0, 1)$ the Conditional tail moment defined by $\text{CTM}_Y(a, \alpha) = \mathbb{E}(Y^a | Y \ge U_{\text{RPD}}(1/\alpha))$. Then, $\text{CTM}_Y(a, \alpha)$ is finite for all $a\gamma < 1$ and is given by

$$\text{CTM}_Y(a, \alpha | \gamma, \rho, \beta) = \exp(c) \alpha^{-a\gamma} \int_0^1 u^{-a\gamma} \exp(-c\alpha^{-\rho}u^{-\rho}) du,$$

with $c = -a\beta\gamma/\rho$.

The next section illustrates how the RPD can be used to design reduced-bias estimators of extreme risk measures such as extreme quantiles or conditional tail moments.
# 3 Reduced-bias estimation of extreme risk measures

## 3.1 Reduced-bias extrapolation

We begin by recalling the well-known link between an extreme and an intermediate quantile.

**Lemma 1** Let $k = k_n$ be an intermediate sequence, i.e. such that $k \to \infty$ and $k/n \to 0$ as $n \to \infty$. Consider $(p_n)$ an extreme probability sequence such that $p_n \in (0, 1)$ and $np_n \to 0$ as $n \to \infty$. Then, $\log U(1/p_n) = \log U(n/k) + \log U_{k/n}(d_n)$, where $d_n := k/(np_n)$ is the so-called extrapolation factor and $U_{k/n}(\cdot)$ is defined in Proposition 2.

In view of Proposition 2, under assumptions (5) and (6), the excesses $Y_{i,k} := X_{n-k+i,n}/X_{n-k,n}$ defined for $i \in \{1, \dots, k\}$ are approximately distributed from a RPD$(\gamma, \rho, \beta_n)$ with $\beta_n = \beta(n/k)^\rho$. Combining this result with Lemma 1, the extreme quantile can thus be approximated through

$$ \overline{\log U}(1/p_n) := \log U(n/k) + \log U_{\text{RPD}}(d_n | \gamma, \rho, \beta_n). $$

The extreme quantile is then estimated by replacing the intermediate quantile $U(n/k)$ by its empirical counterpart $X_{n-k,n}$ and plugging an estimator of $(\gamma, \rho, \beta_n)$ in the above approximation:

$$ \log \hat{U}(1/p_n) = \log X_{n-k,n} + \log U_{\text{RPD}}(d_n | \hat{\gamma}, \hat{\rho}, \hat{\beta}_n). $$

Here, by construction, $U_{\text{RPD}}(1) = 1$ and the estimator can be equivalently rewritten as

$$ \log \hat{U}(1/p_n) = \log X_{n-k,n} + \hat{\gamma} \log(d_n) + \hat{\gamma} \hat{\beta}_n K_{\hat{\rho}}(d_n). \quad (12) $$

One thus finds back the extreme quantile estimator introduced in (Gomes and Pestana, 2007, Equations (3.11) and (3.12)). The RPD framework induces a new interpretation of this estimator together with new estimators $(\hat{\gamma}, \hat{\rho}, \hat{\beta}_n)$ based on the $k$-sample of excesses, as illustrated in the two paragraphs below. The next result establishes that the proposed extrapolation method does not introduce any first-order bias in the estimation.

**Theorem 5** Let $X$ be a random variable with tail quantile function $U$ such that (5) and (6) hold. Let $k = k_n$ be an intermediate sequence such that $\sqrt{k}A(n/k) \to \lambda \in \mathbb{R}$ as $n \to \infty$. Let $(p_n)$ a probability sequence such that $d_n := k/(np_n) \to \infty$ as $n \to \infty$. Let $(\hat{\gamma}, \hat{\rho}, \hat{\beta}_n)$ be some estimators of $(\gamma, \rho, \beta_n := \beta(n/k)^\rho)$ such that $\sqrt{k}(\hat{\gamma} - \gamma) \xrightarrow{d} \xi$, $\hat{\beta}_n/\beta_n = O_P(1)$ and $\hat{\rho} \le \rho^\dagger$ almost surely for some $\rho^\dagger < 0$. Then,

$$ \frac{\sqrt{k}}{\log(d_n)} \left( \log \hat{U}(1/p_n) - \log U(1/p_n) \right) \xrightarrow{d} \xi. $$

It appears that the extreme quantile estimator based on the RPD inherits its asymptotic distribution from $\hat{\gamma}$ without any additional bias term even when $\sqrt{k}A(n/k) \to \lambda \neq 0$. Moreover, this result does not require convergent estimators of the second-order parameters. One can then build an asymptotically unbiased extreme quantile estimator from an estimator of $\gamma$ benefiting from the same property. A similar extrapolation principle can be derived for the CTM.

**Lemma 2** Under the assumptions of Lemma 1 and, if moreover $a\gamma < 1$, then

$$ \log \text{CTM}_X(a, p_n) = a \log U(n/k) + \log \text{CTM}_{k/n}(a, 1/d_n), $$

where $\text{CTM}_{k/n}$ is the conditional tail moment associated with $U_{k/n}$.

Consequently, the extreme CTM can thus be approximated through

$$ \overline{\log \text{CTM}}_X(a, p_n) := a \log U(n/k) + \log \text{CTM}_Y(a, 1/d_n | \gamma, \rho, \beta_n). $$
where $CTMy$ refers to the conditional tail moment associated with the RPD and given in Proposition 4. As previously, the extreme CTM is estimated by replacing the intermediate quantile $U(n/k)$ by its empirical counterpart $X_{n-k,n}$ and by plugging an estimator of $(\gamma, \rho, \beta_n)$ in the above approximation:

$$ \begin{aligned} \log \widehat{CTM}_X(a, p_n) = & a \log X_{n-k,n} \\ & + \log \mathrm{CTM}_Y(a, 1/d_n | \hat{\gamma}, \hat{\rho}, \hat{\beta}_n) \\ = & a \log X_{n-k,n} + \hat{c} + a\hat{\gamma} \log d_n \\ & + \log \int_0^1 u^{-a\hat{\gamma}} \exp(-\hat{c}d_n^{\hat{\rho}}u^{-\hat{\rho}})du, \end{aligned} \quad (13) $$

from Proposition 4 with $\hat{c} = -a\hat{\beta}_n\hat{\gamma}/\hat{\rho}$. Both estimators (12) and (13) are based on preliminary estimators of $(\gamma, \rho, \beta_n)$. The computation of maximum likelihood estimators seems to be out of reach since the density function of the RPD is not available. In contrast, sampling from the RPD is easy and fast because the quantile function is explicit. We thus propose hereafter Bayesian estimators that do not use the likelihood but rather a sampling procedure.

## 3.2 Bayesian estimators

Let $\theta = (\gamma, \rho, \beta_n)$ be the parameter of interest and suppose we are given a prior distribution $\pi(\theta)$ and a sample $\mathcal{Y} := \{Y_1, \dots, Y_k\}$ from RPD($\theta$). The goal of Bayesian inference is to compute $p(\theta | \mathcal{Y})$, the posterior distribution of $\theta$ given $\mathcal{Y}$. When the likelihood $\mathcal{L}_{\theta}(\mathcal{Y})$ is available, Bayes' Theorem states that $p(\theta | \mathcal{Y})$ is proportional to $\pi(\theta)\mathcal{L}_{\theta}(\mathcal{Y})$. Bayesian estimators of $\theta$ as well as credible intervals can then be derived, at least from a theoretical point of view. Here, we limit ourselves to assuming that one is able to sample from $p(\theta | \mathcal{Y})$ to get $M$ realizations $\{\tilde{\theta}_1, \dots, \tilde{\theta}_M\}$, with $\tilde{\theta}_j = (\tilde{\gamma}_j, \tilde{\rho}_j, \tilde{\beta}_{j,n})$, for any $j \in \{1, \dots, M\}$. The empirical posterior mean of $\log U(1/p_n)$ is deduced from (12):

$$ \log \hat{U}_{\text{Bayes}}(1/p_n) = \log X_{n-k,n} + \hat{\gamma}_{\text{Bayes}} \log d_n + \frac{1}{M} \sum_{j=1}^{M} (\tilde{\gamma}_j \tilde{\beta}_{j,n} K_{\tilde{\rho}_j}(d_n)), \quad (14) $$

where $\hat{\gamma}_{\text{Bayes}}$ is the empirical posterior mean of $\gamma$ defined as

$$ \hat{\gamma}_{\text{Bayes}} = \frac{1}{M} \sum_{j=1}^{M} \tilde{\gamma}_j. \quad (15) $$

Note that one could define Bayesian estimators of $\rho$ and $\beta_n$ in a similar way. Moreover, credible intervals are easily obtained by considering the empirical quantiles of the ABC outputs $\log \tilde{U}_j(1/p_n) = \log X_{n-k,n} + \tilde{\gamma}_j \log d_n + \tilde{\gamma}_j \tilde{\beta}_{j,n} K_{\tilde{\rho}_j}(d_n)$, $j \in \{1, \dots, M\}$. We provide in the next paragraph a simple strategy to sample from $p(\theta | \mathcal{Y})$ using an Approximate Bayesian computation (ABC) algorithm, see Karabatsos and Leisen (2018); Marin et al. (2012) for reviews.

## 3.3 Approximation Bayesian computation

ABC is a “likelihood-free” inference method: It enables to perform approximate Bayesian analysis on any statistical model from which it is possible to sample new data, without the need to explicitly evaluate the likelihood function. Let $\mathcal{Y}^{\text{obs}} = \{Y_1^{\text{obs}}, \dots, Y_k^{\text{obs}}\}$ be the observed data and assume that the statistical model for $\mathcal{Y}^{\text{obs}}$ is described by a likelihood function $\mathcal{L}_{\theta}$ with parameter $\theta$ which is to be inferred. The basic scheme of one step of ABC is described in Algorithm 1.

The outcome of repeated uses of Algorithm 1 is a sample of $\theta$ values approximately distributed according to its posterior distribution. The basic (rejection) ABC approach in Step 3 amounts to specifying a tolerance level $\varepsilon > 0$, and then keeping $\theta$ if $d(\mathcal{Y}, \mathcal{Y}^{\text{obs}}) < \varepsilon$ for some distance $d(\cdot, \cdot)$
between samples. Here, we adopt the common approach which consists in selecting the tolerance level $\epsilon$ as a fixed quantile of the distance $d(\mathcal{Y}, \mathcal{Y}^{\text{obs}})$. More specifically, Algorithm 1 is repeated $M'$ times, out of which the $M$ smallest distances $d(\mathcal{Y}, \mathcal{Y}^{\text{obs}})$ are retained, yielding a quantile of order $M/M'$ (Robert, 2021). In the following, the focus is on the 1-Wasserstein distance (Villani, 2008, Chapter 6) denoted by $d_W(\cdot, \cdot)$. It is a distance function defined between probability distributions which can be easily computed on the (univariate) log-excesses, as described in Algorithm 2.

Finally, we introduce PWM estimators based on the results of Proposition 3 that will serve as easy-to-calculate competitors.

## 3.4 PWM estimators

The PWM associated with the excesses {$Y_1, \dots, Y_k$} is estimated by its empirical counterpart: For any $a > -1$,

$$ 
\begin{aligned}
\hat{m}_Y(a) &= \frac{(a+1)^2}{k} \sum_{i=1}^{k} \left(1 - \frac{i-1}{k}\right)^a \log Y_{i,k} \\
&= \frac{(a+1)^2}{k} \sum_{i=1}^{k} \left(\frac{i}{k}\right)^a \log(X_{n-i+1,n}/X_{n-k,n}).
\end{aligned}
\quad (16) $$

Remark that $\hat{m}_Y(a)$ can be interpreted as a weighted Hill estimator of the extreme-value index while $\hat{m}_Y(0)$ is the Hill estimator (Hill, 1975). In particular, $\hat{m}_Y(a) \xrightarrow{P} \gamma$ for all $a > -1$ under the first-order condition (3), see (Hüsler et al., 2006, Theorem 2.1). Replacing in (9)-(11) yields the following family of estimators:

$$ 
\left\{
\begin{aligned}
\hat{\gamma}_{\text{PWM}}(x, y) &= \hat{m}_Y(0) + \frac{(y-x)(\hat{m}_Y(y) - \hat{m}_Y(0))}{y(x+1) - \hat{r}(x, y)x(y+1)}, \\
\hat{\rho}_{\text{PWM}}(x, y) &= y + 1 + \frac{y(y-x)}{\hat{r}(x, y)x-y}, \\
\hat{\beta}_{\text{PWM}}(x, y) &= \frac{xy(1-\hat{r}(x, y))}{\hat{r}(x, y)x-y} \times \frac{(x-y)(\hat{m}_Y(y) - \hat{m}_Y(0))}{\hat{m}_Y(y)(y-x)+x(y+1)(1-\hat{r}(x, y))\hat{m}_Y(0)},
\end{aligned}
\right.
\quad (17) $$

with $\hat{r}(x, y) = (\hat{m}_Y(y) - \hat{m}_Y(0))/(\hat{m}_Y(x) - \hat{m}_Y(0))$, which is defined for all $x \neq y \in (-1, \infty)$. Let us note that $\hat{\gamma}_{\text{PWM}}(y, x) = \hat{\gamma}_{\text{PWM}}(x, y)$, $\hat{\rho}_{\text{PWM}}(y, x) = \hat{\rho}_{\text{PWM}}(x, y)$ and $\hat{\beta}_{\text{PWM}}(y, x) = \hat{\beta}_{\text{PWM}}(x, y)$. The finite sample behavior of the proposed ABC procedure and PWM estimators is assessed on simulated data in the next section.

# 4 Validation on simulated data

## 4.1 Experimental design

The efficiency of the ABC estimator is illustrated on seven heavy-tailed distributions. The first one is the RPD (see Definition 1) with $U_{\text{RPD}}(1) = \beta = 1$. The next six distributions: Burr, Fréchet, Fisher, GPD, Inverse Gamma, and Student $t$ all belong to the Hall-Welsh class (7), they satisfy the second-order condition (4) with (6), see (Daouia et al., 2024a, Table C2) for their definitions and associated values of $\gamma$ and $\rho$. For all seven distributions four tail-index values $\gamma \in \{1/8, 1/4, 1/2, 1\}$ are investigated. The second-order parameter $\rho$ depends on the considered distribution:

* Fréchet distribution: the second-order parameter is fixed to $\rho = -1$.

* Fisher, GPD and Inverse Gamma: the second-order parameter is fixed to $\rho = -\gamma$. In the case of the Fisher distribution, we set $\nu_1 = 3$.

* Student $t$ distribution: the second-order parameter is fixed to $\rho = -2\gamma$.

* For both RPD and Burr distributions, the second-order parameter can be chosen independently from the tail-index, five values are tested $\rho \in \{-1/8, -1/4, -1/2, -1, -2\}$.
Summarizing, this experimental design includes 60 configurations. In each case, $N = 500$ replications of a data set of size $n = 500$ are simulated from the 60 considered parametric models. Finally, we focus on the order of the extreme quantile given by $p_n = 1/(2n) = 10^{-3}$.

## 4.2 Competitors

The first considered estimator is the Weissman estimator Weissman (1978) given by

$$ \log \hat{U}_W(1/p_n) = \log X_{n-k,n} + \hat{\gamma}_{\text{Hill}} \log(d_n), \quad (18) $$

where $\hat{\gamma}_{\text{Hill}}(k)$ is the Hill estimator Hill (1975) defined as the mean of the log-excesses:

$$ \hat{\gamma}_{\text{Hill}}(k) = \hat{m}_Y(0) = \frac{1}{k} \sum_{i=1}^{k} \log Y_{i,k}, \quad (19) $$

see (16). Weissman estimator can be seen as a simplified version of (12) with $\hat{\beta}_n = 0$ which amounts to neglecting the second-order term in the Hall-Welsh model (7).

The second estimator is the so-called Corrected Weissman estimator (CW) introduced in Gomes and Pestana (2007):

$$ \log \hat{U}_{\text{CW}}(1/p_n) = \log X_{n-k,n} + \hat{\gamma}_{\text{CH}}(k) \log d_n + \hat{\gamma}_{\text{CH}}(k) \hat{\beta}_n K_{\hat{\rho}}(d_n), \quad (20) $$

where $\hat{\gamma}_{\text{CH}}(k)$ is the Corrected-Hill estimator (Caeiro et al., 2005):

$$ \hat{\gamma}_{\text{CH}}(k) = \hat{\gamma}_{\text{Hill}}(k) \left( 1 - \frac{\hat{\beta}_n}{1 - \hat{\rho}} \right), \quad (21) $$

with $\hat{\rho}$ and $\hat{\beta}_n$ some estimators of the second-order parameters $\rho$ and $\beta_n = \beta(n/k)^\rho$ which are not detailed here. Both estimators $\hat{\gamma}_{\text{CH}}$ and $\log \hat{U}_{\text{CW}}$ are based on a precise asymptotic analysis of the bias terms involved in $\hat{\gamma}_{\text{Hill}}$ and $\log \hat{U}_W$ respectively. The proposed corrections provide asymptotically unbiased estimators when $\sqrt{k}A(n/k) \to \lambda \neq 0$ leading to improved finite sample performance, as illustrated in Gomes and Pestana (2007). Clearly, (12) and (20) share the same structure, even if they were obtained by different ways.

Third, the proposed Bayesian estimator (14) is combined with the ABC procedure described in Algorithm 2 to sample $M = 100$ data from the RPD using $M' = 2,000$ trials in the ABC procedure, leading to $M/M' = 5\%$. The prior parameters are set to $\rho_{\text{loc}} = \hat{\rho}$, $\beta_{\text{loc}} = \hat{\beta}_n$ (the estimators used in (20)) and $\gamma_{\text{max}} = \max\{\hat{\gamma}_{\text{Hill}}(2), \dots, \hat{\gamma}_{\text{Hill}}(|n/2|)\}$, where $\cdot$ denotes the lower integer part.

Fourth, we consider the PWM estimator obtained by plugging in (12) the estimators $\hat{\gamma}_{\text{PWM}}(1, 2)$, $\hat{\rho}_{\text{PWM}}(1, 2)$ and $\hat{\beta}_{\text{PWM}}(1, 2)$ given in (17) i.e. using only the first two moments.

Finally, the classical GPD estimator is computed by plugging in the approximation (2) the PWM estimator of the parameters (Hosking and Wallis, 1987). To summarize, five extreme quantile estimators are compared in the following. They are denoted by Weissman, CW, ABC, PWM and GPD.

## 4.3 Results

The performance of the extreme quantile estimators is assessed on a log scale using the empirical bias defined for any $k \in \{2, \dots, n-1\}$ by

$$ \text{Bias}_k \left( \hat{U}(1/p_n) \right) = \frac{1}{N} \sum_{i=1}^{N} \log \left( \frac{\hat{U}_k^{(i)} (1/p_n)}{U(1/p_n)} \right), $$
and the Relative mean-square error

$$ \text{RMSE}_k(\hat{U}(1/p_n)) = \frac{1}{N} \sum_{i=1}^{N} \left( \log \left( \frac{\hat{U}_k^{(i)}(1/p_n)}{U(1/p_n)} \right) \right)^2, $$

where, for all $i \in \{1, \dots, N = 500\}$, $\hat{U}_k^{(i)}(1/p_n = 1000)$ is the estimator computed using the $k$ excesses from the $i$th replication of size $n = 500$. We also consider two “minimum” versions: on the one hand $M_1$RMSE is the sample version of the optimal asymptotic mean-square error (de Haan and Ferreira, 2006, pages 77–83):

$$ M_1\text{RMSE}(\hat{U}(1/p_n)) = \min_{k \in \{2, \dots, 0.9n\}} \text{RMSE}_k(\hat{U}(1/p_n)), $$

and on the other hand $M_2$RMSE computed on a neighborhood of the optimal sample fraction associated with each sample

$$ \hat{k}(i) = \underset{k \in \{2, \dots, 0.9n\}}{\operatorname{argmin}} \left( \log \left( \frac{\hat{U}_k^{(i)}(1/p_n)}{U(1/p_n)} \right) \right)^2, $$

and defined as

$$ M_2\text{RMSE}(\hat{U}(1/p_n)) = \frac{1}{N} \sum_{i=1}^{N} \frac{1}{0.1n+1} \sum_{k=\hat{k}(i)-0.05n}^{\hat{k}(i)+0.05n} \left( \log \left( \frac{\hat{U}_k^{(i)}(1/p_n)}{U(1/p_n)} \right) \right)^2. $$

The average on $(0.1n+1)$ values is introduced to account for the uncertainty on the optimal sample fraction in practice. The obtained $M_2$RMSEs are provided in Table 1 for the RPD, in Table 2 for the Burr distribution and in Table 3 for Fréchet, Fisher, GPD, Inverse Gamma, and Student t distributions.

First, one can note that, as expected, for all considered estimators, when $\gamma$ is fixed, the $M_2$RMSE decreases as $|\rho|$ increases. Conversely, when $\rho$ is fixed, the $M_2$RMSE increases as $\gamma$ increases. This pattern is clearly visible in Table 1 and Table 2, where $\gamma$ and $\rho$ can be chosen independently for both RPD and Burr distributions. Additionally, it is noteworthy that, in the two cases where $\gamma = 1$ and $\rho \ge -1/4$ for the Burr distribution, all estimators fail to accurately estimate the logarithm of the extreme quantile, in the sense that $M_2$RMSE $\ge 100\%$. These two cases are the most challenging scenarios where both $\gamma$ and $\rho$ are large leading to very heavy tails with a large deviation from the Pareto model. Second, it clearly appears that the proposed ABC estimator is the most accurate one, delivering the best result in 40 out of 58 cases. When it does not provide the top performance, it consistently ranks second. The ABC estimator outperforms all the considered competitors in every situation for the RPD and Fréchet distributions or for all remaining distributions with large tail-index ($\gamma \ge 1/2$). The GPD estimator yields the best results in 13 cases associated with small tail-indices ($\gamma \le 1/4$). The original Weissman estimator only achieves the best result four times, in the RPD case with $\rho = -2$, a scenario where the bias is low and therefore no bias reduction is necessary, making the first-order approximation precise enough. The CW estimator only ranks first in one case of the Student t distribution, while the PWM estimator does not achieve any top ranking. Similar results are obtained using $M_1$RMSE, see Table B1–B3 in the Supplementary material.

As an illustration, the behavior of $k \mapsto \text{Bias}_k(\hat{U}(1000))$ and $k \mapsto \text{RMSE}_k(\hat{U}(1000))$ associated with the five above mentioned estimators is depicted on the first two rows of Figure 2 (Inverse Gamma distribution with $\gamma = -\rho = 1/2$ and Student t distribution with $2\gamma = -\rho = 1$). See also Figure A1 (RPD and Burr distributions with $\gamma = -\rho = 1/2$) and Figure A2 (Fisher and GPD distributions with $\gamma = -\rho = 1/2$) in the Supplementary material. It is clearly seen that the ABC
estimator enjoys the best performance in terms of bias and RMSE with stable sample paths for all considered distributions, making the choice of *k* less crucial in practice.

Finally, the third row of Figure 2 and Figure A1, A2 displays the ABC estimator as a function of $k \in \{2, \dots, n-1 = 499\}$ and the associated 90% credible interval computed on a single replication. It appears that, for the six considered distributions, the credible intervals are stable with respect to $k \in \{2, \dots, 250\}$, i.e. in the first half of the sample, and contain the theoretical extreme quantile.

# 5 Illustration on two real data sets

Two data sets of insurance claims are considered: the Secura Belgian reinsurance data set of $n = 371$ automobile claims (Beirlant et al., 2006), and the OpenFEMA data set of $n = 555$ residential flood insurance claims (Daouia et al., 2024b). Both data sets are accessible through the R software in the packages CASdatasets (Dutang and Charpentier, 2024) and Expectrem (see https://github.com/AntoineUC/Expectrem) respectively. The associated histograms are displayed in Figure 1. For each data set, the tail-index as well as two risk measures are estimated: the quantile or Value-at-Risk at level $p_n = 1/n$ and the Expected Shortfall at level $p_n = 0.5\%$. These levels are chosen following Allouche et al. (2023); Daouia et al. (2024b); El Methni and Stupfler (2017) making it possible to compare the estimated quantile and Expected Shortfall with their empirical counterparts. Furthermore, the use of risk measures at level $p_n = 0.5\%$ is recommended by major regulators, including those in the European Union, United Kingdom, Australia, and Canada (see Daouia et al. (2024b) for details). In the following, the prior parameters for the ABC algorithm are set to $\rho_{\text{loc}} = \hat{\rho}$, $\beta_{\text{loc}} = \hat{\beta}_n$ and $\gamma_{\text{max}} = \max\{\hat{\gamma}_{\text{CH}}(2), \dots, \hat{\gamma}_{\text{CH}}([n/2])\}$.

<figure><img src="image_4.png" /><figcaption>**Fig. 1:** Histogram of the data sets. Left panel: Secura Belgian automobile claims. Right panel: OpenFEMA flood claims.</figcaption></figure>

## 5.1 Secura Belgian automobile claims

The Secura Belgian reinsurance data set was first introduced and analyzed in Beirlant et al. (2006). It includes $n = 371$ automobile claims from 1998 to 2001, each valued at 1,200,000€ or more, and adjusted for inflation. This data set serves as a classic example of heavy-tailed data and was investigated from the perspective of extreme risk measures in Allouche et al. (2023); El Methni
and Stupfler (2017). First, the estimation of the tail-index is displayed on the top left panel of Figure 3 as a function of k for ABC (15), PWM (17), Hill (19), CH (21) and GPD estimators. It appears that Hill and GPD estimators both suffer from a severe bias when k ≥ 100 while the PWM estimator is very unstable. Both ABC and CH estimators point towards γ ≃ 1 / 4, which is consistent with the results obtained by Allouche et al. (2023) and El Methni and Stupfler (2017).

Our second goal is to estimate the extreme quantile U ( n ) and to compare it to the sample maximum x$_{n,n}$ = 7 , 898 , 639 € . The behavior of k ↦→ log ˆ U ( n ) associated with the five competitors is represented on the middle left panel of Figure 3. All considered estimators share a similar behavior for k ∈ { 2 , . . . , 100 } providing consistent results. Moreover, ABC, CW, PWM and GPD estimators benefit from a nice stable behavior on the wider range k ∈ { 2 , . . . , 300 } . These pointwise results can be complemented by the 1 − α = 90% credible intervals (for ABC) and asymptotic confidence intervals (for Weissman and CW). In the latter cases, the asymptotic confidence interval can be written as

$$
CI_{1-\alpha,k}(U(1/p_n)) = \left[ \log \hat{U}(1/p_n) - \sigma_1(k) z_{\alpha/2} \frac{\log d_n}{\sqrt{k}} ; \log \hat{U}(1/p_n) + \sigma_2(k) z_{\alpha/2} \frac{\log d_n}{\sqrt{k}} \right],
$$

with z$_{α/}$$_{2}$ the quantile of level α/ 2 from the standard Gaussian distribution and

$$
\sigma_1(k) = \sigma_2(k) = \hat{\gamma}_{\text{Hill}}(k) \quad \text{for Weissman estimator,}
$$

$$
\sigma_1(k) = \frac{\hat{\gamma}_{\text{CH}}(k)}{1 + \frac{z_{\alpha/2}}{\sqrt{k}} + \frac{\hat{\beta}(n/k)^{\hat{\rho}}}{1-\hat{\rho}}}
$$

$$
\sigma_2(k) = \frac{\hat{\gamma}_{\text{CH}}(k)}{1 - \frac{z_{\alpha/2}}{\sqrt{k}} + \frac{\hat{\beta}(n/k)^{\hat{\rho}}}{1-\hat{\rho}}}
$$

see Gomes and Pestana (2007) for details. The results are depicted on the right panel of Figure 3. It appears that the ABC credible intervals include the sample maximum for all k ∈ { 2 , . . . , 370 } which is hardly the case for the asymptotic confidence intervals associated with Weissman estimator indicating that the bias reduction is useful. The CW estimator also enjoys a nice behavior with the sample maximum belonging to the corresponding asymptotic confidence intervals for all k ∈ { 2 , . . . , 300 } even though its sample path is less regular than ABC’s. Thanks to its stability, the ABC estimate points towards a VaR around exp(15 . 9) ≃ 8 , 040 , 485 € which is coherent with the ones obtained in Allouche et al. (2023) with two reduced bias estimators, namely Refined Weissman: 8 , 298 , 000 € and CW: 8 , 203 , 000 € . It is worth noticing that ABC gives the closest result to the sample maximum.

Third, we investigate the estimation of the Expected Shortfall at level p$_{n}$ = 5% defined as ES$_{X}$ ( p$_{n}$ ) = CTM$_{X}$ (1 , p$_{n}$ ) = E ( X | X ≥ U (1 /p$_{n}$ )) . The ABC estimator derived from (13) is compared with the empirical version (Brazauskas et al., 2008; Deme et al., 2015) given by

$$
\widehat{\mathrm{ES}}_{\mathrm{EMP}}(p_n) = \frac{1}{[np_n]} \sum_{j=1}^{[np_n]} X_{n-j+1,n}. \qquad (22)
$$

The bottom left panel of Figure 3 shows the estimation of log ES$_{X}$ ( p$_{n}$ ) as a function of k ∈ { 2 , . . . , 370 } using the ABC procedure along with the associated 90% credible intervals. The behav- ior of the estimated value remains very stable across all thresholds with an estimate around exp(16) ≃ 8 , 886 , 111 € . Despite overestimating ˆ ES$_{EMP}$ (0 . 005) = 7 , 692 , 936 € , we obtain a value closer than the two estimates provided in El Methni and Stupfler (2017), namely 9 , 690 , 000 € and 9 , 854 , 000 € .
<caption>**Algorithm 1** General ABC algorithm: Simulation of one single sample</caption>

**Require:** $\mathcal{Y}^{\text{obs}} = \{Y_1^{\text{obs}}, \dots, Y_k^{\text{obs}}\}$
1: Sample $\theta$ from the prior distribution $\pi(\theta)$;
2: Given $\theta$, sample $\mathcal{Y} = \{Y_1, \dots, Y_k\}$ from $\mathcal{L}_{\theta}$;
3: **if** $\mathcal{Y}$ is too different from $\mathcal{Y}^{\text{obs}}$ **then**
4:     discard $\theta$,
5:     **else** keep $\theta$.
6: **end if**

<caption>**Algorithm 2** ABC algorithm for the RPD: Simulation of M samples</caption>

**Require:** $\mathcal{Y}^{\text{obs}} = \{Y_1^{\text{obs}}, \dots, Y_k^{\text{obs}}\}, \gamma_{\text{max}}, \rho_{\text{loc}}, \beta_{\text{loc}}, M', M$
1: **for** $\ell \in \{1, \dots, M'\}$ **do**
2:     a: Sample $\gamma^{(\ell)}$ from the uniform distribution $\mathcal{U}(0, \gamma_{\text{max}})$;
3:     b: Sample $\rho^{(\ell)}$ from the truncated normal distribution $\mathcal{T}\mathcal{N}(a_1, b_1, \mu_1, \sigma_1^2)$ with starting point $a_1 = -\infty$, endpoint $b_1 = 0$, location parameter $\mu_1 = \rho_{\text{loc}}$ and scale parameter $\sigma_1 = 1$, see (Johnson et al., 1995, Section 10.1) for details;
4:     c: Sample $\beta^{(\ell)}$ from $\mathcal{T}\mathcal{N}(a_2, b_2, \mu_2, \sigma_2^2)$ with starting point $a_2 = -1$, endpoint $b_2 = \infty$, location parameter $\mu_2 = \beta_{\text{loc}}$ and scale parameter $\sigma_2 = 1$;
5:     d: Given $(\gamma^{(\ell)}, \rho^{(\ell)}, \beta^{(\ell)})$, sample $\mathcal{Y}^{(\ell)} = \{Y_1^{(\ell)}, \dots, Y_k^{(\ell)}\}$ from RPD$(\gamma^{(\ell)}, \rho^{(\ell)}, \beta^{(\ell)})$;
6:     e: Compute the 1-Wasserstein distance on the logarithm of the excesses:

$$d_W(\mathcal{Y}^{(\ell)}, \mathcal{Y}^{\text{obs}}) = \frac{1}{k} \sum_{i=1}^{k} \left| \log(Y_{i,k}^{(\ell)} / Y_{i,k}^{\text{obs}}) \right|.$$

2: **end for**

3: The M smallest $d_W$ distances are recorded and the associated sets of parameters are denoted by $(\tilde{\gamma}_j, \tilde{\rho}_j, \tilde{\beta}_{j,n})$, $j \in \{1, \dots, M\}$.
## 5.2 OpenFEMA flood claims

The OpenFEMA data set includes $n = 555$ residential flood insurance claims in the USA in 2012, see Daouia et al. (2024b) for details on the construction of the data set. Following the same methodology as in the first real data application, we start by focusing on the estimation of the tail-index which is displayed on the top left panel of Figure 4 as a function of $k \in \{2, \dots, \lfloor n/2 \rfloor = 277\}$. Hill, CH and GPD estimators exhibit the classic behavior of tail-index estimators, with a large variance when $k$ is small, then a range of stability and finally a drift due to the bias when $k$ increases. The PWM estimator is unstable for low values of $k$ and stabilizes around $\gamma \approx 0.4$ for $k \in \{100, \dots, 150\}$. In contrast, the ABC estimator is stable throughout $k \in \{2, \dots, 200\}$ around $\gamma \approx 0.4$ which is in accordance with the results of Daouia et al. (2024b) based on the moment estimator (Dekkers et al., 1989).

Let us now estimate the extreme quantile $U(n)$ and compare it to the sample maximum $x_{n,n} = 634,262 \$. The behavior of $k \mapsto \log \hat{U}(n)$ associated with the five estimators of interest is represented on the middle left panel of Figure 4. All ABC, CW, PWM and GPD estimators share a similar behavior providing consistent values smaller than the sample maximum, while the Weissman estimator deviates as the bias increases. These results should be read in conjunction with the $1 - \alpha = 90\%$ credible intervals (for ABC) and asymptotic confidence intervals (for Weissman and CW) displayed on the right panel of Figure 4. Only the credible intervals associated with the ABC estimator contain the maximum over the entire range of $k$ values.

Finally, the estimation of $\text{ES}_X(0.005)$ with the ABC procedure enjoys a remarkable stability for all $k \in \{50, \dots, 200\}$ pointing towards $\exp(13) \approx 442,413 \$ which is very close to the empirical estimate $\hat{\text{ES}}_{\text{EMP}}(0.005) \approx 466,085 \$. Moreover, this estimate is included in all credibility intervals. Again these results are coherent with the estimation of $500,000 \$ provided in Daouia et al. (2024b) and also based on a reduced-bias approach.

# 6 Conclusion

We introduced the Refined Pareto Distribution (RPD) as a natural alternative to the GPD to approximate the excesses distribution above a high threshold. The RPD verifies the second-order condition in a non-asymptotic way, and as such, it allows the construction of reduced-bias approximations of extreme risk measures. The parameters of the RPD are estimated using a dedicated ABC method, leading to reduced-bias estimators of extreme risk measures (quantiles and conditional tail moments) and the associated credible intervals. The ABC estimator yields good results on several simulated heavy-tailed distributions and is moreover stable over a wide range of $k$ values, facilitating the practical choice of the sample fraction. A future avenue for research would be to extend this work dedicated to heavy-tailed distributions to all maximum domains of attraction.
<table><caption>**Table 1:** M₂RMSE (in %) associated with five estimators of log U(1/p_n = 1000) computed on N = 500 replications of a data set of size n = 500 from a RPD. The best result is emphasized in bold. M₂RMSEs larger than 100% are not reported.</caption><thead><tr><th>RPD</th><th>Weissman</th><th>ABC</th><th>CW</th><th>PWM</th><th>GPD</th></tr></thead><tbody><tr><td colspan="6"><strong>γ = 1/8</strong></td></tr><tr><td>ρ = -1/8</td><td>1.22</td><td><strong>0.83</strong></td><td>3.49</td><td>1.21</td><td>0.85</td></tr><tr><td>ρ = -1/4</td><td>1.08</td><td><strong>0.56</strong></td><td>3.71</td><td>1.02</td><td>0.61</td></tr><tr><td>ρ = -1/2</td><td>0.85</td><td><strong>0.37</strong></td><td>2.22</td><td>0.75</td><td>0.43</td></tr><tr><td>ρ = -1</td><td>0.45</td><td><strong>0.23</strong></td><td>0.84</td><td>0.75</td><td>0.49</td></tr><tr><td>ρ = -2</td><td><strong>0.24</strong></td><td>0.33</td><td>0.56</td><td>0.74</td><td>0.47</td></tr><tr><td colspan="6"><strong>γ = 1/4</strong></td></tr><tr><td>ρ = -1/8</td><td>4.90</td><td><strong>3.30</strong></td><td>13.96</td><td>4.86</td><td>4.31</td></tr><tr><td>ρ = -1/4</td><td>4.79</td><td><strong>2.31</strong></td><td>14.11</td><td>3.96</td><td>2.77</td></tr><tr><td>ρ = -1/2</td><td>3.39</td><td><strong>1.47</strong></td><td>8.89</td><td>2.99</td><td>1.86</td></tr><tr><td>ρ = -1</td><td>1.79</td><td><strong>0.92</strong></td><td>3.36</td><td>3.00</td><td>1.99</td></tr><tr><td>ρ = -2</td><td><strong>0.97</strong></td><td>1.30</td><td>2.23</td><td>2.95</td><td>1.91</td></tr><tr><td colspan="6"><strong>γ = 1/2</strong></td></tr><tr><td>ρ = -1/8</td><td>19.58</td><td><strong>13.22</strong></td><td>55.83</td><td>19.45</td><td>22.42</td></tr><tr><td>ρ = -1/4</td><td>19.14</td><td><strong>9.25</strong></td><td>56.44</td><td>15.85</td><td>12.55</td></tr><tr><td>ρ = -1/2</td><td>13.59</td><td><strong>5.89</strong></td><td>35.55</td><td>11.97</td><td>9.47</td></tr><tr><td>ρ = -1</td><td>7.15</td><td><strong>3.67</strong></td><td>13.45</td><td>12.01</td><td>10.79</td></tr><tr><td>ρ = -2</td><td><strong>3.89</strong></td><td>5.21</td><td>8.93</td><td>11.78</td><td>11.46</td></tr><tr><td colspan="6"><strong>γ = 1</strong></td></tr><tr><td>ρ = -1/8</td><td>78.32</td><td><strong>52.86</strong></td><td>-</td><td>77.80</td><td>-</td></tr><tr><td>ρ = -1/4</td><td>76.57</td><td><strong>36.99</strong></td><td>-</td><td>63.40</td><td>-</td></tr><tr><td>ρ = -1/2</td><td>54.38</td><td><strong>23.58</strong></td><td>-</td><td>47.88</td><td>86.75</td></tr><tr><td>ρ = -1</td><td>28.62</td><td><strong>14.69</strong></td><td>53.79</td><td>48.02</td><td>70.92</td></tr><tr><td>ρ = -2</td><td><strong>15.56</strong></td><td>20.86</td><td>35.73</td><td>47.13</td><td>71.40</td></tr></tbody></table>

<table><caption>**Table 2:** M₂RMSE (in %) associated with five estimators of log U(1/p_n = 1000) computed on N = 500 replications of a data set of size n = 500 from a Burr distribution. The best result is emphasized in bold. M₂RMSEs larger than 100% are not reported.</caption><thead><tr><th>Burr</th><th>Weissman</th><th>ABC</th><th>CW</th><th>PWM</th><th>GPD</th></tr></thead><tbody><tr><td colspan="6"><strong>γ = 1/8</strong></td></tr><tr><td>ρ = -1/8</td><td>7.82</td><td>5.16</td><td>6.26</td><td>5.25</td><td><strong>1.23</strong></td></tr><tr><td>ρ = -1/4</td><td>2.69</td><td>1.48</td><td>1.99</td><td>1.59</td><td><strong>1.02</strong></td></tr><tr><td>ρ = -1/2</td><td>1.11</td><td><strong>0.43</strong></td><td>0.57</td><td>0.58</td><td>0.77</td></tr><tr><td>ρ = -1</td><td>0.49</td><td><strong>0.13</strong></td><td>0.20</td><td>0.34</td><td>0.63</td></tr><tr><td>ρ = -2</td><td>0.25</td><td><strong>0.16</strong></td><td>0.39</td><td>0.56</td><td>0.48</td></tr><tr><td colspan="6"><strong>γ = 1/4</strong></td></tr><tr><td>ρ = -1/8</td><td>31.27</td><td>20.62</td><td>25.03</td><td>21.00</td><td><strong>4.21</strong></td></tr><tr><td>ρ = -1/4</td><td>10.74</td><td>5.92</td><td>7.96</td><td>6.36</td><td><strong>2.79</strong></td></tr><tr><td>ρ = -1/2</td><td>4.44</td><td><strong>1.71</strong></td><td>2.29</td><td>2.32</td><td>2.39</td></tr><tr><td>ρ = -1</td><td>1.99</td><td><strong>0.51</strong></td><td>0.78</td><td>1.35</td><td>2.25</td></tr><tr><td>ρ = -2</td><td>1.05</td><td><strong>0.55</strong></td><td>1.34</td><td>1.93</td><td>1.88</td></tr><tr><td colspan="6"><strong>γ = 1/2</strong></td></tr><tr><td>ρ = -1/8</td><td>-</td><td>82.50</td><td>-</td><td>83.99</td><td><strong>68.75</strong></td></tr><tr><td>ρ = -1/4</td><td>42.97</td><td>23.69</td><td>31.85</td><td>25.43</td><td><strong>7.51</strong></td></tr><tr><td>ρ = -1/2</td><td>17.77</td><td><strong>6.82</strong></td><td>9.17</td><td>9.27</td><td>10.79</td></tr><tr><td>ρ = -1</td><td>7.98</td><td><strong>2.06</strong></td><td>3.13</td><td>5.40</td><td>11.23</td></tr><tr><td>ρ = -2</td><td>4.01</td><td><strong>2.51</strong></td><td>6.19</td><td>8.87</td><td>10.62</td></tr><tr><td colspan="6"><strong>γ = 1</strong></td></tr><tr><td>ρ = -1/8</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>ρ = -1/4</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>ρ = -1/2</td><td>71.07</td><td><strong>27.29</strong></td><td>36.70</td><td>37.07</td><td>87.81</td></tr><tr><td>ρ = -1</td><td>31.94</td><td><strong>8.22</strong></td><td>12.51</td><td>21.61</td><td>71.20</td></tr><tr><td>ρ = -2</td><td>16.04</td><td><strong>10.04</strong></td><td>24.75</td><td>35.55</td><td>72.86</td></tr></tbody></table>
<table><caption>**Table 3:** M₂RMSE (in %) associated with five estimators of log U(1/p_n = 1000) computed on N = 500 replications of data sets of size n = 500 from five heavy-tailed distributions. The best result is emphasized in bold.</caption><thead><tr><th></th><th>Weissman</th><th>ABC</th><th>CW</th><th>PWM</th><th>GPD</th></tr></thead><tbody><tr><td colspan="6">Fréchet (ρ = -1)</td></tr><tr><td>γ = 1/8</td><td>0.33</td><td><strong>0.16</strong></td><td>0.37</td><td>0.51</td><td>0.50</td></tr><tr><td>γ = 1/4</td><td>1.33</td><td><strong>0.65</strong></td><td>1.50</td><td>2.04</td><td>1.90</td></tr><tr><td>γ = 1/2</td><td>5.30</td><td><strong>2.62</strong></td><td>6.00</td><td>8.15</td><td>10.94</td></tr><tr><td>γ = 1</td><td>21.21</td><td><strong>10.47</strong></td><td>24.00</td><td>32.60</td><td>71.89</td></tr><tr><td colspan="6">Fisher (ρ = -γ)</td></tr><tr><td>γ = 1/8</td><td>6.01</td><td>3.67</td><td>4.62</td><td>4.00</td><td><strong>1.06</strong></td></tr><tr><td>γ = 1/4</td><td>8.98</td><td>4.50</td><td>6.17</td><td>5.61</td><td><strong>2.41</strong></td></tr><tr><td>γ = 1/2</td><td>16.46</td><td><strong>6.34</strong></td><td>8.47</td><td>8.47</td><td>13.04</td></tr><tr><td>γ = 1</td><td>30.17</td><td><strong>8.43</strong></td><td>16.60</td><td>25.15</td><td>78.57</td></tr><tr><td colspan="6">GPD (ρ = -γ)</td></tr><tr><td>γ = 1/8</td><td>8.08</td><td>5.41</td><td>6.38</td><td>5.69</td><td><strong>1.18</strong></td></tr><tr><td>γ = 1/4</td><td>10.97</td><td>6.53</td><td>8.13</td><td>6.48</td><td><strong>3.09</strong></td></tr><tr><td>γ = 1/2</td><td>18.33</td><td><strong>7.13</strong></td><td>9.08</td><td>11.05</td><td>12.02</td></tr><tr><td>γ = 1</td><td>32.27</td><td><strong>8.02</strong></td><td>13.19</td><td>24.63</td><td>67.65</td></tr><tr><td colspan="6">Inverse Gamma (ρ = -γ)</td></tr><tr><td>γ = 1/8</td><td>2.19</td><td>0.98</td><td>1.50</td><td>1.23</td><td><strong>0.92</strong></td></tr><tr><td>γ = 1/4</td><td>4.21</td><td><strong>1.53</strong></td><td>2.70</td><td>2.37</td><td>2.07</td></tr><tr><td>γ = 1/2</td><td>9.33</td><td><strong>2.70</strong></td><td>5.41</td><td>7.10</td><td>13.42</td></tr><tr><td>γ = 1</td><td>21.65</td><td><strong>9.03</strong></td><td>21.65</td><td>29.64</td><td>71.30</td></tr><tr><td colspan="6">Student t (ρ = -2γ)</td></tr><tr><td>γ = 1/8</td><td>3.79</td><td>2.42</td><td>2.94</td><td>2.64</td><td><strong>1.09</strong></td></tr><tr><td>γ = 1/4</td><td>5.37</td><td>2.53</td><td>3.37</td><td>3.44</td><td><strong>2.35</strong></td></tr><tr><td>γ = 1/2</td><td>10.54</td><td>3.20</td><td><strong>2.75</strong></td><td>6.07</td><td>10.35</td></tr><tr><td>γ = 1</td><td>18.15</td><td><strong>6.11</strong></td><td>19.47</td><td>23.63</td><td>67.92</td></tr></tbody></table>
<figure><img src="image_5.png" /><figcaption>**Fig. 2:** Illustration on simulated data sets of size *n* = 500 from an Inverse Gamma distribution with *γ* = −ρ = 1/2 (left panel) and a Student *t* distribution with 2γ = −ρ = 1 (right panel). Top: *k* ∈ {2, ..., *n* − 1} ↦→ RMSE<sub>*k*</sub>(Û(1/*p*<sub>*n*</sub> = 1000)) and center: *k* ∈ {2, ..., *n* − 1} ↦→ Bias<sub>*k*</sub>(Û(1/*p*<sub>*n*</sub> = 1000)) computed on *N* = 500 replications associated with Weissman (blue), ABC (black), CW (green), PWM (purple) and GPD (orange) estimators. Bottom (in log scale): 90% credible intervals (blue) associated with the ABC estimator (black) computed on one replication. The theoretical extreme quantile *U*(1/*p*<sub>*n*</sub> = 1000) is depicted by a red horizontal line.</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Fig. 3:** Illustration on the Secura Belgian reinsurance data set. Top left panel: Estimation of the tail-index as a function of $k \in \{2, \dots, n-1 = 370\}$ by Hill (blue), ABC (black), CH (green), PWM (purple) and GPD (orange) estimators. Middle left panel: Estimation of the extreme quantile $\log U(n)$ as a function of $k \in \{2, \dots, n-1\}$ by Weissman (blue), ABC (black), CW (green), PWM (purple) and GPD (orange) estimators, $\log x_{n,n}$ is depicted by a red horizontal line. Bottom left panel: Estimation of $\log \text{ES}_X(p_n = 0.5\%)$ as a function of $k \in \{2, \dots, n-1\}$ by ABC (black) and the 90% credible intervals associated, $\log \hat{\text{ES}}_{\text{EMP}}(p_n)$ is depicted by a red horizontal line. Right panel: 90% credible intervals for $U(n)$ associated with the ABC estimator (top) and 90% asymptotic confidence intervals based on Weissman (middle) and CW (bottom) estimators, $\log x_{n,n}$ is depicted by a red horizontal line.</figcaption></figure>
<figure><img src="image_7.png" /><figcaption>**Fig. 4:** Illustration on the OpenFEMA data set. Top left panel: Estimation of the tail-index as a function of $k \in \{2, \dots, \lfloor n/2 \rfloor = 277\}$ by Hill (blue), ABC (black), CH (green), PWM (purple) and GPD (orange) estimators. Middle left panel: Estimation of the extreme quantile $\log U(n)$ as a function of $k \in \{2, \dots, 277\}$ by Weissman (blue), ABC (black), CW (green), PWM (purple) and GPD (orange) estimators, $\log x_{n,n}$ is depicted by a red horizontal line. Bottom left panel: Estimation of $\log \text{ES}_X(p_n = 0.5\%)$ as a function of $k \in \{2, \dots, 277\}$ by ABC (black) and the 90% credible intervals associated, $\log \hat{\text{ES}}_{\text{EMP}}(p_n)$ is depicted by a red horizontal line. Right panel: 90% credible intervals for $U(n)$ associated with the ABC estimator (top) and 90% asymptotic confidence intervals based on Weissman (middle) and CW (bottom) estimators, $\log x_{n,n}$ is depicted by a red horizontal line.</figcaption></figure>
References

Allouche, M., El Methni, J., and Girard, S. (2023). A refined Weissman estimator for extreme quantiles. *Extremes*, 26(3):545-572.

Allouche, M., El Methni, J., and Girard, S. (2024). Reduced-bias estimation of the extreme conditional tail expectation for Box-Cox transforms of heavy-tailed distributions. *Journal of Statistical Planning and Inference*, 233:106189.

Beirlant, J., Goegebeur, Y., Segers, J., and Teugels, J. L. (2006). *Statistics of extremes: theory and applications*. John Wiley & Sons.

Beirlant, J., Joossens, E., and Segers, J. (2009). Second-order refined peaks-over-threshold modelling for heavy-tailed distributions. *Journal of Statistical Planning and Inference*, 139(8):2800-2815.

Bingham, N. H., Goldie, C. M., and Teugels, J. L. (1989). *Regular variation*, volume 27. Cambridge university press.

Brazauskas, V., Jones, B. L., Puri, M. L., and Zitikis, R. (2008). Estimating conditional tail expectation with actuarial applications in view. *Journal of Statistical Planning and Inference*, 138(11):3590-3604.

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

Carreau, J. and Bengio, Y. (2009). A hybrid Pareto model for asymmetric fat-tailed data: the univariate case. *Extremes*, 12:53-76.

Carreau, J., Naveau, P., and Sauquet, E. (2009). A statistical rainfall-runoff mixture model with heavy-tailed components. *Water Resources Research*, 45(10).

Carreau, J. and Vrac, M. (2011). Stochastic downscaling of precipitation with neural network conditional mixture models. *Water Resources Research*, 47(10).

Daouia, A., Gijbels, I., and Stupfler, G. (2019). Extremiles: A new perspective on asymmetric least squares. *Journal of the American Statistical Association*, 114(527):1366-1381.

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., Stupfler, G., and Usseglio-Carleve, A. (2024a). Bias-reduced and variance-corrected asymptotic Gaussian inference about extreme expectiles. *Statistics and Computing*, 34(4):130.

Daouia, A., Stupfler, G., and Usseglio-Carleve, A. (2024b). A unified theory of extreme Expected Shortfall inference. Toulouse School of Economics Working Papers.

de Haan, L. and Ferreira, A. (2006). *Extreme value theory: an introduction*. Springer.

Dekkers, A. L., Einmahl, J. H., and de Haan, L. (1989). A moment estimator for the index of an extreme-value distribution. *The Annals of Statistics*, 17:1833-1855.

Deme, E. H., Girard, S., and Guillou, A. (2015). Reduced-bias estimator of the conditional tail expectation of heavy-tailed distributions. In *M. Hallin et al., editor, Mathematical Statistics and Limit Theorems: Festschrift in Honour of Paul Deheuvels*, pages 105-123. Springer.

Dutang, C. and Charpentier, A. (2024). *CASdatasets: Insurance datasets*. R package version 1.2-0, https://doi.org/10.57745/P0KHAG.

El Methni, J., Gardes, L., and Girard, S. (2014). Non-parametric estimation of extreme risk measures from conditional heavy-tailed distributions. *Scandinavian Journal of Statistics*, 41(4):988-1012.

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

Gamet, P. and Jalbert, J. (2022). A flexible extended generalized Pareto distribution for tail estimation. *Environmetrics*, 33(6):e2744.

Gomes, M. I. and Pestana, D. (2007). A sturdy reduced-bias extreme quantile (VaR) estimator. *Journal of the American Statistical Association*, 102(477):280-292.

Hall, P. (1982). On some simple estimates of an exponent of regular variation. *Journal of the Royal Statistical Society: Series B*, 44(1):37-42.

Hill, B. M. (1975). A simple general approach to inference about the tail of a distribution. *The Annals of Statistics*, 3:1163-1174.
Hosking, J. R. and Wallis, J. R. (1987). Parameter and quantile estimation for the generalized Pareto distribution. *Technometrics*, 29(3):339–349.

Hüsler, J., Li, D., and Müller, S. (2006). Weighted least squares estimation of the extreme value index. *Statistics & Probability Letters*, 76(9):920–930.

Johnson, N. L., Kotz, S., and Balakrishnan, N. (1995). *Continuous univariate distributions*, volume 2. John Wiley & sons.

Karabatsos, G. and Leisen, F. (2018). An approximate likelihood perspective on ABC methods. *Statistics Surveys*, 12:66–104.

Kottegoda, M. and Rosso, R. (2008). *Statistics, probability and reliability for civil and environmental engineers*. Mc-Graw-Hill Publishing Company, Colombus.

Marin, J.-M., Pudlo, P., Robert, C. P., and Ryder, R. J. (2012). Approximate Bayesian computational methods. *Statistics and Computing*, 22(6):1167–1180.

Naveau, P., Huser, R., Ribereau, P., and Hannart, A. (2016). Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection. *Water Resources Research*, 52(4):2753–2769.

Papastathopoulos, I. and Tawn, J. A. (2013). Extended generalised Pareto models for tail estimation. *Journal of Statistical Planning and Inference*, 143(1):131–143.

Pickands, J. (1975). Statistical inference using extreme order statistics. *The Annals of Statistics*, 3:119–131.

Resnick, S. (2007). *Heavy-Tail Phenomena: Probabilistic and Statistical Modeling*. Springer, New York.

Robert, C. (2021). Approximate bayesian computation, an introduction. In *Statistics for Astrophysics: Bayesian Methodology*, pages 77–112. EDP Sciences.

Rootzén, H. and Tajvidi, N. (1997). Extreme value statistics and wind storm losses: a case study. *Scandinavian Actuarial Journal*, 1:70–94.

Stein, M. L. (2021). A parametric model for distributions with flexible behavior in both tails. *Environmetrics*, 32(2):e2658.

Tencaliec, P., Favre, A.-C., Naveau, P., Prieur, C., and Nicolet, G. (2020). Flexible semiparametric generalized pareto modeling of the entire range of rainfall amount. *Environmetrics*, 31(2):e2582.

Valdez, E. A. (2005). Tail conditional variance for elliptically contoured distributions. *Belgian Actuarial Bulletin*, 5(1):26–36.

Villani, C. (2008). *Optimal transport: old and new*, volume 338. Springer.

Weissman, I. (1978). Estimation of parameters and large quantiles based on the $k$ largest observations. *Journal of the American Statistical Association*, 73(364):812–815.

**Supplementary information.** The supplementary material is a pdf file including the proofs of the main results and additional figures and tables to illustrate the experiments on simulated data.

**Acknowledgements.** This work is partially supported by the French National Research Agency (ANR) under the grant ANR-23-CE40-0009. S. Girard acknowledges the support of the Chair “Stress Test, Risk Management and Financial Steering”, led by the French Ecole Polytechnique and its Foundation and sponsored by BNP Paribas.

## Competing Interests

Authors have no financial or non-financial interests to disclose, directly or indirectly related to this work.
# Supplementary material for 'Approximate Bayesian Computation of reduced-bias extreme risk measures from heavy-tailed distributions'

Section A and Section B provide additional illustrations and numerical results associated with the experiments on simulated data (Section 4 of the main paper). The proofs of the theoretical results are given in Section C.

## A Additional figures

The behavior of $k \mapsto \text{Bias}_k(\hat{U}(1000))$ and $k \mapsto \text{RMSE}_k(\hat{U}(1000))$ associated with the five estimators described in Section 4.2 is depicted on the first two rows of Figure A1 (RPD and Burr distributions with $\gamma = -\rho = 1/2$) and Figure A2 (Fisher and GPD distributions with $\gamma = -\rho = 1/2$). The third row displays the ABC estimator as a function of $k \in \{2, \dots, n-1=499\}$ and the associated 90% credible interval computed on a single replication.
<figure><img src="image_8.png" /><figcaption>**Fig. A1:** Illustration on simulated data sets of size *n* = 500 from a RPD (left panel) and a Burr distribution (right panel) with $\gamma = -\rho = 1/2$ in both cases. Top: $k \in \{2, \dots, n-1\} \mapsto \text{RMSE}_k(\hat{U}(1/p_n = 1000))$ and center: $k \in \{2, \dots, n-1\} \mapsto \text{Bias}_k(\hat{U}(1/p_n = 1000))$ computed on *N* = 500 replications associated with Weissman (blue), ABC (black), CW (green), PWM (purple) and GPD (orange) estimators. Bottom (in log scale): 90% credible intervals (blue) associated with the ABC estimator (black) computed on one replication. The theoretical extreme quantile $\underline{U}(1/p_n = 1000)$ is depicted by a red horizontal line.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. A2:** Illustration on simulated data sets of size *n* = 500 from a Fisher distribution (left panel) and a GPD (right panel) with $\gamma = -\rho = 1/2$ in both cases. Top: $k \in \{2, \dots, n-1\} \mapsto \text{RMSE}_k(\hat{U}(1/p_n = 1000))$ and center: $k \in \{2, \dots, n-1\} \mapsto \text{Bias}_k(\hat{U}(1/p_n = 1000))$ computed on $N = 500$ replications associated with Weissman (blue), ABC (black), CW (green), PWM (purple) and GPD (orange) estimators. Bottom (in log scale): 90% credible intervals (blue) associated with the ABC estimator (black) computed on one replication. The theoretical extreme quantile $\underline{U}(1/p_n = 1000)$ is depicted by a red horizontal line.</figcaption></figure>
B Additional tables

The M₁RMSEs are provided in Table B1 for the RPD, in Table B2 for the Burr distribution and in Table B3 for Fréchet, Fisher, GPD, Inverse Gamma, and Student t distributions.
<table><caption>**Table B1:** M₁RMSE associated with five estimators of log U(1/p_n = 1000) computed on N = 500 replications of a data set of size n = 500 from a RPD. The best result is emphasized in bold. M₁RMSEs larger than 1 are not reported.</caption><thead><tr><th>RPD</th><th>Weissman</th><th>ABC</th><th>CW</th><th>PWM</th><th>GPD</th></tr></thead><tbody><tr><td colspan="6"><strong>γ = 1/8</strong></td></tr><tr><td>ρ = -1/8</td><td>0.0175</td><td><strong>0.0106</strong></td><td>0.0475</td><td>0.0202</td><td>0.0176</td></tr><tr><td>ρ = -1/4</td><td>0.0165</td><td><strong>0.0092</strong></td><td>0.0385</td><td>0.0198</td><td>0.0132</td></tr><tr><td>ρ = -1/2</td><td>0.0107</td><td><strong>0.0059</strong></td><td>0.0254</td><td>0.0140</td><td>0.0113</td></tr><tr><td>ρ = -1</td><td>0.0060</td><td><strong>0.0037</strong></td><td>0.0111</td><td>0.0115</td><td>0.0138</td></tr><tr><td>ρ = -2</td><td><strong>0.0035</strong></td><td>0.0042</td><td>0.0065</td><td>0.0073</td><td>0.0129</td></tr><tr><td colspan="6"><strong>γ = 1/4</strong></td></tr><tr><td>ρ = -1/8</td><td>0.0700</td><td><strong>0.0423</strong></td><td>0.1901</td><td>0.0810</td><td>0.0572</td></tr><tr><td>ρ = -1/4</td><td>0.0638</td><td><strong>0.0349</strong></td><td>0.1359</td><td>0.0740</td><td>0.0425</td></tr><tr><td>ρ = -1/2</td><td>0.0427</td><td><strong>0.0237</strong></td><td>0.1016</td><td>0.0558</td><td>0.0303</td></tr><tr><td>ρ = -1</td><td>0.0240</td><td><strong>0.0150</strong></td><td>0.0445</td><td>0.0459</td><td>0.0454</td></tr><tr><td>ρ = -2</td><td><strong>0.0139</strong></td><td>0.0168</td><td>0.0261</td><td>0.0293</td><td>0.0443</td></tr><tr><td colspan="6"><strong>γ = 1/2</strong></td></tr><tr><td>ρ = -1/8</td><td>0.2801</td><td><strong>0.1691</strong></td><td>0.7603</td><td>0.3239</td><td>0.2076</td></tr><tr><td>ρ = -1/4</td><td>0.2551</td><td>0.1397</td><td>0.5435</td><td>0.2961</td><td><strong>0.1203</strong></td></tr><tr><td>ρ = -1/2</td><td>0.1707</td><td><strong>0.0950</strong></td><td>0.4065</td><td>0.2234</td><td>0.1071</td></tr><tr><td>ρ = -1</td><td>0.0960</td><td><strong>0.0598</strong></td><td>0.1779</td><td>0.1835</td><td>0.1301</td></tr><tr><td>ρ = -2</td><td><strong>0.0555</strong></td><td>0.0671</td><td>0.1045</td><td>0.1171</td><td>0.1599</td></tr><tr><td colspan="6"><strong>γ = 1</strong></td></tr><tr><td>ρ = -1/8</td><td>-</td><td><strong>0.6764</strong></td><td>-</td><td>-</td><td>-</td></tr><tr><td>ρ = -1/4</td><td>-</td><td><strong>0.5588</strong></td><td>-</td><td>-</td><td>-</td></tr><tr><td>ρ = -1/2</td><td>0.6830</td><td><strong>0.3799</strong></td><td>-</td><td>0.8935</td><td>-</td></tr><tr><td>ρ = -1</td><td>0.3840</td><td><strong>0.2393</strong></td><td>0.7117</td><td>0.7340</td><td>0.8531</td></tr><tr><td>ρ = -2</td><td><strong>0.2219</strong></td><td>0.2685</td><td>0.4179</td><td>0.4683</td><td>0.8102</td></tr></tbody></table>

<table><caption>**Table B2:** M₁RMSE associated with five estimators of log U(1/p_n = 1000) computed on N = 500 replications of a data set of size n = 500 from a Burr distribution. The best result is emphasized in bold. M₁RMSEs larger than 1 are not reported.</caption><thead><tr><th>Burr</th><th>Weissman</th><th>ABC</th><th>CW</th><th>PWM</th><th>GPD</th></tr></thead><tbody><tr><td colspan="6"><strong>γ = 1/8</strong></td></tr><tr><td>ρ = -1/8</td><td>0.0593</td><td>0.0563</td><td>0.0585</td><td>0.0775</td><td><strong>0.0264</strong></td></tr><tr><td>ρ = -1/4</td><td>0.0268</td><td><strong>0.0182</strong></td><td>0.0249</td><td>0.0321</td><td>0.0246</td></tr><tr><td>ρ = -1/2</td><td>0.0134</td><td><strong>0.0065</strong></td><td>0.0081</td><td>0.0142</td><td>0.0173</td></tr><tr><td>ρ = -1</td><td>0.0062</td><td><strong>0.0029</strong></td><td>0.0049</td><td>0.0056</td><td>0.0152</td></tr><tr><td>ρ = -2</td><td>0.0035</td><td><strong>0.0027</strong></td><td>0.0055</td><td>0.0032</td><td>0.0132</td></tr><tr><td colspan="6"><strong>γ = 1/4</strong></td></tr><tr><td>ρ = -1/8</td><td>0.2373</td><td>0.2251</td><td>0.2339</td><td>0.3101</td><td><strong>0.0961</strong></td></tr><tr><td>ρ = -1/4</td><td>0.1071</td><td>0.0727</td><td>0.0995</td><td>0.1284</td><td><strong>0.0427</strong></td></tr><tr><td>ρ = -1/2</td><td>0.0537</td><td><strong>0.0259</strong></td><td>0.0324</td><td>0.0568</td><td>0.0566</td></tr><tr><td>ρ = -1</td><td>0.0249</td><td><strong>0.0115</strong></td><td>0.0194</td><td>0.0225</td><td>0.0520</td></tr><tr><td>ρ = -2</td><td>0.0146</td><td><strong>0.0103</strong></td><td>0.0198</td><td>0.0128</td><td>0.0456</td></tr><tr><td colspan="6"><strong>γ = 1/2</strong></td></tr><tr><td>ρ = -1/8</td><td>0.9493</td><td>0.9005</td><td>0.9355</td><td>-</td><td><strong>0.8022</strong></td></tr><tr><td>ρ = -1/4</td><td>0.4286</td><td>0.2908</td><td>0.3978</td><td>0.5136</td><td><strong>0.0626</strong></td></tr><tr><td>ρ = -1/2</td><td>0.2147</td><td><strong>0.1034</strong></td><td>0.1297</td><td>0.2272</td><td>0.1237</td></tr><tr><td>ρ = -1</td><td>0.0994</td><td><strong>0.0460</strong></td><td>0.0778</td><td>0.0900</td><td>0.1756</td></tr><tr><td>ρ = -2</td><td>0.0556</td><td><strong>0.0431</strong></td><td>0.0883</td><td>0.0513</td><td>0.1632</td></tr><tr><td colspan="6"><strong>γ = 1</strong></td></tr><tr><td>ρ = -1/8</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>ρ = -1/4</td><td>-</td><td>-</td><td>-</td><td>-</td><td>-</td></tr><tr><td>ρ = -1/2</td><td>0.8559</td><td><strong>0.4137</strong></td><td>0.5190</td><td>0.9090</td><td>-</td></tr><tr><td>ρ = -1</td><td>0.3978</td><td><strong>0.1841</strong></td><td>0.3112</td><td>0.3602</td><td>0.8302</td></tr><tr><td>ρ = -2</td><td>0.2263</td><td><strong>0.1723</strong></td><td>0.3522</td><td>0.2052</td><td>0.8257</td></tr></tbody></table>
<table><caption>**Table B3:** M₁RMSE associated with five estimators of log U(1/p_n = 1000) computed on N = 500 replications of data sets of size n = 500 from five heavy-tailed distributions. The best result is emphasized in bold.</caption><thead><tr><th></th><th>Weissman</th><th>ABC</th><th>CW</th><th>PWM</th><th>GPD</th></tr></thead><tbody><tr><td>Fréchet (ρ = -1)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>γ = 1/8</td><td>0.0047</td><td><strong>0.0031</strong></td><td>0.0047</td><td>0.0042</td><td>0.0135</td></tr><tr><td>γ = 1/4</td><td>0.0189</td><td><strong>0.0122</strong></td><td>0.0186</td><td>0.0169</td><td>0.0460</td></tr><tr><td>γ = 1/2</td><td>0.0759</td><td><strong>0.0488</strong></td><td>0.0745</td><td>0.0679</td><td>0.1624</td></tr><tr><td>γ = 1</td><td>0.3035</td><td><strong>0.1953</strong></td><td>0.2981</td><td>0.2716</td><td>0.8276</td></tr><tr><td>Fisher (ρ = -γ)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>γ = 1/8</td><td>0.0495</td><td>0.0412</td><td>0.0490</td><td>0.0683</td><td><strong>0.0329</strong></td></tr><tr><td>γ = 1/4</td><td>0.0860</td><td>0.0562</td><td>0.0766</td><td>0.1236</td><td><strong>0.0555</strong></td></tr><tr><td>γ = 1/2</td><td>0.2062</td><td><strong>0.1022</strong></td><td>0.1345</td><td>0.2158</td><td>0.1632</td></tr><tr><td>γ = 1</td><td>0.4578</td><td><strong>0.2051</strong></td><td>0.3587</td><td>0.3849</td><td>0.8800</td></tr><tr><td>GPD (ρ = -γ)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>γ = 1/8</td><td>0.0651</td><td>0.0551</td><td>0.0612</td><td>0.0785</td><td><strong>0.0253</strong></td></tr><tr><td>γ = 1/4</td><td>0.1088</td><td>0.0729</td><td>0.0977</td><td>0.1362</td><td><strong>0.0441</strong></td></tr><tr><td>γ = 1/2</td><td>0.2118</td><td><strong>0.1042</strong></td><td>0.1322</td><td>0.2253</td><td>0.1357</td></tr><tr><td>γ = 1</td><td>0.4030</td><td><strong>0.1823</strong></td><td>0.3030</td><td>0.3492</td><td>0.7971</td></tr><tr><td>Inverse Gamma (ρ = -γ)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>γ = 1/8</td><td>0.0256</td><td><strong>0.0146</strong></td><td>0.0234</td><td>0.0305</td><td>0.0260</td></tr><tr><td>γ = 1/4</td><td>0.0552</td><td><strong>0.0272</strong></td><td>0.0444</td><td>0.0579</td><td>0.0529</td></tr><tr><td>γ = 1/2</td><td>0.1268</td><td><strong>0.0605</strong></td><td>0.0843</td><td>0.1217</td><td>0.1755</td></tr><tr><td>γ = 1</td><td>0.3082</td><td><strong>0.1789</strong></td><td>0.2739</td><td>0.2529</td><td>0.8198</td></tr><tr><td>Student t (ρ = -2γ)</td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>γ = 1/8</td><td>0.0303</td><td><strong>0.0263</strong></td><td>0.0290</td><td>0.0427</td><td>0.0272</td></tr><tr><td>γ = 1/4</td><td>0.0572</td><td><strong>0.0343</strong></td><td>0.0477</td><td>0.0819</td><td>0.0541</td></tr><tr><td>γ = 1/2</td><td>0.1406</td><td>0.0619</td><td><strong>0.0501</strong></td><td>0.1170</td><td>0.1601</td></tr><tr><td>γ = 1</td><td>0.2721</td><td><strong>0.1579</strong></td><td>0.4192</td><td>0.2192</td><td>0.7659</td></tr></tbody></table>
C Proofs

Proof of Proposition 1.

Let us introduce $\varphi(x) = \log U(\exp x)$ for all $x \ge 0$. Replacing in (8) yields the equation

$$ \frac{1}{A(t)}(\varphi(\log t + \log y) - \varphi(\log t) - \gamma \log y) = K_{\rho}(y), $$

for all $t \ge 1$ and $y > 0$, or equivalently, letting $s = \log t$ and $x = \log y$,

$$ \frac{1}{A(\exp s)} \left( \frac{\varphi(s+x) - \varphi(s)}{x} - \gamma \right) = \frac{K_{\rho}(\exp x)}{x}, $$

for all $x \ne 0$ and $s \ge 0$. Letting $x \to 0$ and remarking that $K_{\rho}(\exp x)/x \to \rho$ yield

$$ \varphi'(s) = \gamma + \rho A(\exp s) = \gamma(1 + \beta \exp(\rho s)), $$

in view of (6). Integrating, it follows that, for all $s \ge 0$,

$$ \varphi(s) = \varphi(0) + \gamma(s + \beta K_{\rho}(\exp s)), $$

leading to

$$ U(x) = U(1)x^{\gamma} \exp(\beta\gamma K_{\rho}(x)), $$

for all $x \ge 1$. Conversely, it is easily checked that the above $U$ function is a solution of (8). Finally, the condition $\beta \ge -1$ is required to ensure that $U$ is increasing. ■

Proof of Proposition 2.

(i) The tail quantile function of $Y$ is given for all $y \ge 1$ by

$$ U_{\alpha}(y) = U(y/\alpha)/U(1/\alpha) \quad (23) $$

so that (5) can be rewritten with $t = 1/\alpha$ as

$$ \lim_{\alpha \to 0} \frac{1}{A(1/\alpha)} (\log U_{\alpha}(y) - \gamma \log(y)) = K_{\rho}(y), $$

or equivalently,

$$ \lim_{\alpha \to 0} \frac{1}{A(1/\alpha)} (\log U_{\alpha}(y) - \gamma \log(y) - A(1/\alpha)K_{\rho}(y)) = 0. $$

Taking account of (6), it follows

$$ \lim_{\alpha \to 0} \alpha^{\rho} (\log U_{\alpha}(y) - \gamma \log(y) - \beta \gamma \alpha^{-\rho} K_{\rho}(y)) = 0, $$

leading to

$$ \lim_{\alpha \to 0} \alpha^{\rho} (\log U_{\alpha}(y) - \log U_{\text{RPD}}(y | \gamma, \rho, \beta \alpha^{-\rho})) = 0, $$

if $U_{\text{RPD}}(1) = 1$. Taking the exponential concludes the proof.

(ii) The result is a consequence of (23) and of the identity

$$ \frac{U_{\text{RPD}}(y/\alpha | \gamma, \rho, \beta)}{U_{\text{RPD}}(1/\alpha | \gamma, \rho, \beta)} = (1/\alpha)^{\gamma} \exp \left( \frac{\beta\gamma}{\rho} \alpha^{-\rho} (y^{\rho} - 1) \right) = U_{\text{RPD}}(y | \gamma, \rho, \beta \alpha^{-\rho}) $$

that holds for all $y \ge 1$. ■
### Proof of Proposition 3.

The probability weighted moment of order $a > -1$ is given by

$$m_Y(a) = (a+1)^2 \mathbb{E}((\log Y)(\bar{F}_{\text{RPD}}(Y))^a) = (a+1)^2 \int_{U_{\text{RPD}}(1)}^{+\infty} \log(x)(\bar{F}_{\text{RPD}}(x))^a f_{\text{RPD}}(x) dx,$$

where $\bar{F}_{\text{RPD}}$ and $f_{\text{RPD}}$ denote respectively the survival function and the density function associated with the tail quantile function $U_{\text{RPD}}$ given in Definition 1. The change of variable $y \mapsto x = U_{\text{RPD}}(y)$ yields

$$ \begin{aligned} m_Y(a) &= (a+1)^2 \int_1^{+\infty} \log(U_{\text{RPD}}(y))y^{-a-2} dy, \\ &= (a+1)^2 \int_1^{+\infty} \log(U_{\text{RPD}}(1)y^{\gamma} \exp(\beta\gamma K_{\rho}(y))) y^{-a-2} dy, \\ &= (a+1)^2 \log(U_{\text{RPD}}(1)) \int_1^{+\infty} y^{-a-2} dy + (a+1)^2 \gamma \int_1^{+\infty} \log(y) y^{-a-2} dy \\ &\quad + (a+1)^2 \beta\gamma \int_1^{+\infty} K_{\rho}(y)y^{-a-2} dy. \end{aligned} $$

The first term vanishes given that $U_{\text{RPD}}(1) = 1$. Integrating by parts the second term concludes the proof. ■

### Proof of Proposition 4.

Suppose $U_{\text{RPD}}(1) = 1$, $a\gamma < 1$ and introduce $c = -a\beta\gamma/\rho$. Then, a direct calculation yields

$$ \begin{aligned} \text{CTM}_Y(a, \alpha) &= \frac{1}{\alpha} \int_0^\alpha U_{\text{RPD}}^a (1/v \mid \gamma, \rho, \beta) dv \\ &= \frac{1}{\alpha} \int_0^\alpha v^{-a\gamma} \exp(a\beta\gamma K_\rho(1/v)) dv \\ &= \frac{\exp(c)}{\alpha} \int_0^\alpha v^{-a\gamma} \exp(-cv^{-\rho}) dv \\ &= \exp(c) \alpha^{-a\gamma} \int_0^1 u^{-a\gamma} \exp(-c\alpha^{-\rho}u^{-\rho}) du, \end{aligned} $$

which is the desired result. Finally, note that, if $\beta > 0$, then $c > 0$ too and the CTM can be computed as

$$ \text{CTM}_Y(a, \alpha) = - \frac{\exp(c)}{\alpha \rho} c^{\frac{1-a\gamma}{\rho}} \Gamma_{\ell} \left( \frac{a\gamma - 1}{\rho}, c\alpha^{-\rho} \right), $$

where $\Gamma_{\ell}(\cdot, \cdot)$ is the incomplete lower gamma function. ■

### Proof of Lemma 1.

The result is a straightforward consequence of (23) in the proof of Proposition 2: Letting $\alpha = k/n$ and $y = d_n = k/(np_n)$ yields $U_{k/n}(d_n) = U(1/p_n)/U(n/k)$ and the result is proved. ■

### Proof of Theorem 5.

Let $d_n = k/(np_n)$ be the extrapolation factor. The following expansion holds:

$$ \begin{aligned} \log \hat{U}(1/p_n) - \log U(1/p_n) &= (\log \hat{U}(1/p_n) - \overline{\log U}(1/p_n)) + (\overline{\log U}(1/p_n) - \log U(1/p_n)) \\ &= A_{1,n} + A_{2,n} + A_{3,n} + B_n, \end{aligned} $$
with $A_{1,n} = \log X_{n-k,n} - \log U(n/k)$,

$$A_{2,n} = (\hat{\gamma} - \gamma) \log(d_n),$$

$$A_{3,n} = \hat{\gamma}\hat{\beta}_n K_{\hat{\rho}}(d_n) - \gamma\beta_n K_{\rho}(d_n),$$

$$B_n = \log U_{\text{RPD}}(d_n | \gamma, \rho, \beta(n/k)^\rho) - \log U_{k/n}(d_n).$$

Each term is considered separately. First, (de Haan and Ferreira, 2006, Theorem 2.4.1) yields

$$\frac{\sqrt{k_n}}{\log(d_n)} A_{1,n} = O_p(1/\log(d_n)), \qquad (24)$$

since $\sqrt{k}A(n/k) \to \lambda$ as $n \to \infty$. Second, letting $\xi_n := \sqrt{k}(\hat{\gamma} - \gamma)$, one has

$$\frac{\sqrt{k_n}}{\log(d_n)} A_{2,n} = \xi_n. \qquad (25)$$

Third, let us rewrite $A_{3,n}$ as

$$A_{3,n} = \beta_n \left( \hat{\gamma} \frac{\hat{\beta}_n}{\beta_n} K_{\hat{\rho}}(d_n) - \gamma K_{\rho}(d_n) \right) = \gamma \beta(n/k)^\rho (K_{\hat{\rho}}(d_n) O_P(1) - K_{\rho}(d_n)),$$

since $\hat{\gamma} \xrightarrow{\mathbb{P}} \gamma$ and $\hat{\beta}_n/\beta_n = O_P(1)$ by assumption. Remark that $K_\rho(d_n) \to -1/\rho$ as $n \to \infty$ while $|K_{\hat{\rho}}(d_n)| \le -2/\hat{\rho} \le -2/\rho^\dagger$ almost surely. As a consequence,

$$\frac{\sqrt{k_n}}{\log(d_n)} A_{3,n} = O_P(1/\log(d_n)), \qquad (26)$$

under the assumption $\sqrt{k}A(n/k) \to \lambda$ as $n \to \infty$. Finally, $B_n$ is a non-random term controlled with (de Haan and Ferreira, 2006, Eq. (3.2.7)): For all $\varepsilon > 0$, there exists $t_0$ such that for $t \ge t_0$ and $y \ge 1$,

$$\left| \frac{1}{A(t)} (\log U(ty) - \log U(t) - \gamma \log(y)) - K_{\rho}(y) \right| \le \varepsilon y^{\rho+\varepsilon}.$$

Taking account of (6) and considering $y = d_n \to \infty$ and $t = n/k \to \infty$, it follows

$$|\log U(1/p_n) - \log U(n/k) - \gamma \log(d_n) - \gamma \beta(n/k)^\rho K_\rho(d_n)| \le \varepsilon d_n^{\rho+\varepsilon} |A(n/k)|,$$

or equivalently, in view of Definition 1 and (23) in the proof of Proposition 2, $|B_n| \le \varepsilon d_n^{\rho+\varepsilon} |A(n/k)|$, so that

$$\frac{\sqrt{k_n}}{\log(d_n)} |B_n| \le \varepsilon |\lambda| \frac{d_n^{\rho+\varepsilon}}{\log(d_n)} (1 + o(1)), \qquad (27)$$

under the assumption $\sqrt{k}A(n/k) \to \lambda$ as $n \to \infty$. Combining (24)-(27), it follows that

$$\frac{\sqrt{k}}{\log(d_n)} \left( \log \hat{U}(1/p_n) - \log U(1/p_n) \right) = \xi_n + O_P \left( \frac{1}{\log(d_n)} \right) + O \left( \frac{d_n^{\rho+\varepsilon}}{\log(d_n)} \right),$$

Choosing $\varepsilon < -\rho$ concludes the proof. ■

**Proof of Lemma 2.**

Recall that

$$\text{CTM}_X(a, p_n) = \frac{1}{p_n} \int_0^{p_n} U^a(1/v)dv,$$
with $U(1/v) = U(n/k)U_{k/n}(k/(nv))$ from (23) in the proof of Proposition 2. Replacing, one obtains

$$ \mathrm{CTM}_X(a, p_n) = \frac{1}{p_n} U^a(n/k) \int_0^{p_n} U_{k/n}^a(k/(nv)) dv = d_n U^a(n/k) \int_0^{1/d_n} U_{k/n}^a(1/u) du, $$

with the change of variable $u = nv/k$. Remarking that

$$ d_n \int_0^{1/d_n} U_{k/n}^a(1/u) du = \mathrm{CTM}_{k/n}(a, 1/d_n) $$

concludes the proof.

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