Threshold Selection in Univariate Extreme Value Analysis

TL;DR

This paper introduces two new automated, parameter-free methods for threshold selection in univariate extreme value analysis, improving the estimation of rare event quantiles and adapting to changing data characteristics.

Contribution

The paper presents two novel threshold selection procedures that do not require manual tuning, enhancing the automation and robustness of extreme value analysis.

Findings

Both methods outperform existing procedures in simulations.

The methods provide reliable estimates across diverse distributions.

Application to financial data demonstrates practical utility.

Abstract

Threshold selection plays a key role for various aspects of statistical inference of rare events. Most classical approaches tackling this problem for heavy-tailed distributions crucially depend on tuning parameters or critical values to be chosen by the practitioner. To simplify the use of automated, data-driven threshold selection methods, we introduce two new procedures not requiring the manual choice of any parameters. The first method measures the deviation of the log-spacings from the exponential distribution and achieves good performance in simulations for estimating high quantiles. The second approach smoothly estimates the asymptotic mean square error of the Hill estimator and performs consistently well over a wide range of distributions. The methods are compared to existing procedures in an extensive simulation study and applied to a dataset of financial losses, where the underlying extreme value index is assumed to vary over time. This application strongly emphasizes the importance of solid automated threshold selection.

Tables

Table 1. Table 1: The averages of adaptive γ 𝛾 \gamma estimates and their root mean square error (RMSE) in brackets are presented for thresholds k ^ ^ 𝑘 \hat{k} that are selected using SAMSEE or SAMSEE ρ ^ subscript SAMSEE ^ 𝜌 \mathrm{SAMSEE}_{\hat{\rho}} with the true ρ 𝜌 \rho , ρ = − 1 𝜌 1 \rho=-1 or ρ ^ = ρ ^ ( 1 ) ^ 𝜌 superscript ^ 𝜌 1 \hat{\rho}=\hat{\rho}^{(1)} .

			$𝔼\left[{\hat{\gamma }}_{\hat{k}}\right]$ (RMSE)
	$\gamma$	$\rho$	true $\rho$	$\rho =-1$	$\hat{\rho }$
Student-t(6)	0.17	-1/3	0.21 (0.09)	0.26 (0.12)	0.28 (0.14)
Fréchet(2)	0.50	-1	0.51 (0.07)	0.51 (0.07)	0.51 (0.08)
Cauchy	1.00	-2	1.01 (0.13)	0.97 (0.17)	0.99 (0.16)
Burr(2,1)	2.00	-1	2.05 (0.34)	2.05 (0.34)	2.03 (0.40)

Table 2. Table 2: Efficiency values EFF γ subscript EFF 𝛾 \mathrm{EFF}_{\gamma} based on 2000 samples if n = 500 𝑛 500 n=500 and on 500 samples if n = 5000 𝑛 5000 n=5000 . Lowest (best) efficiency values are highlighted in blue.

$n=500$	SAM	GH	DK	GO	DB	B	RT	sIHS
Student-t(6)	1.07	1.68	1.18	1.38	1.06	1.04	1.04	1.14
Fréchet(2)	1.13	1.15	1.08	1.12	1.60	1.49	2.00	1.41
Cauchy	1.37	1.19	1.32	1.16	2.14	1.85	2.11	1.47
Loggamma	0.98	1.06	1.27	1.11	1.12	1.04	1.32	0.78
Burr(2,1)	1.11	1.22	1.47	1.13	1.68	1.42	1.82	1.14
negBias	1.06	1.13	1.56	1.13	1.07	1.22	1.89	2.27

Table 3. Table 3: Efficiency values EFF q subscript EFF 𝑞 \mathrm{EFF}_{q} for p = 0.001 𝑝 0.001 p=0.001 based on 2000 samples if n = 500 𝑛 500 n=500 and on 500 samples if n = 5000 𝑛 5000 n=5000 . Lowest (best) efficiency values are highlighted in blue.

$n=500$	SAM	GH	DK	GO	DB	B	RT	sIHS
Student-t(6)	1.09	2.30	1.23	1.60	1.01	1.04	1.16	1.10
Fréchet(2)	0.96	1.07	1.01	1.06	1.07	1.16	1.42	0.83
Cauchy	0.89	1.04	1.03	0.95	0.86	1.40	1.59	0.65
Loggamma	0.84	0.95	2.10	1.02	0.88	1.06	1.55	0.50
Burr(2,1)	0.79	2.15	8.60	0.98	0.71	1.43	3.19	0.41
negBias	1.66	1.37	2.32	2.13	0.80	1.98	3.75	8.05

Equations

$$\underset{n\to \infty }{\lim }{F}^n(a_{n}x+b_n)=G_{\gamma }(x):=\exp \left(-x^{-1/\gamma }\right).$$

$$\underset{t\to \infty }{\lim }\frac{1-F(tx)}{1-F(t)}=x^{-1/\gamma },$$

$$\displaystyle\frac{X_{1}}{t}\,\Big{|}\,X_{1}>t\ \overset{\mathcal{D}}{\longrightarrow}\$$

$$\displaystyle\log\left(\frac{X_{1}}{t}\right)\Big{|}\,X_{1}>t\ \overset{\mathcal{D}}{\longrightarrow}\$$

$${\hat{\gamma }}_k:=\frac{1}{k}\sum _{i=1}^k\log \left(\frac{X_{(n-i+1,n)}}{X_{(n-k,n)}}\right)=\frac{1}{k}\sum _{i=1}^k{Y}_{(i,k)}.$$

$$\mathrm{ISE}(k):=\int {\left(h_{\gamma }(x)-h_{{\hat{\gamma }}_k}(x)\right)}^2\mathrm{d}x=\frac{1}{2\gamma }-\frac{2}{\gamma +{\hat{\gamma }}_k}+\frac{1}{2{\hat{\gamma }}_k}.$$

$$\mathrm{MISE}(k)-\frac{1}{2\gamma }:={\mathbb{E}}_{H_0}[\mathrm{ISE}(k)]-\frac{1}{2\gamma }=-\frac{1}{\gamma }C(k)+\frac{k}{2(k-1)\gamma },$$

$${\mathbb{E}}_{H_0}\left[\frac{2}{\gamma +{\hat{\gamma }}_k}\right]\approx \frac{1}{\gamma }={\mathbb{E}}_{H_0}\left[\frac{k-1}{k{\hat{\gamma }}_k}\right].$$

$${\mathbb{E}}_{H_0}\left[{\left(\frac{w}{{\hat{\gamma }}_k}-\frac{2}{{\hat{\gamma }}_k+\gamma }\right)}^2\right]\approx \frac{w^2{k}^2}{{\gamma }^2(k-1)(k-2)}-\frac{2wk}{{\gamma }^2(k-1)}+\frac{1}{{\gamma }^2}.$$

$$\mathrm{IHS}(k):=\frac{1}{2{\hat{\gamma }}_k}-\frac{k-2}{{\hat{\gamma }}_{k}k}=\frac{4-k}{2{\hat{\gamma }}_{k}k}$$

$${\hat{k}}_{\mathrm{IHS}}:={\mathrm{argmin}}_{1<k<n}\mathrm{IHS}(k).$$

$${\mathrm{IHS}}^{-}(k):=\frac{4+k}{2{\hat{\gamma }}_{k}k}\text{and}{\hat{k}}_{{\mathrm{IHS}}^{-}}:={\mathrm{argmin}}_{1<k<n}{\mathrm{IHS}}^{-}(k).$$

$$\mathrm{IHS}(k)=\mathbb{E}[\mathrm{IHS}](k)+\sigma {ϵ}_k,k=1,\dots ,n,$$

$$\mathrm{IHS}+\frac{1}{2\gamma }=\mathrm{ISE}+\frac{2}{k{\hat{\gamma }}_k}+\frac{{\hat{\gamma }}_k-\gamma }{{\hat{\gamma }}_k({\hat{\gamma }}_k+\gamma )}.$$

$$\underset{t\to \infty }{\lim }\frac{\frac{U(tx)}{U(t)}-x^{\gamma }}{A(t)}=x^{\gamma }\frac{x^{\rho }-1}{\rho },$$

$$\sqrt{k}({\hat{\gamma }}_k-\gamma )\stackrel{\mathcal{D}}{⟶}\mathcal{N}\left(\frac{\lambda }{(1-\rho )},{\gamma }^2\right),$$

$$\sqrt{k}\left(\frac{1}{{\hat{\gamma }}_k}-\frac{1}{\gamma }\right)\stackrel{\mathcal{D}}{⟶}\mathcal{N}\left(\frac{-\lambda }{(1-\rho ){\gamma }^2},\frac{1}{{\gamma }^2}\right).$$

$$k_{\mathrm{opt}}=\left[{\left(\frac{{\gamma }^2(1-\rho {)}^2}{-2\rho c^2}\right)}^{1/(1-2\rho )}{n}^{-2\rho /(1-2\rho )}\right].$$

$$k_{\mathrm{IHS}}$$

$$\approx {\mathrm{argmin}}_k\left\{\frac{2}{\gamma k}+\frac{A(n/k)}{2{\gamma }^2(1-\rho )}\right\}=\left[{\left(\frac{4\gamma (1-\rho )}{-\rho c}\right)}^{1/(1-\rho )}{n}^{-\rho /(1-\rho )}\right].$$

$$\frac{k_{\mathrm{IHS}}}{k_{\mathrm{opt}}}\approx {\left(\frac{-\rho }{32}\cdot k_{\mathrm{IHS}}\right)}^{-1/(1-2\rho )}\approx d\cdot n^{\frac{\rho }{(1-2\rho )(1-\rho )}}⟶0,$$

$${\hat{q}}_k(p)=X_{(n-k,n)}{\left(\frac{k}{np}\right)}^{{\hat{\gamma }}_k}.$$

$$\mathrm{MSEQ}:=\mathbb{E}\left[{\left(\frac{{\hat{q}}_k(p)}{q(p)}-1\right)}^2\right]/\log \left(\frac{k}{np}\right),$$

$$\mathbb{AE}[({\hat{\gamma }}_k-\gamma {)}^2]=\frac{{\gamma }^2}{k}+\frac{A(n/k{)}^2}{(1-\rho {)}^2}.$$

$$M_{n,k}:=\frac{1}{k}\sum _{i=1}^k{Y}_{(i,k)}^2,{\hat{\gamma }}_{\mathrm{V},k}:=\frac{M_{n,k}}{2{\hat{\gamma }}_k},\text{ and}{\hat{\gamma }}_k^{\mathrm{GJ}}:=2{\hat{\gamma }}_{\mathrm{V},k}-{\hat{\gamma }}_k,$$

$${\bar{\gamma }}_k:=\frac{1}{k}\sum _{i=1}^k{\hat{\gamma }}_i\text{and}{\bar{\gamma }}_{\mathrm{up},K,k}:=\frac{1}{K-k+1}\sum _{i=k}^K{\hat{\gamma }}_i,$$

$${\bar{b}}_{\mathrm{up},K,k}:={\bar{\gamma }}_{\mathrm{up},K,k}-{\bar{\gamma }}_K.$$

$$\mathbb{AE}[{\bar{b}}_{\mathrm{up},K,k}]=\frac{-\rho A(n/k)}{(1-\rho {)}^2}=\frac{1}{2}\frac{A(n/k)}{(1-\rho )}$$

$$\mathbb{AE}[{\hat{\gamma }}_k]=\mathbb{AE}[{\hat{\gamma }}_{\mathrm{V},k}+{\bar{b}}_{\mathrm{up},K,k}].$$

$$E^2(K):=\frac{1}{K}\sum _{k=1}^K{\left({\hat{\gamma }}_{\mathrm{V},k}+{\bar{b}}_{\mathrm{up},K,k}-{\hat{\gamma }}_k\right)}^2$$

$$\displaystyle K^{*}:=\underset{K}{\mathop{\rm arg\min}}\left\{\sum_{L=K-2}^{K+2}\Big{(}E^{2}(K)-E^{2}(L)\Big{)}^{2}\right\}.$$

Full text

The methods are compared to existing procedures in an extensive simulation study and applied to a dataset of financial losses, where the underlying extreme value index is assumed to vary over time. This application strongly emphasizes the importance of solid automated threshold selection.

AMS 2010 Subject Classification:

Primary 62G32; secondary 62G05, 62F12, 97M30

Keywords:

Extreme value statistics, peak-over-threshold approach, power laws, Hill estimator, tuning parameter selection, bias estimation

1 Introduction

Extreme value analysis of heavy-tailed distributions is an important model in various applications. In seismology and climatology, for example, statistics of extremes is used to study earthquakes (Beirlant et al., 2018) or heavy precipitation (Carreau et al., 2017). Another important field of research is analysing high financial losses, which becomes particularly interesting if the losses depend on covariates (Chavez-Demoulin et al., 2016; Hambuckers et al., 2018). In this situation an automated threshold selection procedure could bring additional benefits by enabling the selection of the threshold depending on a covariate. We will discuss this possibility in more detail in Section 5.

To mathematically investigate the behaviour of heavy tails, we consider random variables from the domain of attraction (DoA) of a Fréchet distribution. Let $X_{1},\dots,X_{n}$ be independent identically distributed (i.i.d.) random variables with distribution function $F$, where $F$ is in the DoA of an extreme value distribution (evd) $G_{\gamma}$ with extreme value index $\gamma>0$. This means there exist sequences $a_{n}>0$ and $b_{n}$ real, s.t.

[TABLE]

In this situation the following first order condition holds,

[TABLE]

i.e. the survival function $1-F$ is regularly varying with index $-1/\gamma$. Distributions fulfilling this condition are called Pareto-type distributions, because they only differ from the Pareto distribution by a slowly varying function $\ell_{F}(x)$, i.e. $1-F(x)=x^{-1/\gamma}\ell_{F}(x)$.

We can interpret the quotient in (1) as a conditional probability, and it follows directly that

[TABLE]

Thus, for a sufficiently large threshold $t$ the data above this threshold can be modelled by a Pareto or an exponential distribution. In this article we concentrate on the exponential approximation and utilize it for inference on the extreme value index. It is common to consider the threshold $t=X_{(n-k,n)}$ and choose the sample fraction $k$ instead of $t$, where $X_{(1,n)}\leq\cdots\leq X_{(n,n)}$ denote the order statistics of a sample of size $n$. In this case, a natural estimator for $\gamma$ under the exponential approximation of the log-spacings $Y_{(i,k)}:=\log(X_{(n-i+1,n)})-\log(X_{(n-k,n)})$ is their mean, the Hill estimator (Hill, 1975),

[TABLE]

The Hill estimator is still among the most popular and well-known estimators for the extreme value index, although its sample path as a function in $k$ can be highly unstable and estimation therefore crucially depends on the choice of the sample fraction $k$. This dependence highlights the difficulties in estimating $\gamma$: even from univariate i.i.d. observations from $F\in\mathrm{DoA}(G_{\gamma})$, estimation is hard, since only few observations contain information about the extreme value distribution $G_{\gamma}$. To select a threshold above which the data can be used for statistical inference about the tail is one of the most fundamental problems in the field of extreme value analysis.

Due to the importance of this task, the appropriate choice of the threshold has been discussed extensively in extreme value research over the last decades, and suggested solutions cover a variety of methodologies. We give a short summary on different types of approaches and stress the specific difficulties that arise. We mainly concentrate on methods we compare in our simulation study in Section 4. More comprehensive reviews about threshold selection can be found in Scarrott and MacDonald (2012) and Dey and Yan (2016).

One basic concept in threshold selection is data visualisation, which is also discussed more deeply in Kratz and Resnick (1996) and Drees et al. (2000). Popular graphical diagnostics used in this context are the Zipf plot, Hill plot, QQ-plot or the mean-excess plot to name a few. A major drawback of these methods is their subjectivity due to the necessarily personal interpretation of the plot. Further, it is a burden to choose each threshold manually, especially in high dimensional settings or when analysing many samples. Easier ways to select the sample fraction are rules-of-thumb such as using the upper 10% of the data (DuMouchel, 1983) or $k=\sqrt{n}$ (Ferreira et al., 2003). However, these suggestions are neither theoretically justified nor data driven. Reiss and Thomas (2007) present a procedure that tries to find a region of stability among the estimates of the extreme value index. Their method depends on a tuning parameter, whose choice is further analysed in Neves and Fraga Alves (2004). To our knowledge no theoretical analysis exists for this approach.

Besides these and similar heuristic approaches, there is a class of theoretically motivated procedures that target the optimal sample fraction for specific estimation tasks, such as quantile estimates (Ferreira et al., 2003), estimation of high probabilities (Hall and Weissman, 1997) or the Hill estimator, see below. We also mention two other methodologies. First, there are suggestions that utilize comparing the empirical distribution to the fitted generalized Pareto distribution (GPD) via goodness-of-fit tests (Bader et al., 2018) or by minimizing the distance between them (Pickands, 1975; Gonzalo and Olmo, 2004; Clauset et al., 2009), where the latter approach is theoretically analysed by Drees et al. (2018). Further, Goegebeur et al. (2008) propose a family of kernel statistics to test for exponentiality in order to select a threshold.

Of particular interest to us are methods that aim to estimate the sample fraction $k_{\mathrm{opt}}$ which minimizes the asymptotic mean square error (AMSE) of the Hill estimator. To construct an estimator for $k_{\mathrm{opt}}$, Drees and Kaufmann (1998) utilize the Lepskii method and an upper bound on the maximum random fluctuation of $\hat{\gamma}_{k}$ around $\gamma$. To apply their approach it is necessary to choose several tuning parameters and to obtain consistent initial estimates for $\gamma$ and a second order parameter $\rho$. They recommend specific choices of the parameters based on a numerical study and we employ their proposals in our simulations. However, the choice of these parameters is not data-driven. In Guillou and Hall (2001), a test statistic $Q_{k}$ is constructed based on an accumulation of log-spacings, which takes values around 1 as long as the bias of the Hill estimator is not significantly large. Their statistic depends on a tuning parameter as well, and a critical value to test $Q_{k}$ against has to be chosen. Again we adopt the parameter choice suggested in their simulation study. Danielsson et al. (2001) introduce a double bootstrap approach to estimate the optimal sample fraction. They need to choose the number of bootstrap samples and a parameter $n_{1}$. For $n_{1}$, a data-driven but computationally expensive selection method is provided, where the whole bootstrap procedure is repeated for various possible values of $n_{1}$. Another estimator for $k_{\mathrm{opt}}$ is given by Beirlant et al. (2002), which employs least squares estimates from an exponential regression approach. The method depends on an estimate for $\rho$ and a sample fraction $k_{0}$. To avoid the choice of $k_{0}$ they suggest taking the median of the estimates over a range of values, e.g. $k_{0}\in\{3,\dots,n/2\}$. A different approach is taken by Goegebeur et al. (2008), who use the properties of a test statistic regarding bias estimation to construct an estimator for the AMSE$/\gamma$ and minimize it with respect to $k$. If one fixes $\rho=-1$, as they suggest in their simulations chapter, there is no further tuning parameter to be chosen. However, no result about consistency of $\hat{k}$ in the sense of $\hat{k}/k_{\mathrm{opt}}\overset{\mathbb{P}}{\rightarrow}1$ is known in contrast to the approaches in Drees and Kaufmann (1998), Guillou and Hall (2001), Danielsson et al. (2001) and Beirlant et al. (2002).

In this paper we contribute to the problem of threshold selection by introducing two new methods. The first one presented in Section 2 is inspired by the idea of testing the exponential approximation. We estimate the integrated square error (ISE) of the exponential density under the assumption that the log-spacings are indeed exponentially distributed. The error functional we obtain, denoted as inverse Hill statistic (IHS), is very easy to compute and does not depend on any tuning parameters. Since this criterion is variable for small $k$, it can be additionally smoothed to improve the performance. The minimizing sample fraction of IHS is asymptotically smaller than $k_{\mathrm{opt}}$, as it is stricter against deviation from the exponential approximation. This estimator performs remarkably well for adaptive quantile estimation on finite samples, as illustrated in our simulation study.

In our second approach we suggest a smooth estimator for the AMSE of the Hill estimator, called SAMSEE (smooth AMSE estimator). This estimator is constructed by a preliminary estimate of $\gamma$ using the generalized Jackknife approach in Gomes et al. (2000) and a bias estimator for the Hill estimator introduced in Section 3. By minimizing SAMSEE we estimate the optimal sample fraction $k_{\mathrm{opt}}$. For estimation, the choice of a large sample fraction $K$ is necessary, for which we present a data-driven selection procedure in Section 3. SAMSEE utilizes the idea of fixing $\rho=-1$, which is justified by good performance in simulations and leads to a simpler and more robust estimator. However, the estimator can also be adjusted to any $\rho$ by including a consistent estimator $\hat{\rho}$, as described in Section 3.1.

After introducing our two novel threshold selection methods in Sections 2 and 3 we compare these methods to various other approaches in an numerical analysis in Section 4. In Section 5 the importance of automated threshold selection procedures is illustrated in an application, where we non-parametrically estimate an extreme value index that varies over time. The proof of Theorem 3, which describes the asymptotic behaviour of our bias estimator, and auxiliary theoretical results can be found in Appendix A.

2 IHS – The inverse Hill statistic

In this section we introduce the first threshold selection procedure by analysing the integrated square error (ISE) between the exponential density $h_{\gamma}$ and its parametric estimator $h_{\hat{\gamma}_{k}}$ employing the Hill estimator,

[TABLE]

The first term of ISE is constant and thus plays no role for selecting $k$. The last term of ISE is known, but the second term is not. Therefore, we cannot minimize ISE directly. Instead, we want to estimate and minimize its expectation under the exponential approximation. This is based on the idea of considering the hypothesis $H_{0}$ that the log-spacings $Y_{(i,k)}$ are indeed exponentially distributed. Under $H_{0}$ the Hill estimator is gamma distributed, see Lemma 1, and the mean of ISE (MISE) can be calculated explicitly. We observe that MISE is a decreasing function in $k$ under the exponential approximation,

[TABLE]

where $C(k):=2\exp(k)k^{k}\Gamma(1-k,k)$ and $\Gamma(a,b)$ denotes the upper incomplete gamma function. The function $C(k)$ converges to 1 very fast, s.t. we obtain

[TABLE]

This provides us with an unbiased estimator for the first term in (4) under $H_{0}$. However, due to the high variability for small $k$, we instead want to find an estimator of the form $w/\hat{\gamma}_{k}$ for some $w$ depending on $k$ that minimizes the MSE under the exponential approximation. To do so, we approximate its MSE in the following way,

[TABLE]

The approximation depends on similar functions as $C(k)$, which quickly become constant. The MSE in (5) is minimized for $w=(k-2)/k$. Thus, we suggest the inverse Hill statistic

[TABLE]

to estimate $\mathrm{MISE}(k)-(2\gamma)^{-1}$ and the threshold selected via minimizing IHS,

[TABLE]

By minimizing IHS we select a sample fraction where IHS starts increasing and contradicts $H_{0}$ by behaving contrarily to MISE under the exponential approximation. This criterion can be compared to hypothesis testing with a large significance level $\alpha$, which implies seeking high confidence when deciding to not reject $H_{0}$. Further properties of $\hat{k}_{\mathrm{IHS}}$ are analysed theoretically in Section 2.1 and for finite samples in a numerical study in Section 4.

Note that the performance of IHS depends on the bias of the Hill estimator being positive and increasing, see Section 2.1. However, the bias can be negative for some non-standard distributions. In case of a negative bias, we instead suggest to use,

[TABLE]

The two cases can easily be distinguished by analysis of the Hill estimator for large $k$. Both $\mathrm{IHS}$ and $\mathrm{IHS}^{-}$ are justified by asymptotic results in Section 2.1.

Figure 1 illustrates that IHS is highly varying for small $k$, which makes automatic threshold choices more variable. To control this problematic behaviour we smooth the IHS. More specifically, we want to estimate $\mathbb{E}[\mathrm{IHS}]$ by considering the regression problem

[TABLE]

where $\sigma>0$ and $\mathbb{E}[\epsilon_{k}]=0$. Due to the structure of the Hill estimator, the random variables $\epsilon_{k}$ are highly dependent, which needs to be taken into account in estimation. In our simulations, we apply a Bayesian non-parametric procedure introduced by Serra et al. (2018) which simultaneously estimates mean and covariance and is available in the R-package eBsc. The approach provides a smooth estimator for the expectation of IHS – denoted as sIHS – comprising less variation for small $k$. This way we can improve the performance by selecting a more suitable threshold, as illustrated in Figure 1. Of course, one can also use other smoothing procedures suitable for dependent data (Opsomer et al., 2001; Krivobokova and Kauermann, 2007; Lee et al., 2010).

We finally want to remark on the relation between IHS and ISE, which is given by

[TABLE]

This equation points out that minimizing IHS does not minimize ISE, as IHS takes an additional bias term into account. If the bias of the Hill estimator is positive, IHS selects smaller $k$ (larger thresholds) than ISE. This is not surprising, because we estimate the expectation of the ISE under the hypothesis that the exponential approximation holds. This is a much more conservative error functional, meaning it is more strict against deviation from the exponential distribution.

In conclusion, with IHS we do not aim to estimate $k_{\mathrm{opt}}$ but to find a sample fraction where we can be very certain that the exponential approximation still holds. The impact of this consideration is illustrated in simulations and an application in Sections 4 and 5.

2.1 Theorectical analysis of IHS

In order to understand the IHS asymptotically we consider the second order condition,

[TABLE]

for $x>0$ and with second order parameter $\rho<0$. Here, $A(t)$ denotes a function converging to zero as $t$ goes to infinity and $|A|$ is regularly varying with index $\rho$. Further, $U$ is defined by $U(x):=F^{\leftharpoonup}\!\left(1-\frac{1}{x}\right)$, where $F^{\leftharpoonup}$ denotes the left inverse of the distribution function $F$. In this setting the following asymptotic normality statements for the Hill estimator $\hat{\gamma}_{k}$ hold.

Theorem 1** **(Theorem 3.2.5 in de Haan and

Ferreira (2006)).

Let $X_{1},\dots,X_{n}$ be i.i.d. random variables with distribution function $F\in\mathrm{DoA}(G_{\gamma})$ for $\gamma>0$. If (7) holds and $k$ is an intermediate sequence, i.e. $k\rightarrow\infty$ and $k/n\rightarrow 0$ as $n\rightarrow\infty$, then

[TABLE]

with $\lambda:=\underset{k\rightarrow\infty}{\lim}\sqrt{k}A(n/k)$.

Theorem 2**.**

Under the conditions of Theorem 1, it holds that

[TABLE]

Proof.

Applying the delta method to Thm. 1. ∎

Following the reasoning in de Haan and Ferreira (2006), page 78, the minimizing point of the AMSE can be found explicitly if considering $A(t)=ct^{\rho}$ with $c\neq 0$. In this special case the minimizing sample fraction can be expressed as

[TABLE]

Under the same assumption we can calculate the minimizing point $k_{\mathrm{IHS}}$ of the asymptotic expectations of $\mathrm{IHS}$ and $\mathrm{IHS}^{-}$. Let $\mathbb{AE}$ denote the asymptotic expectation referring to the expectation of the limiting distribution in Thm. 2. Then

[TABLE]

It is easy to check that the same formula holds for $\mathrm{IHS}^{-}$ if $c$ is replaced by its absolute value. Further note that by Lemma 2 it is sufficient to consider intermediate sequences when determining the minimizing sequence. Comparing $k_{\mathrm{opt}}$ and $k_{\mathrm{IHS}}$ for a fixed $\rho>-\infty$ we obtain that

[TABLE]

as $n\rightarrow\infty$ and for a constant $d$ depending on $\rho$, $\gamma$ and $c$. This supports what equation (6) already suggested: minimizing $\mathrm{IHS}$ gives asymptotically a smaller $k$ than $k_{\mathrm{opt}}$. Thus, $k_{\mathrm{IHS}}$ asymptotically performs suboptimally for the Hill estimator but still leads to a consistent sequence of estimates. For finite samples the ratio crucially depends on $\rho$, and $k_{\mathrm{IHS}}$ can be even larger than $k_{\mathrm{opt}}$, as illustrated in Figure 2. The graphic presents the quotient of the two sample fractions as a function in the second order parameter $\rho$ for different samples sizes. The parameters $c$ and $\gamma$ are fixed to 1, as they have a weaker impact on the proportion. It also holds that $k_{\mathrm{IHS}}/k_{\mathrm{opt}}\rightarrow 1$, as $\rho\rightarrow-\infty$, since both sample fractions converge to $n$ in this case.

Although $k_{\mathrm{IHS}}$ is of smaller order than $k_{\mathrm{opt}}$ asymptotically, the simulation study in Section 4 shows that $\hat{k}_{\mathrm{IHS}}$ works remarkably well when used for quantile estimation. We consider the following quantile estimator for the $(1-p)$-quantile,

[TABLE]

The sample fraction $k_{\mathrm{opt}}$ also minimizes the asymptotic relative MSE of $\hat{q}_{k}(p)$, see e.g. Theorem 4.3.8 in de Haan and Ferreira (2006). For finite samples however, the quantile estimator seems to benefit from $k_{\mathrm{IHS}}$. This has different reasons, two of which are illustrated by Figure 3. On the left we see the empirical expectation of IHS, the empirical versions of the MSE of $\hat{\gamma}_{k}$ and the relative MSE of the quantile estimator,

[TABLE]

as used in Theorem 4.3.8 in de Haan and Ferreira (2006). We observe that $k_{\mathrm{IHS}}$ (blue dot) is indeed smaller than $k_{\mathrm{opt}}$ (black) but so is the minimizer of MSEQ (pink) as well.

On the right we see a plot of the empirical $\mathbb{E}[\mathrm{IHS}]$ and MSE of $\hat{\gamma}_{k}$ for Loggamma distributed samples of size 5000. This graphic highlights the similarities between MSE and IHS for the boundary case $\rho=0$.

These observations indicate why $\hat{k}_{\mathrm{IHS}}$ outperforms other methods that try to minimize the MSE of the Hill estimator when adaptively estimating $q(p)$ by (10) on most of our exemplary distributions and sample sizes $n=500$ and $n=5000$, see Section 4.

3 SAMSEE - The smooth AMSE estimator

In this section we illustrate a way to smoothly estimate the AMSE of the Hill estimator. Via minimizing this AMSE estimator, called SAMSEE, we obtain an estimator for $k_{\mathrm{opt}}$. By this means, we extend previous methods which also estimate $k_{\mathrm{opt}}$ by estimating the AMSE itself. From Thm. 1 it is easy to see that the AMSE, which is the asymptotic variance plus the asymptotic squared bias, equals

[TABLE]

Thus, to estimate the AMSE as a function in $k$ we employ two estimators, one for $\gamma$ and one for the bias term as a combination of $\rho$ and $A$. First we explain how we estimate $\gamma$ and then we define the bias estimator. This bias estimator has a quite smooth sample path in $k$, and it depends on the choice of a large sample fraction $K$, for which we afterwards provide a data-driven selection procedure.

Note that, for the moment, we assume that the second order parameter $\rho$ is equal to $-1$ to motivate the construction of the AMSE estimator. The idea of misspecifying $\rho$ to simplify estimation – via avoiding the additional uncertainty through estimating $\rho$ or selecting an influential tuning parameter – was already used, for example, by Gomes et al. (2000), Drees and Kaufmann (1998) and Goegebeur et al. (2008). It is also motivated by the simulations in Section 3.1.

For $\gamma$ we consider the generalized Jackknife estimator $\hat{\gamma}_{k}^{\mathrm{GJ}}$ introduced by Gomes et al. (2000) as $\gamma_{n}^{\mathrm{G}_{1}}$. This estimator is defined by

[TABLE]

where $Y_{(i,k)}$ denotes the log-spacings as in equation (3). Note, that $\hat{\gamma}_{V,k}$ is the de Vries estimator introduced under this name in de Haan and Peng (1998) and $\hat{\gamma}_{k}$ is the Hill estimator as above. The generalized Jackknife estimator has a reduced bias compared to the Hill estimator and is even asymptotically unbiased if $\rho=-1$, see (2.11) in Gomes et al. (2000). This property is useful here, since the bias estimator $\bar{b}_{\mathrm{up},K,k}$ defined in the following performs optimally for $\rho=-1$ as well. Furthermore, the same large sample fraction $K$ can be used for $\hat{\gamma}_{K}^{\mathrm{GJ}}$ and $\bar{b}_{\mathrm{up},K,k}$.

To construct this bias estimator, we study the following averages of Hill estimators,

[TABLE]

where $k<K$. Plotting these averages illustrates how they smoothly frame the sample path of the Hill estimator. Especially the upper mean $\bar{\gamma}_{\mathrm{up},K,k}$ seems to contain a lot of structural information about the underlying asymptotic bias of the Hill estimator when choosing the upper bound $K$ appropriately, see Figure 4.

This similarity between the upper mean and the bias of the Hill estimator inspires the definition

[TABLE]

The estimator $\bar{b}_{\mathrm{up},K,k}$ is indeed a sensible estimator for a bias function, since

[TABLE]

follows for $\rho=-1$ from Theorem 3.

Danielsson et al. (2001) use $(\hat{\gamma}_{\mathrm{V},k}-\hat{\gamma}_{k})$ to access the bias of $\hat{\gamma}_{k}$ and apply a double bootstrap procedure to stabilize this highly varying estimate. We use the difference of two estimators for $\gamma$ as well, but now consider averaging to smooth the bias estimate. The idea to average the Hill estimator in order to smooth the Hill plot and decrease the variance is also studied in Resnick and Stǎricǎ (1997).

It remains to choose an appropriate $K$ in order to complete SAMSEE and to estimate the optimal sample fraction $k_{\mathrm{opt}}$. We need $K$ to be large enough to allow for minimization over all relevant $k$ and small enough to be an intermediate sequence itself (see Theorem 3 for this condition). To find such a $K$ we use the following relation between the estimators,

[TABLE]

This provides us with a relatively stable function in $k$, $\hat{\gamma}_{\mathrm{V},k}+\bar{b}_{up,K,k}$, that has the same asymptotic expectation as the highly non-smooth Hill estimator. We want to find an intermediate sequence $K$ for which (16) holds and thus define

[TABLE]

to measure the deviation from approximation (16) uniformly over all $k\leq K$. Based on this, we suggest to choose

[TABLE]

In this way we select a $K^{*}$ where the asymptotic approximation (16) is most stable, since we minimize the local variation of $E^{2}(K)$. Simulations suggest that this criterion is not sensitive to slightly increasing the region of stability $\{K-h,\dots,K+h\}$ from $h=2$ to $h=5$ or $10$ depending on the sample size.

Now we finally combine the previously described estimators to approach the AMSE in (12) under the assumption that $\rho=-1$. With $K^{*}$ in (18) and the property of $\bar{b}_{\mathrm{up},K,k}$ in (15), we obtain an estimator for the AMSE of the Hill estimator and for $k_{\mathrm{opt}}$ by

[TABLE]

Figure 5 illustrates how such a smooth estimate of the AMSE can look like. On the left, SAMSEE is displayed for a Fréchet sample with parameters $\gamma=1/2$ and $\rho=-1$. On the right, the Hill plot of the same sample is presented for all $k\leq K^{*}=388$.

This smooth estimate of the AMSE can be useful beyond the context of threshold selection. For extreme value mixture models or Bayesian threshold selection approaches, SAMSEE could be used to construct a transition function between bulk and tail distribution or an empirical prior for the threshold, respectively, see Scarrott and MacDonald (2012) for a review on mixture models.

3.1 SAMSEE if $\rho\neq-1$

We next want to analyse SAMSEE in the broader context of an unknown second order parameter $\rho$. The first thing to note is that the generalized Jackknife estimator is no longer unbiased in this situation. Secondly, the behaviour of our bias estimator changes, as it is described in the following Theorem.

Theorem 3**.**

Under the conditions of Theorem 1 and for $k/K\rightarrow c$ with $0<c<1$ as $n\rightarrow\infty$, it holds for $\bar{b}_{\mathrm{up},K,k}$ in (14) that

[TABLE]

where $\delta_{\rho}(c)=(c^{\rho}-1)/(-\rho(c^{-1}-1))$ and $\nu(c)=2c^{2}/(1-c)^{2}\cdot(1-c+c\log(c))$ with $0\leq\nu(c)\leq 1$.

Proof.

The proof can be found at the end of Appendix A. ∎

From Theorem 3 follows that

[TABLE]

For $\rho=-1$ the function $\delta_{-1}(c)$ is equal to 1. If $\rho\neq-1$, we can observe that $\delta$ bends our bias estimator and it will therefore increase slightly too fast or too slow. We can still apply SAMSEE in this situation and select $K^{*}$ from (18). However, approximation (16) does not hold anymore and instead the following holds,

[TABLE]

The absolute value of the error described by (20) is high if $\delta_{\rho}$ strongly differs from $1$ and the bias term $A$ is large. If $\rho\neq-1$, $\delta_{\rho}$ indeed deviates from $1$ and we minimize the error by minimizing the bias. This is why applying (18) in this case leads to a small $K^{*}$. On the other hand, if $\delta_{\rho}=1$, the approximation stays valid for an increasing bias and $K^{*}$ will typically be larger.

An alternative to fixing $\rho=-1$ is to incorporate a consistent estimator $\hat{\rho}$ of the second order parameter. This can be done via

[TABLE]

and

[TABLE]

In this way we can construct an estimator for $k_{\mathrm{opt}}$ in the general setting of Pareto-type distributions.

In Table 1, we present the results of a simulation study indicating for which distributions it is beneficial to use $\hat{\rho}$ instead of $\rho=-1$. We estimate $\rho$ using the estimator $\hat{\rho}^{(1)}$ suggested in Theorem 1 in Drees and Kaufmann (1998). The results indicate that, in general, it is sensible to fix $\rho=-1$ in SAMSEE, since only for the Cauchy distribution using $\hat{\rho}$ performs slightly better regarding bias and RMSE. This confirms the observations already made by others (Gomes et al., 2000; Drees and Kaufmann, 1998; Goegebeur et al., 2008), that it is often recommendable to select $\rho=-1$ instead of allowing for further variability by including an additional estimator.

4 Simulation study

In the following we numerically analyse the performance of eight threshold selection methods on heavy-tailed distributions with very different tail behaviour. The simulation study is based on the following distributions:

• •

the Student-t distribution with 6 degrees of freedom, which corresponds to $\gamma=1/6$ and $\rho=-1/3$,

• •

the Fréchet distribution with parameter $\alpha=2$ and distribution function $F(x)=\exp(-x^{-\alpha})$ for $x>0$, which implies $\gamma=1/2$ and $\rho=-1$,

• •

the standard Cauchy distribution leading to a tail behaviour with $\gamma=1$ and $\rho=-2$,

• •

the Loggamma distribution with $\gamma=1$ and $\rho=0$ and density function

[TABLE]

• •

the Burr distribution with a parametrisation such that $\gamma=2$, $\rho=-1$ and distribution function

[TABLE]

• •

a logarithmically perturbed Pareto distribution of the random variable $g(U)$ with $\gamma=1$ and $\rho=-1$, where $U\sim\mathrm{Unif}(0,1)$ and $g(x)=x^{-1}/\log(x^{-1})$. This distribution is denoted as negBias due to its negative bias in the Hill estimator.

On these distributions we evaluate the methods by their root mean square error (RMSE) when adaptively estimating $\gamma$ with the Hill estimator relative to the RMSE obtained using $k_{\mathrm{opt}}$,

[TABLE]

where $\mathbb{E}_{n}$ denotes the empirical expectation. These efficiency quotients are also used by, e.g., Guillou and Hall (2001), Gomes et al. (2000) and Drees and Kaufmann (1998). The smaller the quotient the better the threshold selection procedure performs compared to the asymptotically optimal sample fraction $k_{\mathrm{opt}}$. Furthermore, we study the efficiency in quantile estimation with the estimator defined in (10) for $p=0.001$,

[TABLE]

Since we do not know the true minimizer $k_{\mathrm{opt}}$ of the AMSE, we utilize an empirical version suggested by Gomes et al. (2000). Following their approach we approximate $k_{\mathrm{opt}}$ by the mean of 20 independent replicates of $\bar{k}_{\mathrm{opt}}$, which is the minimizer of the empirical MSE based on $1000$ samples, i.e. $\bar{k}_{\mathrm{opt}}=\underset{k}{\mathrm{argmin}}\ \mathbb{E}_{n=1000}[(\hat{\gamma}_{k}-\gamma)^{2}].$

We compare these efficiency values for eight different threshold selection methods. Most of the considered approaches are constructed for adaptive estimation of $\gamma$ applying the Hill estimator. This includes one procedure that looks for a stable region among the Hill estimates, while the others aim to estimate $k_{\mathrm{opt}}$. The only exception is the IHS approach discussed in Section 2, which is motivated to minimize the deviation from the exponential approximation. We still evaluate the performance of this procedure in the same simulations, although it is not primarily tailored for the specific applications. In total, the following methods are considered:

sIHS:

IHS smoothed by using the eBsc package, see Section 2,

SAM:

SAMSEE procedure with $\rho=-1$ as defined by (19) in Section 3,

GH:

method by Guillou and Hall (2001) utilizing $c_{\mathrm{crit}}=1.25$ and $p=1$,

DK:

procedure by Drees and Kaufmann (1998) with fixed $\rho=-1$,

GO:

approach by Goegebeur et al. (2008) defined in their equation (3.3) with fixed $\rho=-1$,

DB:

double bootstrap approach by Danielsson et al. (2001) with the choice $n_{1}=120$ if $n=500$ and $n_{1}=1000$ if $n=5000$,

B:

method by Beirlant et al. (2002) with $\rho=-1$,

RT:

method by Reiss and Thomas (2007) with $\beta=0$ as suggested by Neves and Fraga Alves (2004).

When looking at the results for estimating $\gamma$ adaptively for $n=500$ and $n=5000$ in Table 2, we observe a very diverse picture of methods performing best. Overall we get the impression that SAMSEE together with the approach by Goegebeur et al. (2008) performs most stable over the variety of distributions. This is interesting, because those are the methods which depend least on tuning parameters. The performance of the approaches GH, DK and B is comparable, but we obtain from Table 2 that on average over all distributions the SAMSEE procedure is superior.

For estimating a high quantile SAMSEE also performs convincingly, see Table 3, but additionally sIHS and the approach by Danielsson et al. (2001) show very good efficiency values. They are closely followed by B and GO. Looking at the average performance over all distributions, SAMSEE performs best again. However, if we exclude the negBias distribution, sIHS works superior on average.

In conclusion, we can see that SAMSEE performs very efficiently and comparable to $k_{\mathrm{opt}}$ over all exemplary distributions. It works especially well for estimating a high quantile. Only in the case of estimating $\gamma$ for the Cauchy distribution it performs worse than DK and GO, but still better than most other approaches. Recalling the results of the simulation on the influence of $\rho$ in Table 1 in Section 3.1, it is not very surprising that SAMSEE performs slightly weaker in this situation. There, the Cauchy distribution is the only example we considered that benefits from estimating $\rho$ instead of fixing it to $-1$.

From Table 3 we furthermore observe that sIHS is a strong choice when estimating high quantiles from small samples with $n$ up to $5000$. However, the performance when estimating $\gamma$ is quite variable and it seems that sIHS does not perform particularly well for distributions with a small second order parameter ($\rho\leq-1$). This behaviour is already discussed in Section 2.1 and highlighted in Figure 2: sIHS selects smaller $k$ than optimal for the Hill estimator, especially if $\rho$ is in the regime between $-1$ and $-8$.

The reason why some approaches perform worse on quantiles than they do on $\gamma$ is that the estimator $\hat{q}_{k}(p)$ defined in (10) depends on $\hat{\gamma}_{k}$ in the exponent and is thus very sensitive to overestimation in case of $\gamma>1$. Hence, when estimating a high quantile, an estimate $\hat{\gamma}_{k}$ that is too large will lead to an even stronger overestimation of the quantile. This is why a few outliers among the $\gamma$ estimates can already cause much higher $\mathrm{EFF}_{q}$ values.

5 Application to varying extreme value index

In this section we analyse our new procedures in a financial application, where we study operational losses of a bank. We are, of course, particularly interested in the distributional properties of very high losses. It has been discussed before that it is reasonable to assume the distribution of such extreme losses being heavy-tailed (Chavez-Demoulin et al., 2016; Moscadelli, 2004) and to change with the financial market over time (Hambuckers et al., 2018; Cope et al., 2012). In this context, we want to estimate how the extreme value index changes depending on the univariate covariate time. For this task, we utilize the approaches presented in Sections 2 and 3 for locally optimal selection of a threshold.

The observations of interest are operational losses from the Italian bank UniCredit from 2005 to 2014. In Hambuckers et al. (2018) the data is analysed in a regularized generalized Pareto regression approach including several firm-specific, macroeconomic and financial indicators as covariates. This approach describes the dependence of the GPD parameters on various covariates via parametric functions.

We consider an easier and more direct approach to study the temporal dependence of the extreme value index without taking into account possible interference by other covariates. Our aim is to estimate the time dependent extreme value index $\gamma(t)$ non-parametrically with a simple ad hoc estimator that extends the estimator from de Haan and Zhou (2017) by employing our threshold selection procedures sIHS and SAMSEE. We present the estimator in Section 5.1 and the results we obtain when applying this estimator to the dataset of operational losses in Section 5.2.

5.1 Estimating a varying extreme value index

In de Haan and Zhou (2017), the authors already discussed estimating a trend in the extreme value index non-parametrically. They consider $n$ independent random variables $X_{i}\sim F_{i/n}$, where $F_{s}\in\mathrm{DoA}(G_{\gamma(s)})$ for $s\in[0,1]$. To address this problem, they introduce the following estimator for $\gamma(s)$, which locally applies the Hill estimator and is based on a global sample fraction $k$,

[TABLE]

where $I_{n}(s)$ is the $h$-neighbourhood of $s$, i.e. $I_{n}(s):=\{i:|i/n-s|\leq h\}$. This estimator depends on the choice of the bandwidth $h$ and the global sample fraction $k$, which is then rescaled to $2kh$ for the individual regions $I_{n}(s)$. A small bandwidth $h$ leads to very high variability in $\hat{\gamma}_{k}(s)$ and a large value of $h$ smooths out all interesting features. Thus, the choice of $h$ should balance these two effects.

We suggest a modification of their estimator, where we locally estimate an optimal threshold $\hat{k}(s)$, i.e.

[TABLE]

To compare these two approaches, we repeat the simulation presented in Figure 2 (i) in de Haan and Zhou (2017) on samples of size $n=5000$ with $X_{i}\sim\text{Fr\'{e}chet}(1/\gamma(i/n))$ and $\gamma(s)=1+s$. Figure 6 illustrates the benefits of locally optimizing the threshold via SAMSEE from Section 3, as it strongly tightens the empirical confidence interval around the average, which is obtained from the $2.5$% and $97.5$% quantiles among $1000$ estimates.

5.2 Functional extreme value index of operational losses

The operational losses in the dataset of UniCredit are grouped by the type of event that caused the specific loss. We consider the event type CPBP, which provides sufficient observations for our local estimation approach. The CPBP losses are caused by clients, products and business practices related to derivatives or other financial instruments.

First we want to test if the extreme value index is constant over time. Using the test T4 from Einmahl et al. (2016), we can reject the null hypotheses with a $p$-value that is virtually zero and thus are confident that the extreme value index of the losses is indeed varying over time.

We apply the new methodology from (23) to these losses via estimating $\hat{k}$ with sIHS from Section 2 and the SAMSEE approach from Section 3. Figure 7 shows the estimates we obtain for the event type CPBP. It is clearly visible that both procedures yield similar estimates for most time points and that the simple ad hoc estimators recover an increase of the severity of high losses during the financial and Euro crisis from 2008 to 2011. A similar overall trend in the extreme value index can also be identified in the estimates of Hambuckers et al. (2018) for CPBP.

For a more extensive discussion of the data and results of the more complex model including further covariates we refer to Hambuckers et al. (2018).

Appendix A Theoretical results and proof of Theorem 3

Lemma 1** (Distribution of the Hill estimator).**

Let $X_{1},\dots,X_{n}\overset{i.i.d.}{\sim}F$ for $F\in\mathrm{DoA}(G_{\gamma})$ with $\gamma>0$. Then the following distributional representation for the Hill estimator holds,

[TABLE]

where $G_{k}\sim\Gamma(k,\gamma/k)$ and for $b_{n,k}$ it holds that

[TABLE]

for some $b_{c}\in\mathbb{R}$.

Proof.

From the first order condition

[TABLE]

follows that $U=F^{\leftharpoonup}(1-1/x)$ is regularly varying with index $\gamma$ and there exists a slowly varying function $\ell_{U}$, such that $U(x)=x^{\gamma}\ell_{U}(x)$. Let $P_{1},P_{2},\dots$ be i.i.d. random variables with distribution function $1-1/y$. Note that $U(P_{i})\overset{\mathcal{D}}{=}X_{i}$. We define $Y_{(k-i,k)}:=\log(X_{(n-i,n)})-\log(X_{(n-k,n)})$, for which it follows that

[TABLE]

Note that $\log(P_{i})$ is standard exponentially distributed. By Lemma 3.2.3 in de Haan and Ferreira (2006) follows for i.i.d. standard exponential random variables $E_{1},E_{2},\dots$ that $\left\{E_{(n-i,n)}-E_{(n-k,n)}\right\}_{i=1}^{k-1}\overset{\mathcal{D}}{=}\left\{E_{(k-i,k)}\right\}_{i=0}^{k-1}.$

Hence, we obtain for the Hill estimator that

[TABLE]

where $G_{k}\sim\Gamma(k,\gamma/k)$ as the sum of i.i.d. exponentials and $b_{n,k}$ denotes the second average.

If $k/n\rightarrow 0$, $P_{(n-k,n)}\rightarrow\infty$ almost surely by Lemma 3.2.1 in de Haan and Ferreira (2006). Since $\ell_{U}$ is slowly varying, $b_{n,k}$ converges to zero almost surely.

If $k/n\rightarrow c\in(0,1]$, $P_{(n-k,n)}\rightarrow 1/c$ in probability by Cor. 2.2.2 in de Haan and Ferreira (2006). Thus, by the weak law of large numbers

[TABLE]

∎

Lemma 2**.**

Let $X_{1},\dots,X_{n}\overset{i.i.d.}{\sim}F$ for $F\in\mathrm{DoA}(G_{\gamma})$ with $\gamma>0$. Then the following holds for $\mathrm{SKK}=(4-k)/(2\hat{\gamma}_{k}k)$ and $\mathrm{SKK}_{-}=(4+k)/(2\hat{\gamma}_{k}k)$ depending on the sample fraction $k$.

If k is finite,

• •

then $\mathbb{E}[\mathrm{SKK}]+(2\gamma)^{-1}\rightarrow\frac{3}{2\gamma(k-1)}>0$, as $n\rightarrow\infty$.

• •

then $\mathbb{E}[\mathrm{SKK}^{-}]-(2\gamma)^{-1}\rightarrow\frac{5}{2\gamma(k-1)}>0$, as $n\rightarrow\infty$. 2. 2.

If $k\rightarrow\infty,\ k/n\rightarrow 0$,

• •

then $\mathbb{E}[\mathrm{SKK}]+(2\gamma)^{-1}\rightarrow 0$, as $n\rightarrow\infty$.

• •

then $\mathbb{E}[\mathrm{SKK}^{-}]-(2\gamma)^{-1}\rightarrow 0$, as $n\rightarrow\infty$. 3. 3.

If $k\rightarrow\infty,\ k/n\rightarrow c>0$

• •

and $b_{c}\geq 0$, then $\mathbb{E}[\mathrm{SKK}]+(2\gamma)^{-1}\rightarrow\frac{b_{c}}{2\gamma(\gamma+b_{c})}>0$, as $n\rightarrow\infty$.

• •

and $b_{c}\leq 0$, then $\mathbb{E}[\mathrm{SKK}^{-}]-(2\gamma)^{-1}\rightarrow\frac{-b_{c}}{2\gamma(\gamma+b_{c})}>0$, as $n\rightarrow\infty$.

Proof.

We proof these three alternative statements to obtain that the minimizing sequence $k=k_{n}$ is an intermediate sequence.

Let $k$ be finite, by Lemma 1 holds that $\hat{\gamma}_{k}\overset{\mathcal{D}}{=}G_{k}+b_{k,n}$, where $G_{k}\sim\Gamma(k,\gamma/k)$ and $b_{k,n}\rightarrow 0$ as $n\rightarrow\infty$. Thus, as $n\rightarrow\infty$ it follows that

[TABLE] 2. 2.

The statement follows from the consistency of the Hill estimator and the continuous mapping theorem. 3. 3.

Let $k\rightarrow\infty$ and $k/n\rightarrow c>0$, then it holds by Lemma 1 that $\hat{\gamma}_{k}\overset{\mathcal{D}}{=}G_{k}+b_{k,n}$, where $G_{k}\overset{\mathbb{P}}{\rightarrow}\gamma$ and $b_{k,n}\rightarrow b_{c}$ as $n\rightarrow\infty$. Thus, as $n\rightarrow\infty$

[TABLE]

∎

Lemma 3**.**

Let $E_{1},\dots,E_{n}$ be i.i.d. standard exponential random variables and $k$ s.t. $1\leq k\leq n$ and $k\rightarrow\infty$ as $n\rightarrow\infty$. We define the following random variables,

[TABLE]

where $e_{i}^{k}:=\sum_{l=i}^{k}l^{-1}=\mathbb{E}[E_{(k-i+1,k)}]$. Then it holds for $n\rightarrow\infty$ that

[TABLE]

Proof.

We use the Cramér-Wold device that gives us a joint normal limit distribution if all linear combinations have an univariate normal limit distribution. For $a_{1},a_{2},a_{3}\in\mathbb{R}$ we study

[TABLE]

To prove asymptotic normality for the sum in (24) we use Liapounov’s central limit theorem (CLT), see Theorem 7.1.2. in Chung (1974). We consider a sum $S_{n}:=\sum_{i=1}^{k}X_{i,k}$ of independent random variables fulfilling the following three conditions,

$\mathbb{E}[X_{i,k}]=0,\ \forall k\ \forall i$, 2. 2)

$\sum_{i=1}^{k}\mathrm{Var}(X_{i,k})=\sigma^{2}$, 3. 3)

$\Gamma(k)=\sum_{i=1}^{k}\mathbb{E}[|X_{i,k}|^{3}]\longrightarrow 0,\ \text{as }k\rightarrow\infty,$

Then the CLT proves a standard normal limit for $S_{n}$. We define

[TABLE]

for $i=1,\dots,k$ where $e_{k+1}^{k}:=0$, such that $\sum_{i=1}^{k}X_{i,k}\approx$(24), where $c_{k}\approx c$ if $c_{k}\rightarrow c$ as $k\rightarrow\infty$ and the approximation error is due to

[TABLE]

Now we have to check the three conditions. Condition 1) follows immediately from $\mathbb{E}[E_{i}]=1$ and $\mathbb{E}[E_{i}^{2}]=2$. For condition 2) we need to calculate the variance

[TABLE]

since $\mathrm{Var}(E_{i})=1$, $\mathrm{Var}(E_{i}^{2})=20$ and $\operatorname{Cov}(E_{i},E_{i}^{2})=\mathbb{E}[E_{i}^{3}]-\mathbb{E}[E_{i}^{2}]\mathbb{E}[E_{i}]=4$. With the approximation

[TABLE]

follows that

[TABLE]

Condition 3) holds with

[TABLE]

as $k\rightarrow\infty$ and for a constant $c>0$, since the exponential distribution has finite moments and $\sum_{i=1}^{k}(e_{i+1}^{k})^{3}/k\approx 6$. Thus, we obtain that

[TABLE]

This is the limiting distribution of the sum in (24) and also follows from the joint normal distribution. ∎

Lemma 4**.**

Let $X_{1},\dots,X_{n}\overset{i.i.d.}{\sim}F$ for $F\in\mathrm{DoA}(G_{\gamma})$ with $\gamma>0$ and $P_{1},P_{2},\dots$ be i.i.d. random variables with distribution function $1-1/y$. We define

[TABLE]

If the second order condition

[TABLE]

holds for $\rho<0$ and $x>0$, then

[TABLE]

Proof.

The results in (26) and (27) are already stated in the proof of Theorem 1 in de Haan and Peng (1998).

To prove (28) we follow the proof of the asymptotic normality of the Hill estimator in de Haan and Ferreira (2006). Let $A_{0}$ be such that $A(t)/A_{0}(t)\rightarrow 1$, as $t\rightarrow\infty$. Then, for each $\epsilon>0$ there exists a $t_{0}$ such that for $t\geq t_{0}$ and $x\geq 1$ the inequality in Theorem B.2.18 in de Haan and Ferreira (2006) holds. For $t=P_{n-k,n}$ and $x=P_{(n-i,n)}/P_{(n-k,n)}$ we obtain that

[TABLE]

The second term can be approximated by

[TABLE]

and for the third term holds

[TABLE]

as $k\rightarrow\infty$. Note that for $E_{1},\dots,E_{n}$ i.i.d. standard exponential random variables follows by Rényi’s representation that

[TABLE]

This distributional equality enables the following transformations,

[TABLE]

Thus,

[TABLE]

Combining the above arguments as in the proof of Theorem 3.2.5 in de Haan and Ferreira (2006) gives (28). ∎

Lemma 5**.**

For $k\rightarrow\infty$, $k/n\rightarrow 0$ and $k/K\rightarrow c$ with $0<c<1$,

[TABLE]

where $R_{k,n}$ and $P_{k,n}$ are defined in Lemma 3.

Proof.

Let $E_{1},E_{2},\dots$ be i.i.d. standard exponential random variables, where $\operatorname{Cov}(E_{i},E_{j})$ is equal to 1 if $i=j$ and 0 otherwise. Then

[TABLE]

where $c_{k}\approx c$ again denotes that $c_{k}\rightarrow c$ as $k\rightarrow\infty$.

In the same way we obtain

[TABLE]

∎

Theorem 4**.**

Let $\bar{\gamma}_{k}:=\frac{1}{k}\sum_{i=1}^{k}\hat{\gamma}_{i}$ denote the average over Hill estimates. Further let $X_{1},\dots,X_{n}$ be i.i.d. random variables with distribution function $F\in\mathrm{DoA}(G_{\gamma})$, $\gamma>0$. If

[TABLE]

with $U(x):=F^{\leftharpoonup}\left(1-\frac{1}{x}\right)$ holds and $k\rightarrow\infty$ and $k/n\rightarrow 0$ as $n\rightarrow\infty$,

[TABLE]

with $\lambda:=\underset{k\rightarrow\infty}{\lim}\sqrt{k}A(n/k)$.

Proof.

First we have to rewrite the average over the Hill estimator,

[TABLE]

where $\overline{YE}$ is defined in Lemma 4. Following the proof of Lemma 4 it holds that the last term above is in distribution equal to

[TABLE]

From Corollary 2.2.2 in de Haan and Ferreira (2006) follows that

[TABLE]

Thus, \big{(}\log(X_{(n,n)})-\log(X_{(n-k,n)})\big{)}/k=O_{p}(\log(k)/k)+O_{p}(A(n/k)/k), and by Lemma 4 follows

[TABLE]

∎

Theorem 5**.**

Let $\bar{\gamma}_{\mathrm{u}p,K,k}:=\frac{1}{K-k+1}\sum_{i=k}^{K}\hat{\gamma}_{i}$ be the upper mean and $k$ and $K$ intermediate sequences, i.e. $k\rightarrow\infty$ and $k/n\rightarrow 0$, as $n\rightarrow\infty$. Further, let $\sqrt{k}A(n/k)\rightarrow\lambda$ and $k/K\rightarrow c$ with $0<c<1$. Under the conditions of Theorem 4, it holds that

[TABLE]

Proof.

We can write the upper mean as a combination of two averaged Hill estimators and apply Theorem 4,

[TABLE]

We approximate $k/K$ by $c$ and obtain

[TABLE]

Now we need the covariance between $R_{n,k}$ and $R_{n,K}$, see Lemma 5, and apply the following property of regular varying functions,

[TABLE]

This leads to

[TABLE]

∎

Proof of Theorem 3.

The bias estimator is defined as $\bar{b}_{\mathrm{up},K,k}=\bar{\gamma}_{\mathrm{up},K,k}-\bar{\gamma}_{K}$ in equation (14). Thus, we can utilize the asymptotic normality results for $\bar{\gamma}_{k}$ and $\bar{\gamma}_{\mathrm{up},K,k}$ in Theorem 4 and 5. Following the proofs of these theorems it holds that

[TABLE]

Here the random variable $R_{k,n}$ is defined in Lemma 3 and we know that $R_{k,n}$ has a normal limit distribution. With Lemma 3 and Lemma 5 we obtain the following variance,

[TABLE]

The bias term of the normal limit is

[TABLE]

as $n\rightarrow\infty$, which follows from the regular variation of $A$ and due to $\sqrt{k}A(n/k)\rightarrow\lambda$. Since

[TABLE]

the statement of the theorem follows immediately. ∎

Acknowledgement

Support of the DFG RTG 2088 (B4) is gratefully acknowledged. We are grateful to Julien Hambuckers for providing us the dataset of operational losses from UniCredit.