2.1. Proper scoring rules

Let $\mathcal{Y}=\mathbb{R}$ be equipped with the Borel $\sigma$ -algebra $\mathcal{A}$ , and let $\mathcal{F}$ denote a convex class of distribution functions on $\mathcal{Y}$ . Before defining the score function, we introduce the term quasi-integrable.

With this setup, a scoring rule is defined as

With this notation, we write

(1) \begin{align} S(F,G)\,:\!=\,\mathbb{E}[S(F,{Y})]=\int S(F,y)\mathop{}\!\mathrm{d} G(y) \end{align}

where $F,G\in \mathcal{F}$ , and Y is a random variable with distribution function G.

A natural requirement of a scoring rule is that it is proper, which has the following definition.

Definition 2.3 (Proper scoring rule). Scoring rule S is proper relative to $\mathcal{F}$ if

\begin{align*} S(G,G)\geq S(F,G) \end{align*}

for all $F,G\in \mathcal{F}$ . It is strictly proper if the inequality holds with equality if and only if $F=G$ .

Accordingly, using proper scoring rules incentivizes the forecaster to report their honest belief about the distribution. Therefore, the forecaster has no incentive to distort or strategically alter their forecast. Definition 2.3 is stated under the positive score orientation, meaning that larger score values correspond to better predictive performance.

Two widely used scoring rules are the logarithmic score (LogS) and the continuous ranked probability score (CRPS). The former is defined as

\begin{align*}\mathrm{LogS}(F,y)\,:\!=\, \log {f}(y)\ \end{align*}

where f denotes the density of F. CRPS is defined as

\begin{align*} CRPS(F,y)\,:\!=\,-\int_{-\infty}^{\infty}\left(F(x)-\mathbf{1}\{y\le x\}\right)^2\,dx.\end{align*}

CRPS can be viewed as a special case of the Energy score and, unlike the logarithmic score, it does not require the existence of a density. More generally, the Energy score with parameter $\beta \in (0,2)$ is defined as

\begin{align*}ES_\beta(F,y)\,:\!=\,\frac12 \mathbb{E}|X-X'|^\beta - \mathbb{E}|X-y|^\beta,\end{align*}

where $X,X' \sim F$ are independent random variables. The CRPS coincides with $\mathrm{ES}_\beta$ for $\beta = 1$ . In this paper, we present results for the Energy score, of which the CRPS is a special case.

Both the logarithmic score and the Energy score with $\beta\in(0,2)$ are strictly proper; see, for example, Gneiting and Raftery (Reference Gneiting and Raftery2007).

2.2. Extreme value theory framework

A central result in extreme value theory is the Fisher–Tippett–Gnedenko theorem, which states that suitably normalized maxima converge in distribution to one of three possible types: Gumbel, Fréchet, or Weibull. These domains of attraction correspond to fundamentally different tail behaviors in the underlying distribution. See De Haan and Ferreira (Reference De Haan and Ferreira2006) for more information.

The Fréchet domain, also referred to as the heavy-tailed domain, encompasses distributions whose survival function decays according to a power law. Such distributions have infinite right endpoints and exhibit polynomially decreasing tails, implying that extreme events occur with non-negligible probability. Formally, a distribution function G is said to belong to the Fréchet domain of attraction if its survival function is regularly varying. That is, there exists a constant $\gamma \gt 0$ such that

(2) \begin{align}\lim_{t \to \infty} \frac{1 - G(tx)}{1 - G(t)} = x^{-1/\gamma}, \quad x \gt 0.\end{align}

The parameter $\gamma$ , called the tail index, quantifies the heaviness of the tail. When (2) holds, we say that the survival function $1-G$ is regularly varying with index $-1/\gamma$ , and write $1-G \in \mathcal{R}_{-1/\gamma}$ .

An equivalent characterization can be expressed in terms of the tail quantile function

\begin{align*}U(x) = \inf\{ y \in \mathbb{R} \,:\, G(y) \ge 1 - 1/x \}, \qquad x \gt 1.\end{align*}

If $1-G \in \mathcal{R}_{-1/\gamma}$ , then the associated tail quantile function U is regularly varying with index $\gamma$ , that is, $U\in \mathcal{R}_{\gamma}$

In this paper, we let $(Y_1,Y_2, \ldots , Y_n)$ be an i.i.d. sample with common distribution function G. Throughout, we assume that $1-G\in \mathcal{R}_{-1/\gamma_G}$ . This implies that Y is in the Fréchet domain with tail index $\gamma_G$ . We now introduce the tail counterpart of G, which is given by

\begin{align*}G^t(x)=\frac{G(tx)-G(t)}{1-G(t)}\quad \text{for $t\geq 1$.}\end{align*}

We denote $G^\circ$ as the limit distribution of $G^t$ and we call it the tail distribution. $Y^\circ$ denotes a random variable with distribution $G^\circ$ . As Y is in the Fréchet domain, we obtain

\begin{align*}1-G^\circ(x)=\lim_{t\to \infty}1-G^t(x)=x^{-1/\gamma_G},\end{align*}

i.e., $G^t$ converges to a Pareto distribution with parameter $1/\gamma_G$ .

4.1. Scoring with normalized upper order statistics

We begin with consistency. Consistency and asymptotic normality follow directly from Segers (Reference Segers2001) (Theorems 2.1 and 4.5). Here, the assumptions and theorem are restated in our notation. Although formulated below in terms of scoring rules, these results also apply to scoring functions $s(x,\cdot)$ , for any fixed action x, provided the stated regularity conditions are satisfied (see Gneiting, Reference Gneiting2011 for further discussion of scoring functions).

This assumption is mild and does not impose a substantial restriction, since commonly used scoring rules are typically continuous. The domination condition provides the necessary tail control.

Theorem 4.2. Let Assumption 4.1 hold. In addition, assume that $k,n \to \infty$ and $k/n \to 0$ . Then,

\begin{align*}S_k(F)\;\rightarrow\;\mathbb{E}\!\left[S(F, Y^\circ)\right]\end{align*}

in probability, where $Y^\circ$ has distribution $G^\circ$ .

In practical terms, this result ensures that empirical tail scores converge to their population targets when k is chosen in the standard EVT regime.

We next turn to asymptotic normality, which follows as a direct consequence of Theorem 4.5 in Segers (Reference Segers2001). To establish this result, it is necessary to control the growth rate of k. We therefore introduce the following notation. Let $\xi \in \mathbb{R}$ , define the function $h_\xi \,:\, (0,\infty) \to \mathbb{R}$ by

\begin{align*}h_\xi(x)= \int_1^x y^{\xi-1}\,dy=\begin{cases}\displaystyle \frac{x^\xi - 1}{\xi}, & \text{if } \xi \neq 0, \\[1ex]\log x, & \text{if } \xi = 0.\end{cases}\end{align*}

We may now restate Theorem 4.5 of Segers (Reference Segers2001) in our notation.

Theorem 4.3. Let G be a distribution function with tail quantile function U satisfying

(4) \begin{equation}\lim_{t \to \infty} \frac{U(tx)/U(t) - x^{\gamma_G}}{a(t)}= \pm x^{\gamma_G} h_{\xi}(x), \qquad x \ge 1,\end{equation}

for some $\gamma_G \gt 0$ , $\xi \le 0$ and $a \in \mathcal{R}_{\xi}$ satisfying $a(t) \to 0$ as $t \to \infty$ . Let $F\in\mathcal{F}$ be a fixed candidate distribution used in the score comparison. Let $S \,:\,\mathcal{F} \times [1,\infty) \to \mathbb{R}$ be an absolutely continuous function for $F\in \mathcal{F}$ with $\frac{\partial S}{\partial y}$ such that

\begin{align*}\bigg|\frac{\partial S}{\partial y}\bigg| \le A y^{\rho-1}, \qquad y \ge 1,\end{align*}

for some $A \gt 0$ and some $\rho \lt (2\gamma_G)^{-1}$ . If k satisfies $\sqrt{k}\,a{\left(\frac{n+1}{k+1}\right)} \to 0$ when $k,n \to \infty$ and $k/n \to 0$ , then

(5) \begin{equation} \sqrt{k}\bigl(S_k(F) - \mathbb{E}[S(F,Y^\circ)]\bigr)\xrightarrow{d} \mathcal{N}\!\left(0,\operatorname{Var}(S(F,Y^\circ))\right).\end{equation}

The derivative condition is a regularity condition ensuring that the scoring rule behaves smoothly in the tail. As $y\mapsto S(F,y)$ is assumed absolute continuous, the condition also ensures that the variance exists. This normal approximation provides a direct basis for uncertainty quantification when comparing score levels across competing tail models.

The following corollary shows under which conditions Theorems 4.2 and 4.3 are applicable for LogS and ES.

Proof. We verify the conditions separately for LogS and $\mathrm{ES}_\beta$ .

First consider LogS. If we denote the density of $F_\gamma$ as $f_\gamma$ , we have

\begin{align*} \mathrm{LogS}(F_\gamma,y)=\log\!(\kern1.5pt f_\gamma(y))=\log(1/\gamma)-\left(\frac{1}{\gamma}+1\right)\log(y). \end{align*}

$\mathrm{LogS}(F_\gamma,y)$ is absolute continuous for $y\gt1$ and for any $\delta\in (0,1)$ , it is dominated by $Ay^{(1-\delta)/\gamma_G}$ , where A is some constant. Further, we have

\begin{align*} \frac{\partial \mathrm{LogS}}{\partial y}=-\left(\frac{1}{\gamma}+1\right)y^{-1} \end{align*}

For any $\rho \in (0,(2\gamma_G)^{-1})$ , this derivative is dominated by $Ay^{\rho-1}$ .

Next consider $\mathrm{ES}_\beta$ :

\begin{align*}\mathrm{ES}_\beta(F_\gamma,y)&= {\frac{1}{2} \int_1^\infty\int_1^\infty |x - x'|^\beta \, \mathop{}\!\mathrm{d} F_\gamma(x) \mathop{}\!\mathrm{d} F_\gamma(x')-\int_1^\infty |x - y|^\beta \, \mathop{}\!\mathrm{d} F_\gamma(x)}.\end{align*}

For $y\geq 1$ , this expression is absolutely continuous. We first bound the first integral. For $\beta\lt\frac{1}{\gamma}$ , we have that

(6) \begin{align} \int_1^\infty |x - y|^\beta \, dF(x)&\leq A_1 \left(\int_1^y(y-1)^\beta x^{-1/\gamma -1}\mathop{}\!\mathrm{d}\kern1pt x + \int_y^\infty x^{\beta-1/\gamma-1}\mathop{}\!\mathrm{d}\kern1pt x\right) \\ &{\leq A_2 (y^{\beta}+y^{\beta-1/\gamma-1})}{\,\le A_3y^\beta}\nonumber\end{align}

for some constants $A_1,A_2,{A_3\gt0}$ . Hence, we can conclude that (6) is finite. Furthermore, for $\beta\lt\frac{1}{\gamma_G}$ , there exist some $\delta\in(0,1)$ and A such that it is dominated by $Ay^{(1-\delta)/\gamma_G}$ . Next, for the second integral, it is enough to verify finiteness since it does not depend on y:

\begin{align*} &\int_1^\infty\int_1^\infty |x - x'|^\beta \, \mathop{}\!\mathrm{d} F_\gamma(x) \mathop{}\!\mathrm{d} F_\gamma(x')\nonumber\leq C_\beta\left(\int_1^\infty x^\beta\mathop{}\!\mathrm{d} F_\gamma(x)+\int_1^\infty (x')^{\beta}\mathop{}\!\mathrm{d} F_\gamma(x') \right)\end{align*}

for some $\beta$ dependent constant $C_\beta$ . Each of these integrals is finite for $\beta\lt\frac{1}{\gamma}$ . Hence, we can conclude that $\mathrm{ES}_\beta$ satisfies Assumption 4.1 for $\beta\lt\min\left\{1/\gamma,1/\gamma_G\right\}$ . Lastly, we have

\begin{align*} \left|\frac{\partial \mathrm{ES}_\beta}{\partial y}\right|\leq {\beta \int_1^\infty |x-y|^{\beta-1}\mathop{}\!\mathrm{d} F(x)} \leq A_4y^{{\beta-1}}, \end{align*}

where $A_4$ is a constant. If $\beta \lt (2\gamma_G)^{-1}$ , choose $\rho$ with $\beta \lt \rho \lt (2\gamma_G)^{-1}$ . Since $y^{\beta-1} \leq y^{\rho-1}$ for $y \geq 1$ , the derivative condition is satisfied.

Therefore, both scoring rules satisfy the required conditions, which proves the claim.

This shows that the derivative condition mainly affects $\mathrm{ES}_\beta$ , where it imposes an additional restriction on the score parameter $\beta$ . For LogS, the condition is automatically satisfied and therefore imposes no further restriction.

Taken together, Theorems 4.2–4.3 and Corollary 4.4 justify the following workflow: compute empirical tail scores over a range of (k), compare candidates by the score level, and use the normal approximation to assess whether observed score differences are practically meaningful.

4.2. Tail index estimation via score optimization

Beyond their use for ranking predictive distributions, we also investigate whether scoring rules can be employed for estimating the tail index. In this section, we assume that the predictive distribution is Pareto, denoted by $F_\gamma$ , where $\gamma$ is the candidate tail index. Let $\Gamma \,:\!=\, [\gamma_L,\gamma_U]$ , with $0\lt\gamma_L\lt\gamma_G\lt\gamma_U\lt\infty$ . For shorthand, define

\begin{align*}S_k(\gamma)\,:\!=\,\frac{1}{k}\sum_{i=1}^k S\left(F_\gamma, \frac{Y_{n,n-i+1}}{Y_{n,n-k}}\right).\end{align*}

Restricting the predictive distributions to the Pareto family also gives a scoring-function interpretation of the procedure. In particular, any scoring rule S induces a scoring function for the tail index through

\begin{align*}s(\gamma,y)\,:\!=\,S(F_\gamma,y).\end{align*}

Thus, the estimation procedure below may equivalently be viewed as optimization of an induced scoring function over $\gamma$ . Moreover, strict propriety of S implies that the expected induced score is uniquely maximized at the true tail index within the Pareto family.

The following theorem establishes consistency of the resulting tail-index estimator.

Theorem 4.7. Assume that S is continuous on $\Gamma\times[1,\infty)$ and a strictly proper scoring rule. Further assume Assumption 4.1 and that for $\gamma\in \Gamma$ , there exists a measurable function $D\,:\,[1,\infty)\to[0,\infty)$ such that

(7) \begin{align}\sup_{\gamma\in\Gamma}|S(F_\gamma,z)|\le D(z),\qquad z\ge 1,\end{align}

with $D(z)\le A z^{(1-\delta)/\gamma_G}$ for some $A\gt0$ and $\delta\in(0,1)$ . In particular, $E[|S(F_\gamma,Y^\circ)|]\lt\infty$ for all $\gamma\in\Gamma$ . Let $k\to\infty$ with $n/k\to\infty$ . Then,

(8) \begin{align} \hat{\gamma}_k(S)=\mathrm{argmax}_{\gamma\in \Gamma} \frac{1}{k}\sum_{i=1}^k S\left(F_\gamma, \frac{Y_{n,n-i+1}}{Y_{n,n-k}}\right) \xrightarrow{{\mathbb{P}}} \gamma_G. \end{align}

We first establish the following uniform convergence lemma.

Lemma 4.8. Assume the same conditions as in Theorem 4.7. Then, as $k\to\infty$ and $n/k\to\infty$ ,

\begin{align*}\sup_{\gamma\in\Gamma}\left|S_k(\gamma)-E[S(F_\gamma,Y^\circ)]\right|\xrightarrow{P}0 .\end{align*}

Proof of Lemma 4.8. Let

\begin{align*} u(z,\theta,\tau)\,:\!=\, \sup_{\gamma{\in \Gamma}:|\gamma-\theta|\leq \tau}|S(F_\gamma,z)-S(F_\theta,z)|. \end{align*}

First, we argue that $z\mapsto u(z,\theta,\tau)$ is measurable. For fixed $\theta\in\Gamma$ and $\tau\gt0$ , let

\begin{align*} \Gamma_{\theta,\tau}\,:\!=\,\{\gamma\in\Gamma\,:\,|\gamma-\theta|\le \tau\}. \end{align*}

Since $\Gamma_{\theta,\tau}$ is compact and $\gamma\mapsto S(F_\gamma,z)$ is continuous for each fixed z, the supremum over $\Gamma_{\theta,\tau}$ may be taken over any countable dense subset $Q_{\theta,\tau}\subset \Gamma_{\theta,\tau}$ . Hence,

\begin{align*}u(z,\theta,\tau)=\sup_{\gamma\in Q_{\theta,\tau}}|S(F_\gamma,z)-S(F_\theta,z)|.\end{align*}

For each fixed $\gamma$ , the map $z\mapsto S(F_\gamma,z)$ is measurable, since S is a scoring rule. Therefore, $u(\cdot,\theta,\tau)$ is the pointwise supremum of countably many measurable functions and is consequently measurable.

By the continuity assumption of S, we have that $u(Y^\circ,\theta,\tau)\to 0$ for $\tau \to 0$ from above. Since S is dominated by D, by the triangle-inequality, we have that $u(Y^\circ,\theta,\tau)\leq 2D(Y^\circ)$ . Furthermore, since $Y^\circ$ is Pareto with tail index $\gamma_G$ , the bound $D(z)\leq A z^{(1-\delta)/\gamma_G}$ implies $E[D(Y^\circ)]\lt\infty$ . Hence, for a specific $\theta$ , we have by the dominated convergence theorem, $\mathbb{E} [u(Y^\circ,\theta,\tau)]\to 0$ as $\tau\to 0$ .

Now use compactness of the parameter space to cover it by finitely many balls $B(\theta_m,\tau_m)$ , $m=1,\dots,M$ , where $\theta_m\in \Gamma$ . By choosing M sufficiently large, we can ensure that for each ball

(9) \begin{align}\mu_m\,:\!=\,\mathbb{E}\big[u(Y^\circ,\theta_m,\tau_m)\big] \lt \varepsilon \end{align}

for any given $\varepsilon \gt 0$ . For a given $\gamma\in B_m$ we have

(10) \begin{align} |S_k(\gamma)-\mathbb{E}[S(F_\gamma,Y^\circ)]| &\leq |S_k(\gamma)-S_k(\theta_m)| \end{align}

(11) \begin{align} \qquad\qquad\qquad\qquad\qquad\qquad\qquad +|S_k(\theta_m)-\mathbb{E}[S(F_{\theta_m},Y^\circ)]| \end{align}

(12) \begin{align} \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad + |\mathbb{E}[S(F_{\theta_m},Y^\circ)]-\mathbb{E}[S(F_\gamma,Y^\circ)]| . \end{align}

For (10), note that

\begin{align*} |S_k(\gamma)-S_k(\theta_m)| \leq \left|\frac{1}{k}\sum_{i=1}^k u\left(\frac{Y_{n,n-k+i}}{Y_{n,n-k}},\theta_m,\tau_m\right)-\mu_m\right|+\mu_m. \end{align*}

As argued earlier, we have

\begin{align*} u(z,\theta_m,\tau_m) \leq 2Az^{(1-\delta)/\gamma_G}, \end{align*}

and by continuity of S and compactness of $\Gamma_{\theta,\tau}$ , we have that $z\mapsto u(z,\theta,\tau)$ is continuous. Hence, we can use Theorem 4.2 on $z\mapsto u(z,\theta_m,\tau_m)$ to conclude that

\begin{align*} \frac{1}{k}\sum_{i=1}^k u\left(\frac{Y_{n,n-k+i}}{Y_{n,n-k}},\theta_m,\tau_m\right) \to \mu_m, \end{align*}

in probability. Since $\mu_m$ can be made arbitrarily small by selecting a large M, (10) can be made arbitrarily small. Likewise, (11) converges to 0 by Theorem 4.2, while (12) can be made arbitrarily small using

\begin{align*} |\mathbb{E}[S(F_{\theta_m},Y^\circ)]-\mathbb{E}[S(F_\gamma,Y^\circ)]|\leq \mathbb{E}\big[u(Y^\circ,\theta_m,\tau_m)\big] \end{align*}

and (9).

Altogether, this shows that for any $\gamma \in B(\theta_m,\tau_m)$ , the quantity

\begin{align*}|S_k(\gamma)-\mathbb{E}[S(F_\gamma,Y^\circ)]|\end{align*}

can be made arbitrarily small by choosing M sufficiently large. Since this holds for each of the finitely many balls in the cover, it also holds after taking the supremum over $\gamma \in \Gamma$ for M sufficiently large, which concludes the proof.

Proof of Theorem 4.7. By Lemma 4.8, we have

(13) \begin{align} \sup_{\gamma \in \Gamma}|S_k(\gamma)-\mathbb{E}[S(F_\gamma,Y^\circ)]|\xrightarrow{{\mathbb{P}}}0 \end{align}

as $k,n/k\to\infty$ . Using (13), we get that

(14) \begin{align} S_k(\hat{\gamma})\geq S_k(\gamma_G)\geq \mathbb{E}[S(F_{\gamma_G},Y^\circ)]-o_{\mathbb{P}}(1). \end{align}

Using (14) and (13), we get that

(15) \begin{align} \mathbb{E}[S(F_{\gamma_G},Y^\circ)] - \mathbb{E}[S(F_{\hat{\gamma}},Y^\circ)] \leq S_k(\hat{\gamma})+o_{\mathbb{P}}(1)-\mathbb{E}[S(F_{\hat{\gamma}},Y^\circ)] \end{align}

(16) \begin{align} \qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad \leq \sup_{\gamma \in \Gamma} |S_k(\gamma)-\mathbb{E}[S(F_\gamma,Y^\circ)]|+o_{\mathbb{P}}(1) \xrightarrow{{\mathbb{P}}}0 \end{align}

as $k,n/k\to \infty$ . Since $Y^\circ\sim F_{\gamma_G}$ , strict propriety of S implies that

\begin{align*}\mathbb{E}[S(F_\gamma,Y^\circ)]\lt\mathbb{E}[S(F_{\gamma_G},Y^\circ)]\end{align*}

whenever $F_\gamma\neq F_{\gamma_G}$ . By identifiability of the Pareto family, this holds for every $\gamma\neq\gamma_G$ . Moreover, the map $\gamma\mapsto \mathbb{E}[S(F_\gamma,Y^\circ)]$ is continuous by dominated convergence. Hence, for every $\varepsilon\gt0$ , the compact set $\{\gamma\in\Gamma:|\gamma-\gamma_G|\geq\varepsilon\}$ admits a maximizer. Since this maximizer differs from $\gamma_G$ , we obtain

\begin{align*}\sup_{\gamma:|\gamma-\gamma_G|\geq\varepsilon}\mathbb{E}[S(F_\gamma,Y^\circ)]\lt\mathbb{E}[S(F_{\gamma_G},Y^\circ)].\end{align*}

Hence, for every $\varepsilon\gt0$ , there exists an $\eta\gt0$ such that

\begin{align*} \mathbb{E}[S(F_{\hat{\gamma}},Y^\circ)]\lt\mathbb{E}[S(F_{\gamma_G},Y^\circ)]-\eta \quad \text{when } |\hat\gamma-\gamma_G|\geq \varepsilon . \end{align*}

Therefore,

\begin{align*} \{|\hat{\gamma}-\gamma_G|\geq \varepsilon\} \subset\{\mathbb{E}[S(F_{\hat{\gamma}},Y^\circ)]\lt\mathbb{E}[S(F_{\gamma_G},Y^\circ)]-\eta\}, \end{align*}

and as the right-hand side converges to 0 in probability by (16), the left-hand side also converges to 0.

When the empirical maximizer lies on the boundary of $\Gamma$ , the estimate is a boundary optimum induced by the restricted parameter space. It may therefore reflect that the unconstrained empirical optimum lies outside $\Gamma$ , so such a solution should be interpreted with care rather than as direct evidence that the true tail index coincides with the boundary point. For the logarithmic score, this boundary effect can be described explicitly: the unconstrained estimator is the Hill estimator, while the constrained estimator over $\Gamma$ is its projection onto $\Gamma$ .

Remark 4.9. Let $F_\gamma^\circ$ denote the Pareto distribution on $[1,\infty)$ with tail index $\gamma$ , and let

\begin{align*}\widehat\gamma_H=\frac{1}{k}\sum_{i=1}^k\log\left(\frac{Y_{n,n-i+1}}{Y_{n,n-k}}\right)\end{align*}

denote the Hill estimator. Then, the logarithmic-score estimator over the parameter space $\Gamma$ is the projection of the Hill estimator onto this interval, i.e.,

\begin{align*}\mathop{\mathrm{argmax}}\limits_{\gamma\in\Gamma}\frac{1}{k}\sum_{i=1}^k\log\left[f_\gamma^\circ\left(\frac{Y_{n,n-i+1}}{Y_{n,n-k}}\right)\right]=\Pi_{\Gamma}(\widehat\gamma_H),\end{align*}

where $\Pi_{\Gamma}(x)=\min\{\gamma_U,\max\{\gamma_L,x\}\}$ . In particular, if $\widehat\gamma_H\in\Gamma$ , then the logarithmic-score estimator coincides with the Hill estimator.

Remark 4.10. Segers (Reference Segers2001) shows that, under mild conditions on h, the estimator of the form

\begin{align*}\frac{1}{k}\sum_{i=1}^k h\!\left(\frac{Y_{n,n-i+1}}{Y_{n,n-k}}\right)\end{align*}

attains the smallest asymptotic variance when $h=\log$ . This suggests that it is inherently difficult to construct a scoring-rule-based estimator $\hat{\gamma}_k$ that outperforms the Hill estimator in terms of asymptotic efficiency.

The following corollary verifies the assumptions of Theorem 4.7 for LogS and $ES_\beta$ .

For the second statement, $\beta\lt1/\gamma_U$ guarantees the required bound simultaneously for all $\gamma\in\Gamma$ . For a fixed candidate $\gamma$ , however, the weaker condition $\beta\lt1/\gamma$ may suffice.

Proof. For LogS, the map $(\gamma,z)\mapsto \mathrm{LogS}(F_\gamma,z)$ is continuous on $\Gamma\times[1,\infty)$ , and LogS is strictly proper. Moreover, the calculations in the proof of Corollary 4.4 show that

\begin{align*}\sup_{\gamma\in\Gamma}|\mathrm{LogS}(F_\gamma,z)|\leq Az^{(1-\delta)/\gamma_G}.\end{align*}

Hence, the assumptions of Theorem 4.7 are satisfied for LogS.

For $ES_\beta$ , the map $(\gamma,z)\mapsto ES_\beta(F_\gamma,z)$ is continuous on $\Gamma\times [1,\infty)$ , and $ES_\beta$ is strictly proper. By the calculations in the proof of Corollary 4.4, for each fixed $\gamma$ with $\beta\lt1/\gamma$ , there exists a constant $C_\gamma\gt0$ such that

\begin{align*} |ES_\beta(F_\gamma,z)|\leq C_\gamma(1+z^\beta).\end{align*}

Since $\beta\lt1/\gamma_U$ , we have $\beta\lt1/\gamma$ for every $\gamma\in\Gamma$ . Moreover, it is seen from the calculations in the proof of Corollary 4.4 that the constant $C_\gamma$ can be chosen as a continuous function of $\gamma$ . Since $\Gamma$ is compact, $C(\gamma)$ is bounded on $\Gamma$ . Hence, there exists $C\gt0$ such that

\begin{align*} \sup_{\gamma\in\Gamma}|ES_\beta(F_\gamma,z)|\leq C(1+z^\beta).\end{align*}

\begin{align*}|ES_\beta(F_\gamma,z)|\leq C_\gamma\left(1+z^\beta\right).\end{align*}

Since $\gamma_G\in\Gamma$ , the condition $\beta\lt1/\gamma_U$ implies $\beta\lt1/\gamma_G$ . Hence, we may choose $\delta\in(0,1)$ such that

\begin{align*}\beta\lt\frac{1-\delta}{\gamma_G}.\end{align*}

It follows that for $z\geq 1$ ,

\begin{align*}C(1+z^\beta)\leq A z^{(1-\delta)/\gamma_G}\end{align*}

for some $A\gt0$ . Therefore, the uniform envelope condition in Theorem 4.7 is satisfied, and the result follows.

5.1. Ranking tail models: baseline setting

In the baseline setting, the purpose of the simulation study is to assess how well a LogS rule can identify the true tail index of a heavy-tailed distribution and how this depends on the threshold level and the sample size.

Samples of size $n \in \{10^3, 10^4, 10^5\}$ are generated from heavy-tailed distributions with true tail index $\gamma_G = 1$ . We consider two parametric families of distributions, namely the Fréchet and the Burr distributions. The Fréchet distribution is specified by the distribution function

\begin{align*}F(x) = \exp\!\left(-x^{-s}\right), \qquad x \gt 0,\end{align*}

where the shape parameter is set to $s = 1/\gamma_G$ . The Burr distribution is given by

\begin{align*}F(x) = 1 - \bigl(1 + x^{c}\bigr)^{-t}, \qquad x \gt 0,\end{align*}

with shape parameters chosen as $c = 1/\gamma_G$ and $t = 1$ . With these parameterizations, both distributions are characterized by the same tail index $\gamma_G$ .

Figure 1.

Empirical logarithmic scores in (3) (vertical axis) plotted against the number of upper order statistics k (horizontal axis) for candidate tail indices $\gamma \in \{0.8, 1, 1.2, 1.5\}$ . Left panels use Fréchet data-generating distributions, right panels use Burr distributions, and rows correspond to $n=10^3,10^4,10^5$ . Higher curves indicate better tail fit; the true value is $\gamma_G=1$ .

Figure 1. Long description

The image contains six line graphs. Panel A shows Fréchet DGP with n = 10^3, Panel B shows Burr DGP with n = 10^3, Panel C shows Fréchet DGP with n = 10^4, Panel D shows Burr DGP with n = 10^4, Panel E shows Fréchet DGP with n = 10^5, and Panel F shows Burr DGP with n = 10^5. Each graph plots empirical logarithmic scores (vertical axis) against the number of upper order statistics k (horizontal axis) for candidate tail indices γ = 0.8, γ = 1, γ = 1.2, and γ = 1.5. Higher curves indicate better tail fit, with the true value being γ = 1.

Within each sample, the threshold parameter k ranges from 50 to $n/4$ and is evaluated on an evenly spaced grid of 100 points. For each value of k, four logarithmic scores are computed by taking F to be a Pareto distribution with tail index

\begin{align*}\gamma \in \{0.8,\,1,\,1.2,\,1.5\}.\end{align*}

Using this setup, the empirical criterion is computed via (3). In Figure 1, the resulting scores of a single realization are plotted as functions of k, with one curve corresponding to each value of $\gamma$ . If the scoring rule is well aligned with the true tail behavior, the candidate value $\gamma = \gamma_G = 1$ should yield the highest log score across a broad range of k. To formalize this comparison, for each fixed k, the candidate values are ranked according to their empirical logarithmic scores. In particular, the highest-ranked value is

\begin{align*} \mathop{\mathrm{argmax}}\limits_{\gamma \in \{0.8,1,1.2,1.5\}} \mathrm{LogS}_k(\gamma),\end{align*}

with the remaining candidates ordered according to their corresponding scores. For the smaller sample sizes $n=10^3$ and $n=10^4$ , and in particular for the Burr model, the curve corresponding to $\gamma=1$ is not systematically above the competing curves; instead, the candidate values $\gamma=0.8$ and $\gamma=1.2$ often attain similar or larger scores over substantial ranges of k. By contrast, clearer signals appear in the remaining cases, especially for $n=10^5$ and for the Fréchet model, where the curve corresponding to $\gamma=1$ tends to dominate over a noticeable region of smaller k. This pattern is consistent with the expected bias-variance trade-off in tail inference: with larger n, a wider lower-k region supports stable discrimination of the correct tail index. This visual impression is further supported by a Monte Carlo experiment reported in Figure 2. For each configuration, 100 independent samples are generated, and the empirical logarithmic score is computed for the candidate values $\gamma \in \{0.8,1,1.2,1.5\}$ across the same range of thresholds. For each value of $k/n$ , we record whether the score is maximized at the true value $\gamma=1$ and plot the resulting proportion of such occurrences across the simulations. In both the Fréchet and Burr cases, when $n=10^3$ , there is no clear tendency for the proportion to approach one as $k/n$ goes to 0. As the sample size increases, however, the proportion rises in both models. For Fréchet data-generating distributions, this increase is particularly strong. When $n=10^5$ , the proportion is equal to one across the entire range of $k/n$ , indicating that the correct tail index is almost always selected. In the Burr case, the increase is observed only for smaller values of k/n, indicating that the correct tail index is consistently favored only further into the tail.

Figure 2.

Proportion of simulations (based on 100 Monte Carlo replications) in which the empirical logarithmic score in (3) is maximized at the true tail index $\gamma=1$ among candidate values $\gamma \in \{0.8, 1, 1.2, 1.5\}$ , plotted against the relative number of upper order statistics $k/n$ . The left plot uses Fréchet data-generating distributions and the right plot uses Burr distributions. Curves correspond to sample sizes $n=10^3,10^4,10^5$ . Higher values indicate that the logarithmic score more frequently identifies the true tail index.

It should be noted that both the Fréchet and Burr samples are evaluated against Pareto limit distributions. Although $\gamma=1$ gives the correct asymptotic tail index in both cases, finite-threshold exceedances are not exactly Pareto, so larger values of $k/n$ may still reflect pre-asymptotic distributional features. Hence, it is expected that the correct tail-index signal is clearest for smaller values of $k/n$ .

5.2. Ranking tail models under heterogeneous scaling

A practically relevant situation is one in which the observations are not identically distributed. For example, covariates may affect the scale of the underlying distribution. Although such settings go beyond the asymptotic theory developed in Section 4, we include the following simulation experiment to examine the robustness of the proposed approach under systematic scale variation. Specifically, we simulate independent baseline heavy-tailed variables Z from Burr and Fréchet distributions and multiply each realization by a scaling factor. That is,

\begin{align*}Y_i = X_i Z_i .\end{align*}

We follow the same setup as in Section 5.1 to generate the baseline heavy-tailed variable Z. We consider two different structures for the scaling factors X:

1. (i) $X_i^1={1+}\frac{i}{n}$ .

2. (ii) $X_i^2=1.5+0.5\sin\!\left(6\frac{i}{n}\pi\right)$ .

Both scaling structures imply that $X_i\in[1,2]$ . These two structures represent qualitatively different forms of scale heterogeneity and should be viewed as simple heuristic stress tests. They are chosen to separate two effects that may arise in practice: persistent changes in scale and recurrent scale fluctuations. The monotone structure captures the former, with the largest scale occurring only at the end of the sample. The periodic structure captures the latter, producing repeated high-scale regions so that several parts of the sample contain observations with the largest scale.

Importantly, in this simulation design the tail index is constant across all observations. The only source of heterogeneity arises from the scaling component. Consequently, the tail heaviness of the distribution remains unchanged, while only the scale varies across i. Under this setup, the model corresponding to $\gamma = 1$ correctly specifies the tail index of the data-generating process. Therefore, we would expect the logarithmic score (LogS) to attain its highest value at $\gamma = 1$ .

Figures 3 and 4 present the results of a single realization, while Figure 5 reports the proportion results. The figures are constructed in the same manner as those in the previous section and are therefore not described in further detail. The same k-grid and candidate set as in Section 5.1 are used, so changes in behavior can be attributed to the scaling mechanism. Even in the presence of heterogeneous scaling, the qualitative conclusions of Section 5.1 remain valid. In particular, the logarithmic score continues to favor the correct tail index for sufficiently large sample sizes and appropriate choices of k, and the most reliable inference is still concentrated in a stable lower-k region.

Figure 3.

Empirical logarithmic scores (vertical axis) versus k (horizontal axis) for Fréchet baseline samples with heterogeneous scaling, $Y_i=X_iZ_i$ , where $Z_i$ has tail index $\gamma_G=1$ . Candidate tail indices are $\gamma \in \{0.8, 1, 1.2, 1.5\}$ . Left panels use $X_i^1={1+}i/n$ , right panels use $X_i^2=1.5+0.5\sin(6\pi i/n)$ , and rows correspond to $n=10^3,10^4,10^5$ .

Figure 3. Long description

The image contains six line graphs. Each graph shows empirical logarithmic scores on the vertical axis and k on the horizontal axis. The graphs compare different tail indices (gamma values: 0.8, 1, 1.2, 1.5) for Fréchet baseline samples with heterogeneous scaling. Panel A shows linear scaling with n = 10^3. Panel B shows sinusoidal scaling with n = 10^3. Panel C shows linear scaling with n = 10^4. Panel D shows sinusoidal scaling with n = 10^4. Panel E shows linear scaling with n = 10^5. Panel F shows sinusoidal scaling with n = 10^5. Each graph displays multiple lines representing different gamma values, illustrating how the empirical logarithmic scores vary with k for each scaling method and sample size.

Figure 4.

Empirical logarithmic scores (vertical axis) versus k (horizontal axis) for Burr baseline samples with heterogeneous scaling, $Y_i=X_iZ_i$ , where $Z_i$ has tail index $\gamma_G=1$ . Candidate tail indices are $\gamma \in \{0.8, 1, 1.2, 1.5\}$ . Left panels use $X_i^1={1+}i/n$ , right panels use $X_i^2=1.5+0.5\sin(6\pi i/n)$ , and rows correspond to $n=10^3,10^4,10^5$ .

Figure 5.

Proportion of simulations (based on 100 Monte Carlo replications) in which the empirical logarithmic score in (3) is maximized at the true tail index $\gamma=1$ among candidate values $\gamma \in \{0.8, 1, 1.2, 1.5\}$ , plotted against the relative number of upper order statistics $k/n$ . The top panels correspond to Fréchet data-generating distributions and the bottom panels to Burr distributions. The left panels use a linearly varying scale, while the right panels use a sinusoidally varying scale. Curves correspond to sample sizes $n=10^3,10^4,10^5$ . Higher values indicate that the logarithmic score more frequently identifies the true tail index.

Figure 5. Long description

Four line graphs compare the proportion of simulations where the empirical logarithmic score is maximized at the true tail index among candidate values against the relative number of upper order statistics. Panel A: Fréchet DGP, Linear scaling Xi = Xi^1. The y-axis represents the proportion, ranging from 0 to 1, and the x-axis represents the relative number of upper order statistics, ranging from 0.05 to 0.25. Panel B: Fréchet DGP, Sinusoidal scaling Xi = Xi^2. The y-axis represents the proportion, ranging from 0 to 1, and the x-axis represents the relative number of upper order statistics, ranging from 0.05 to 0.25. Panel C: Burr DGP, Linear scaling Xi = Xi^1. The y-axis represents the proportion, ranging from 0 to 1, and the x-axis represents the relative number of upper order statistics, ranging from 0.05 to 0.25. Panel D: Burr DGP, Sinusoidal scaling Xi = Xi^2. The y-axis represents the proportion, ranging from 0 to 1, and the x-axis represents the relative number of upper order statistics, ranging from 0.05 to 0.25. Each panel includes three lines representing different sample sizes: n = 1000, n = 10000, and n = 1e+05. Higher values indicate that the logarithmic score more frequently identifies the true tail index.

5.3. Finite-sample performance of score-optimization estimators

We now examine the finite-sample behavior of tail-index estimators based on the Energy score, $\hat{\gamma}_k(ES_\beta)$ , to illustrate the finite-sample behavior of Energy-score-based estimation. Their performance is compared with that of the classical Hill estimator over a range of tail indices. Recall that the parameter $\beta$ must satisfy $\beta \lt 1/\gamma^U$ .

Samples of size $n = 10^4$ are generated from Fréchet distributions. Companion Burr runs were also conducted but are not shown, as they display the same qualitative sensitivity in k and $\beta$ and were omitted for brevity. As in the previous simulation studies, the remaining parameters are chosen analogously. The number of upper order statistics k is varied over 25 equally spaced values between 50 and $n/4$ . The Energy score was evaluated over the compact parameter space

\begin{align*}\Gamma=[\gamma_L,\gamma_U]=[\,0.8\gamma_G,\,2\gamma_G\,],\end{align*}

using a grid of 150 equidistant candidate values. The interval is a heuristic choice, included to contain the true value while also allowing the estimator to select substantially heavier-tailed alternatives. To study the sensitivity of the estimator to the parameter $\beta$ in the score, three different values were considered,

\begin{align*}\beta_1 = \frac{1}{{\gamma_U}} - 0.001, \qquad\beta_2 = 0.8\beta_1, \qquad\beta_3 = 0.7\beta_1.\end{align*}

For each tail index $\gamma_G$ , the estimated values $\hat{\gamma}_k$ obtained from the three Energy-score-based estimators and the Hill estimator are plotted against k. The true value $\gamma_G$ is shown as a horizontal reference line.

Figure 6 reports the Fréchet results for a single realization. Across all panels, the choice of $\beta$ has only a limited impact on the resulting Energy-score estimator. Moreover, the three Energy-score estimators are almost indistinguishable from the Hill estimator over most of the considered range of k. Thus, in this setting, optimizing the Energy score leads to essentially the same finite-sample behavior as the classical Hill estimator. This should be viewed in light of Remarks 4.9 and 4.10. These results give the logarithmic-score estimator, and hence the Hill estimator, a natural benchmark role. Hence, one should not expect the Energy-score estimators to provide a clear efficiency improvement in this experiment. At the same time, the near-overlap of the curves indicates that using the Energy-score does not lead to a noticeable deterioration of the estimates. Table 1 supports the same conclusion. The reported bias and variance are estimated from 100 independent simulation runs using the usual sample mean and sample variance of the estimates. For all four values of $\gamma_G$ and all three choices of k, the Energy-score estimators and the Hill estimator have very similar bias and variance. The differences between the three choices of $\beta$ are also small, indicating that the finite-sample performance is not particularly sensitive to the precise value of $\beta$ in this simulation design.

Figure 6.

Tail-index estimates $\hat{\gamma}_k$ (vertical axis) as functions of the number of upper order statistics k (horizontal axis) for Fréchet samples with $n=10^4$ . Panels vary the true tail index $\gamma_G$ . Each panel compares the Hill estimator with three Energy-score estimators based on $\beta_1,\beta_2,\beta_3$ ; the horizontal reference line marks $\gamma_G$ .

Table 1.

Bias and variance of the Hill estimator and three Energy-score estimators corresponding to $\beta_1,\beta_2,$ and $\beta_3$ , computed from 100 simulated samples drawn from a Fréchet distribution with sample size $n = 10^4$ . Results are reported for different values of the extreme value index $\gamma$ and for three choices of the number of upper order statistics k.

Table 1. Long description

A table with bias and variance of the Hill estimator and three Energy-score estimators. The table has 12 rows and 8 columns. Columns are labeled with k values and corresponding variance and bias. Rows are labeled with Hill, beta 1, beta 2, and beta 3 for different gamma values. Each row contains values for variance and bias under different conditions. The table shows detailed numerical comparisons of bias and variance for different estimators and conditions.

Figure 7.

USAutoBI claim severity analysis ( $n=1{,}340$ ). Panel (a) compares theoretical Pareto upper-tail quantiles (horizontal axis) with empirical quantiles (vertical axis). Panels (b)–(c) plot logarithmic score (vertical axis) against the number of upper order statistics k (horizontal axis) for five Pareto predictive distributions with $\gamma\in\{0.3,0.5,0.8,1,1.3\}$ . Panel (c) is a zoom of panel (b) over the lower 25% of the k-range. Panels (d)–(f) show score differences relative to the $\gamma=1$ model with pointwise 95% confidence intervals.

Figure 7. Long description

Panel A: A scatter plot compares theoretical Pareto upper-tail quantiles on the horizontal axis with empirical quantiles on the vertical axis. Panel B: A line graph plots logarithmic score on the vertical axis against the number of upper order statistics k on the horizontal axis for five Pareto predictive distributions with different gamma values. Panel C: A line graph similar to Panel B but zoomed in on the lower 25% of the k-range. Panel D: A line graph shows score differences relative to the model with pointwise 95% confidence intervals for gamma equals 0.5. Panel E: A line graph shows score differences relative to the model with pointwise 95% confidence intervals for gamma equals 0.8. Panel F: A line graph shows score differences relative to the model with pointwise 95% confidence intervals for gamma equals 1.3.