1 Introduction↩︎

In many scientific fields one has to predict rare events of the type \(\{Y>y_0\}\), where an unobserved response random variable \(Y\) exceeds a large threshold \(y_0\). Assuming that one observes a vector of covariates \(X = (X_i)_{i=1}^d\), such binary predictors can always be written in the form \(\{h(X)> h_0\}\), for some measurable function \(h\). The optimal predictors, in a certain natural sense, can be explicitly characterized via density ratios, as in Neyman-Pearson’s Lemma (see Section 2.1). Nevertheless, this is the beginning rather than the end of the story. Except in a handful of cases where closed-form solutions exist (Examples 1–3), the optimal predictors need to be estimated. While the estimation of density ratios has received some attention [1] it remains a methodologically and practically challenging problem [2], [3]. In this paper, we are interested in the regime of extremes, where the events \(\{Y>y_0\}\) are rare or extremely rare. An equivalent perspective to this problem is from the point of view of classification, where there is an extreme imbalance between the two classes of labels \(\{Y\le y_0\}\) and \(\{Y>y_0\}\) in the training data. Imbalanced classification problems arise naturally and are ubiquitous in science, medicine, engineering, etc [4]–[6]. Imbalanced classification has received considerable theoretical attention [3], [7], [8]. Nevertheless, the theoretical study of the optimal extreme event prediction or, equivalently, binary classification in the extremely imbalanced regime has received relatively less attention. To the best of our knowledge, [9] is the first paper that has addressed this problem by arguing that the traditional empirical risk minimization approaches do not enjoy provable guarantees in the regime of extremes. They proposed a principled solution based on extreme value theory by considering risk minimization after conditioning on asymptotically extreme events, which has led to an active area of research [10]–[13]. In this paper, we offer an alternative take on the problem from the perspective of optimal extreme event prediction rather than imbalanced classification.

An important starting point of this investigation is the framing of the notion of optimality. Colloquially, prediction of the unobservable event \(\{Y>y_0\}\) via the observed covariates amounts to raising an alarm when an event \(\{h(X)> h_0\}\) occurs. Let the thresholds \(y_0:=F_Y^\leftarrow(p)\) and \(h_0:=F_{h(X)}^\leftarrow(q)\) be generalized quantiles of the distributions of \(Y\) and \(h(X)\), respectively, for some \(p,q\in (0,1)\). Then, under mild regularity conditions, the alarm rate equals \(1-q = \mathbb{P}[h(X)> F_{h(X)}^\leftarrow(q)]\) and the event rate is \(1-p = \mathbb{P}[Y>F_Y^\leftarrow(p)]\). Given certain alarm and event rates, the optimal predictors are the ones that maximize the precision: \[\mathbb{P}[ Y> F_Y^\leftarrow(p)\, |\, h(X)> F_{h(X)}^\leftarrow(p)].\] It is often natural to seek calibrated predictors where \(p=q\), i.e., the alarm and event rates coincide. Then, in the regime of extremes, as \(p\uparrow 1\), the asymptotic precision is precisely the tail-dependence coefficient \[\label{e:lambda-intro} \lambda(Y,h(X)) = \lim_{p\uparrow 1} \mathbb{P}[ Y> F_Y^\leftarrow(p)\, |\, h(X)> F_{h(X)}^\leftarrow(p)],\tag{1}\] which is a well-known quantity in extreme value theory [14].

In the general setting when \(Y\) and \(X\) are jointly regularly varying, and the predictors \(h\) are continuous and positive homogeneous functions, the tail dependence coefficient \(\lambda(Y,h(X))\) is well-defined, and it can be expressed as an integral functional with respect to the angular measure (cf Propositions 2 and 3). Specifically, \[\label{e:lambda-expression-intro} \lambda(Y,h(X)) = \frac{1}{\mathbb{E}[U]} \mathbb{E}[ U \wedge (1-U) h(\Theta) ],\tag{2}\] where the random variable \(U\in [0,1]\) and the random vector \(\Theta = (\Theta_i)_{i=1}^d\) characterize the joint extremal behavior of \(Y\) and \(X\) as follows: \[\frac{(Y,X)}{\tau(Y,X)} \, | \, \{\tau(Y,X)>t\} \stackrel{d}{\to} (U,(1-U)\Theta),\;\; as t\to\infty.\] Here \(\tau(y,x):= y_+ + \tau_X(x)\) is a positive \(1\)-homogeneous gauge function codifying the notion of multivariate regular variation (Definition 4) and \(\Theta\) is a random vector taking values in the generalized unit sphere \(S_{\tau_X} =\{x \in \mathbb{R}^d\, :\, \tau_X(x) = 1\}\).

In view of 1 , we seek homogeneous predictors that maximize the tail-dependence coefficient \(\lambda(Y,h(X))\) and refer to them as asymptotically optimal. Given the joint distribution of \((U,\Theta)\), Relation 2 shows that finding asymptotically optimal homogeneous predictors amounts to solving a variational problem. Namely, finding a homogeneous function \(h\) in an infinite-dimensional class of functions, that maximizes the integral functional in 2 , under the asymptotic calibration constraints: \[\mathbb{P}[ Y>t ] \sim \mathbb{P}[h(X)>t],\;\; as t\to\infty.\] Here and throughout the paper by \(a_t\sim b_t\) we mean asymptotic equivalence in the sense that \(a_t/b_t \to 1\) as \(t\to\infty\).

The first main contribution of the paper is a solution of the aforementioned variational problem (Theorem 2). Then, in Section 3.3, we demonstrate how this solution yields asymptotically optimal homogeneous predictors in the context of the general class of Breiman type models. It turns out that homogeneous optimal predictors can be expressed in terms of the homogenous function \[g^{\rm (opt)}(x):= \tau_X(x) \cdot \frac{q_\alpha(x/\tau_X(x))}{ 1- q_\alpha(x/\tau_X(x))},\] where \(\theta \mapsto q_\alpha(\theta)\) is the conditional \(\alpha\)-level quantile of the tilted conditional probability distribution \(b(\theta)^{-1} (1-u) p_{U|\Theta}(du|\theta)\), where \(p_{U|\Theta}(du |\theta)\) is the regular conditional distribution of \(U\) given \(\Theta\), and \(b(\theta)\) is a normalization constant.

All homogeneous functions of the above form are solutions to the optimal variational problems with constraints governed by the parameter \(\alpha\in (0,1)\) (Theorem 2). Thus, to find an optimal and calibrated predictor, in practice, one merely needs to estimate the quantile functions \(q_\alpha\) and find \(\alpha\) so that \(\mathbb{P}[g^{\rm (opt)}(X)>t] \sim \mathbb{P}[Y>t],\;t\to\infty\). Inference theory for the optimal homogeneous predictors in the peaks-over-threshold framework developed in Section 4. In this setting, by using a contiguity argument, we established the universal consistency of these optimal homogeneous predictors in Theorem 3 under mild regularity conditions. This result may be viewed as a counterpart of the Stone universal consistency theory in the context of extremes [15].

In Section 5, we discuss the implementation of the estimators and illustrate their performance in comparison with optimal oracle predictors, both in the context of discrete and continuous angular measures. A general-purpose, practical, and nearly tuning free implementation of the optimal homogeneous predictor methodology was built by using powerful quantile regression forests [16], [17], which were used to estimate the quantile functions \(q_\alpha\). The resulting predictors nearly match the performance of the oracle predictors in a variety of simulated models (Figures 1 and 2). Furthermore, their application to the open challenge of predicting extreme X-class solar flare events achieves state-of-the-art results without sophisticated feature engineering and tuning (see Section 11). R code and R Shiny Apps illustrating the methodology are freely available on [18], [19].

The paper contains a variety of additional contributions of independent interest. First, in Proposition 1 we provide a characterization of the optimal (non-asymptotic) predictors for Archimedean copula models; A general contiguity argument justifying the peaks-over-threshold based inference was developed (Section 4.2); A flexible family of Pareto-Breiman models was studied where regular variation holds in a stronger, total-variation sense (Proposition 5). The case of Pareto-Dirichlet models has a particularly elegant solution (Proposition 20), while a large class of spectrally discrete models offer insights to further phenomena and applications (Section 5.2). A curious phenomenon of perfect asymptotic precision arises for certain spectrally discrete models, when all underlying latent factors influencing \(Y\) are involved in forming the covariates \(X\) (see Section 9).
The rest of the paper is structured as follows. In Section 2, we begin with the formulation of the optimal prediction framework and present a general characterization results of optimal predictors in terms of density ratios (Theorem 1). This naturally leads to several examples where closed-form optimal (non-asymptotic) predictors can be derived (i.e., Sections 2.1 and 2.2). In Section 2.3, we conclude the preliminaries by relating the notion of asymptotic optimality and tail dependence.

In Section 3, we present the formulation and solution to the variational problem on optimal homogeneous prediction, where the main result is Theorem 2. Section 4 presents an inference theory for these optimal predictors in the context of peaks-over-threshold, where we carefully treat the distribution shift problem via contiguity argument. The main result is Theorem 3, which shows that the so-constructed estimators achieve the asymptotically optimal precision. This can be seem as a counterpart to the celebrated universal consistency results of Stone. Section 5 briefly discusses methodological aspects and illustrates the non-asymptotic, practical performance of the optimal homogeneous predictors with simulations. Further technical details and some lengthy proofs are given in the Supplement.

2 Preliminaries↩︎

In this section, we introduce the notion of optimal event prediction considered in this paper. A general characterization result and several examples of closed-form optimal predictors are discussed, which may be of independent interest and can serve as benchmarks. The extremal counterpart to optimal event prediction is then formulated and related to asymptotic tail dependence. This sets the stage for the problem formulation discussed in Section [sec:problem].

2.1 Optimal event prediction: A Neyman-Pearson perspective↩︎

Consider a vector of covariates \(X= (X_i)_{i=1}^d\) taking values in \(\mathbb{R}^d\) and an unobserved positive response random variable \(Y\). We would like to optimally predict the event \(\{Y > y_0\}\) in terms of \(X\). All such predictors are of the form \(\{h(X)>h_0\}\), for some Borel measurable function \(h\) and a suitable calibration threshold \(h_0\in \mathbb{R}\).

The following result is a restatement of Theorem 2.1 in [20]. It provides a general characterization of the optimal calibrated predictors.

Theorem 1 leads to closed form optimal predictors in a handful but important cases.

The above example leads to explicit characterizations of the optimal event-predictors as well as their precisions in the general class of linear models [20]. In particular, since Gaussian vectors can always be cast in the form of ?? , where \(g\) is a linear function and \(\sigma\) is non-random, Example 1 immediately implies the following.

The above example is rather natural since max-stable laws are positively dependent and hence the extremes of \(Y\) go with the extremes of \(X\). See also Example 6, below.

2.2 Further examples: Archimedean copula↩︎

The closed-form expressions of optimal predictors as in the previous examples are rare. In this section, we present such an expression for the broad class of Archimedean copula with completely monotone generators. It can be of independent interest and shall serve as a benchmark for our optimal prediction inference methodology.

Suppose that the random vector \((Y,X)=(Y,X_1,\cdots,X_d)\) has continuous marginal cdf’s \(F_Y(t) = \mathbb{P}[Y\le t]\) and \(F_{X_i}(t):= \mathbb{P}[X_i\le t],\;t\in \mathbb{R},\;i=1,\cdots,d\) and the copula \[C(y,x):= \mathbb{P}[F_Y(Y)\le y, F_{X_i}(X_i)\le x_i,\;i=1,\cdots,d],\;\;y\in (0,1), x=(x_i)_{i=1}^d\in (0,1)^d,\] The copula \(C\) is said to be Archimedean if \[\label{e:archimedean} C(y,x) = \psi\Big( \psi^{-1}(y)+ \psi^{-1}(x_1)+\cdots+\psi^{-1}(x_d) \Big),\tag{3}\] where \(\psi:(0,\infty) \to (0,1)\) is a decreasing function referred to as the generator of the copula.

It is known that \(C\) is a valid copula provided \(\psi\) is a \((d+1)\)-monotone function [21]. Here, we focus on the important particular case of Archimedean copula, where the generator \(\psi\) is completely monotone and hence \((d+1)\)-monotone for all \(d\). Such functions have the representation \[\label{e:Bernstein} \psi(u) = \int_0^\infty e^{- ux } \mu(dx),\;\;u>0\tag{4}\] for some unique Borel measure \(\mu\) on \([0,\infty)\) [22]. Since we deal with copulas, the measure \(\mu\) will not charge \(\{0\}\).

The proof is given in Section 7.1.2, below. This result furnishes a complete characterization of optimal event prediction in an interesting family of copula models. We will illustrate it further with an important Archimedean copula arising in extreme value theory.

2.3 Tail dependence and optimal extreme event prediction↩︎

In this paper, we are primarily interested in predicting extreme events of the form \[\{Y> F_Y^\leftarrow(p)\},\] as the level \(p\) becomes extreme, i.e., as \(p\uparrow 1\). For simplicity, assume that \(p\in F_Y(\mathbb{R}) = \{F_Y(y),\;y\in \mathbb{R}\}\) and focus on the natural case of balanced predictors, where one has \[\label{e:balanced-predictors} 1-p = \mathbb{P}[ Y> F_Y^\leftarrow(p) ] = \mathbb{P}[ h(X) > F_{h(X)}^\leftarrow(p) ].\tag{5}\]

In this case, the conditional probability of a true positive alarm, referred to as the precision of the predictor, is: \[\label{e:lambda-p} \lambda_p(Y,h(X)) := \mathbb{P}[Y>F_Y^\leftarrow(p)\, |\, h(X)>F_{h(X)}^\leftarrow (p)].\tag{6}\] Observe that, as the level \(p\) becomes extreme, i.e., \(p\uparrow 1\), the above precision recovers the well-known notion of tail-dependence coefficient: \[\label{e:lambda-tdep-def} \lambda(Y, h(X)) := \lim_{p\uparrow 1} \mathbb{P}[Y>F_Y^\leftarrow(p)\, |\, h(X)>F_{h(X)}^\leftarrow (p)],\tag{7}\] whenever the limit exists [14], [23]. The random variables \(Y\) and \(h(X)\) will be referred to as tail-independent if \(\lambda(Y, h(X))=0\); otherwise we will sometimes call them tail-dependent.

Relation 7 motivates the following notion of extremal optimality [20].

Relation 8 indicates that \(\lambda^{{\rm (opt)}}(Y,X)\) is the best asymptotic precision that can be achieved among all balanced calibrated predictors of the extreme events \(\{Y> F_Y^\leftarrow(p)\}\). Theorem 1 implies that the density ratios in ?? yield fixed-\(p\) as well as asymptotically optimal predictors. This is the gold-standard in extreme event prediction, which can in fact be achieved in a number of interesting models (see the Examples 1-4 above and their follow-ups in 5–7 below).

In general, however, the density ratios are hard to compute and estimate [1]. Thus, it is challenging to characterize the optimal extremal predictors and their extremal precisions in full generality. Therefore, it is of interest to limit the notion of extremal optimality to a class of predictors.

We conclude this section by following-up on Examples 1 – 4, where the optimal extremal predictors and their precisions can be derived. The next section focuses on \({\cal G}\)-optimal prediction for the class of homogeneous functions under the framework of joint regular variation for \(Y\) and \(X\).

While elegant, the above examples are the only a handful of cases we know of, where the optimal (extremal) predictors and their (asymptotic) precisions can be obtained in closed form. This motivates the main pursuit of the present paper on characterizing optimal extremal predictors over a smaller but still very rich class of homogeneous predictors.

3 A characterization of the optimal homogeneous predictors↩︎

In this section, we present the core contribution of this paper. We start with the main framework that relies on an assumption of joint regular variation between the response and the predictors. In this context, we cast the main problem of finding asymptotically optimal homogeneous predictors as a calculus of variations problem relative to the angular measure. The general solution of this problem is presented in Theorem 2. The section concludes with a rejoinder, where the solutions to the general variational problem are related to optimal homogeneous predictors in the context of the generalized Breiman models.

3.1 Joint regular variation and problem formulation↩︎

Let \(Z\) be a random vector taking values in \(\mathbb{R}^k\) and let \(\tau :\mathbb{R}^k\to \mathbb{R}_+ :=[0,\infty)\) be a non-negative continuous \(1\)-homogeneous function, i.e., \(\tau(c\cdot x) = c\tau(x),\;\forall c\ge 0,\;x\in \mathbb{R}^k\).

For more details and an equivalent definition, see Section 8.1 below. In particular, the first convergence in Relation ?? holds if and only if \(x\mapsto \mathbb{P}[\tau(Z) > x]\) is a regularly varying function at infinity with exponent \(-\alpha\). Therefore, ?? is equivalent to \[\label{e:d:tau-rv-b40t41} b(t)\mathbb{P}[ \tau(Z) > t ] \longrightarrow c_Z t^{-\alpha},\;(n\to\infty) \; and \; \frac{Z}{\tau(Z)} \, \vert\, \tau(Z)>r \stackrel{d}{\longrightarrow} \Theta_Z,\;\;(r\to\infty)\tag{9}\] where \(b(t)\sim L(t)t^\alpha\), as \(t\to\infty\), for some slowly varying function \(L(\cdot)\) such that \(b(a_n)\sim n\), as \(n\to\infty\). Thus, for an for an \(X\) as in Definition 4, the function \(\{b(t)\}\) in 9 is unique modulo asymptotic equivalence [25]. When we need 9 , we shall equivalently write \(X\in {\rm RV}_\alpha(\mathbb{R}^k,\{a_n\},b(\cdot),\;c_Z,\tau,\sigma)\).

The following is a key result that will allow us to formulate the optimal extreme event prediction problem in terms of a calculus of variations problem for the angular measure. See also Proposition 2.5 in [23] and Theorem 2.1 in [30], for similar results.

For the sake of completeness, the proof of Proposition 2 is given in Section 8.3.1.

In the rest of the section, we will specialize the above regular variation concept to the setting considered in this paper by considering a specific \(\tau\). Namely, let \[Z = (Y,X),\] where \(Y\) is a positive response variable and \(X = (X_i)_{i=1}^d\) is a \(d\)-dimensional vector of covariates. We will assume that \(Y\) and \(X\) are jointly regularly varying, which will be made precise next. For simplicity and without loss of generality, we will assume the exponent of regular variation equals one, i.e.: \[\alpha=1.\]

Suppose now that \(X\) and \(Y\) are jointly regularly varying in the sense of Assumption 1 and let \[\label{e:tilde-U-and-tilde-Theta} \widetilde{U} := \frac{Y}{\tau(Y,X)},\;\; and \;\;\widetilde{\Theta}:= \frac{X}{\tau_X(X)}.\tag{10}\] Observe that with the chosen polar coordinate parameterization \(\tau(Y,X) = Y+\tau_X(X)\), we have: \[\label{e:Z47tau40Z41-via-tilde-U-and-tilde-Theta} \frac{(Y,X)}{\tau(Y,X)} = \Big( \widetilde{U}, (1-\widetilde{U}) \widetilde{\Theta}\Big).\tag{11}\] Hence, in view of Definition 4, we obtain that for all \(t>0\), as \(n\to\infty\) and \(r\to\infty\): \(n\mathbb{P}[ \tau(Y,X) > ta_n] \longrightarrow t^{-1},\) and \[\label{e:U-THETA} (\widetilde{U}, (1-\widetilde{U}) \widetilde{\Theta})\, | \, \tau(Y,X)>r \stackrel{d}\longrightarrow\Theta_{Y,X}=:(U,(1-U)\Theta).\tag{12}\]

In the sequel, we shall use the \((U,(1-U)\Theta)\) parameterization of the angular vector \(\Theta_{Y,X}\) to formulate the optimal homogeneous prediction problem. Observe that if \(U<1\), then \(\Theta=W/(1-U)\) can be identified from \((U,W) := \Theta_{Y,X}.\) Whenever \(U=1\), however, \(\Theta\) is undefined, though unimportant since \((1-U)\Theta = 0\). In this case, by convention, we shall define \(\Theta := \theta_*\), for some arbitrary fixed \(\theta_*\in S_{\tau_X}\) such that \(\sigma_X(\{\theta_*\}) = 0\), where \(\sigma_X\) is the angular measure of \(X\) (cf Remark 2).

Note that (ii) is a technical condition, which is always satisfied if \({S_X}\) is compact.

The proof of this lemma is given in Section 8.3.1. The next proposition leads us to the formulation of the optimal homogeneous prediction problem as a calculus of variations problem.

The claim follows from Proposition 2 and Lemma 1, but its proof is given Section 8.3.1, for the sake of completeness.
Observe that every non-negative homogeneous function \(h:\mathbb{R}^d\to\mathbb{R}_+\), support-dominated by \(\tau_X\), is uniquely determined by its values on the unit sphere \(S_{\tau_X}\), since for all \(\tau_X(x)>0\), we have \(h(x) = \tau_X(x) \cdot h(\theta)\), where \(\theta:= x/\tau_X(x)\). Thus, Proposition 3 above entails that the optimal extremal predictors via homogeneous functions can be obtained by the functions \(g:{S_X}\to \mathbb{R}_+\) that minimize ?? under the constraints ?? . This leads us to the following optimization problem.

Using the identity \(|a-b| = a+b -2(a\wedge b)\), valid for all \(a,b\in \mathbb{R}\), Problem 1 can be equivalently formulated as a minimization problem:

In the following sub-section we provide a general solution to Problem 2 (equivalently, 1) for all possible spectral measures arising under Assumption 1.

3.2 The general solution to the variational problems↩︎

In this section, we will focus on the abstract optimization Problem 2 in the general case of non-negative measurable functions \(g\), where the pair \((U,\Theta)\) takes values in \([0,1]\times {S_X}\) and has probability distribution \(p_{U,\Theta}\). To exclude trivialities, we will assume that \[\label{e:condition-on-U} \mathbb{P}[U=0]<1\;\; and \;\;\mathbb{P}[U=1]<1,\tag{14}\] which corresponds precisely to ?? in the joint regular variation condition on \(Y\) and \(X\). For the purposes of this section, however, the variables \(Y\) and \(X\), will play no role.

The joint distribution \(p_{U,\Theta}(dud\theta)\) can always be written as \[\label{e:p95U44Theta} p_{U,\Theta}(dud\theta) = p(du|\theta) p_\Theta(d\theta),\tag{15}\] where \(p(du|\theta)\) is a regular conditional probability i.e., \(p(du|\theta)\) is a probability transition kernel [31].
We proceed with some convenient notation. Let \[\label{e:b40theta41} b(\theta):= \int_{[0,1]} (1-u)p(du |\theta)\tag{16}\] and observe that \(p_\Theta( \{0<b(\theta)<1\} ) = 1\), by 14 .

The functional \(I(\cdot)\) in 17 can be written as: \[\label{e:I40g41} I(g) := \int_{{S_X}} \Big\{ \int_{[0,1)} | u - (1-u)g(\theta)| p(du|\theta) \Big\} p_{\Theta}(d\theta).\tag{17}\] In view of 16 , by an application of the triangle inequality and Fubini’s theorem, we see that \(I(h)<\infty\), for all measurable \(h\), such that \[\label{e:C40g41} C(h) := \mathbb{E}[ (1-U) |h(\Theta)|] = \int_{{S_X}} |h(\theta)| b(\theta) p_\Theta(d\theta)<\infty.\tag{18}\] We will denote this class of functions by \(L^1(b\cdot p_\Theta)\) and the non-negative such functions by \(L_+^1(b\cdot p_\Theta)\).

The functional \(C(\cdot)\) in 18 will be referred to as the constraint functional. This terminology is motivated by Problem 2, which amounts to minimizing \(I(g)\) over all \(g\in L_+^1(b\cdot p_\Theta)\) such that \(C(g) = \mathbb{E}[U] = \sigma_Y\). Note also that the constraint functional \(C(g)\) is precisely the weighted \(L^1\)-norm of \(g\) in \(L^1(b\cdot p_\Theta)\), so that for \(C(g)<\infty\) is equivalent to \(g\in L^1(b\cdot p_\Theta)\).
Consider now the family of cumulative probability distribution functions indexed by \(\theta\in {S_X}\): \[\label{e:F-theta} F_\theta(t ):= \frac{1}{b(\theta)} \int_{ \{0\le u\le t\}} (1-u)p(du |\theta),\;\; t\in [0,1].\tag{19}\] Observe that the \(F_\theta\)’s are the distribution functions of probability measures concentrated on \([0,1)\). Notably, these measures, also denoted by \(F_\theta\), have no atoms at \(1\), i.e., \(F_\theta(1) = F_\theta(t-) = 1\) and \[\label{e:no-atom-at-1} F_\theta(\{1\}) = F_\theta(1) -F_\theta(1-) =0,\;\; for all \theta\in {S_X}.\tag{20}\] The following lemma establishes a precise formula for the Gâteaux directional derivative of the functional \(I\), which will be key. Its proof is given in Section 8.3.2.

We proceed with some further notation and auxiliary results needed to formulate and prove the main result. It turns out one can obtain an essentially explicit solution to Problem 2 in terms of the quantiles of the weighted conditional probability distribution \(F_\theta\) in 19 . Let \[\label{e:q-alpha} q_\alpha(\theta):= \inf\Big\{ t \ge 0\, :\, F_\theta(t ) \ge \alpha \Big\},\;\;\alpha\in [0,1]\tag{21}\] be the left-continuous generalized inverse of \(F_\theta(\cdot)\), with the small modification that above infimum is restricted to \(t\in [0,\infty)\) since the distributions \(F_\theta\) are supported on \([0,1)\). This implies that \(q_0(\theta) = 0\). We shall also need the right-continuous generalized inverse \[\label{e:q-alpha-plus} q_{\alpha + }(\theta) = \inf\Big\{ t \ge 0\, :\, F_\theta(t ) > \alpha \Big\},\tag{22}\] where by convention we set \(q_{1+}(\theta) = q_1(\theta) =\sup\{ t \in \mathbb{R}\, :\, F_\theta(t) <1\}\le 1\). Note that \[q_{\alpha+}(\theta) = \lim_{\beta>\alpha,\;\beta\downarrow \alpha} q_{\beta}(\theta),\;\; for all 0\le \alpha <1.\] The fact that the distributions \(F_\theta\)’s have no atoms at \(1\) (cf 20 ), implies the important observation that \[\label{e:q-alpha-43601} q_\alpha(\theta) \le q_{\alpha+}(\theta) < 1,\;\; for all 0\le \alpha<1 and \theta\in {S_X}.\tag{23}\]

The following general result will be useful. Its proof is given in Section 8.3.2.

The next lemma provides a technical ingredient needed for our first result.

The following proposition characterizes a broad class of solutions to Problem 2 under arbitrary positive values of the constraint functional. It shows that the solutions can be obtained from functions sandwiched between the left- and right-continuous quantile functions of the \(F_\theta\) distributions.

Proposition 4 provides a rich family of solutions to Problem 2, where the constraint \(\sigma_Y\) is replaced by a potentially arbitrary finite value \(C(g)\). Thus, to complete the solution to Problem 2, it remains to demonstrate that it is always possible to choose such a function \(g\) that meets the constraint \(\sigma_Y = \mathbb{E}[U] = C(g).\) To this end, the conditional quantile functions \(q_\alpha(\theta)\) and \(q_{\alpha+}(\cdot)\) will play an important role.

For all \(0<\alpha<1\) and \(\lambda \in [0,1]\), introduce the functions \[\label{e:g-alpha-lambda} g_\alpha^{(\lambda)} (\theta):= \frac{q_\alpha^{(\lambda)}(\theta)}{1- q_\alpha^{(\lambda)}(\theta)},\tag{31}\] where \[\label{e:q-alpha-lambda} q_\alpha^{(\lambda)}(\theta):= (1-\lambda) q_\alpha(\theta) + \lambda q_{\alpha+}(\theta).\tag{32}\] Recall that by 23 , we have \(q_{\alpha+}^{(\lambda)}(\theta)<1\) and hence \(g_\alpha^{(\lambda)}(\theta) <\infty\), for all \(0\le \alpha <1\) and \(\lambda \in [0,1]\). Note, however, that it is often possible to have \(q_1(\theta) = q_{1+}(\theta) = 1\) in which case \(g_1(\theta)=\infty\), by convention.

The following lemma, proved in Section 8.3.2 below, is essential. It establishes that for all \(\alpha \in [0,1)\), the predictors \(g_\alpha^{(\lambda)}(\cdot) ,\;0\le \lambda\le 1\) belong to \(L_+^1(b\cdot p_\Theta)\) and shows other useful properties of the constraint functional.

Now, we are ready to formulate and prove our solution to the equivalent Problems 1 & 2.

3.3 The generalized Breiman models↩︎

In this section, we discuss the broad class of so-called Breiman’s models, which realize all possible tail dependence measures. We will demonstrate that in the context of these models, the solutions to the variational problems in Theorem 2 yield asymptotically optimal homogeneous predictors. We begin with a sharp version of the classical Breiman’s lemma, which is a reformulation of Lemma 2.3 in [32].

In the case when \(\eta\) in ?? is replaced by a random vector, one obtains a rich family of multivariate regularly varying models [25], [33], which realize all possible asymptotic tail dependence regimes.

The next result shows a curious fact that the multivariate Breiman-type models are, in fact, regularly varying in a stronger total variation sense. This fact is of independent interest, which we will revisit in Section 4.

The proof of Proposition 5 is given in Section 8.3.3.
Pareto-Breiman models. In the rest of this section, we return to the optimal prediction context. We shall consider the important special case of the Breiman models, where the heavy-tailed factor has asymptotically Pareto tails, i.e., ?? holds and for simplicity \(\alpha=1\).

Specifically, we let \[\label{e:Y95X95Breiman} Z = (Y,X) = \xi\cdot (V,W),\tag{37}\] be given by a Breiman-type model as in Proposition 5 such that \(V\ge 0\), \((V,W) = (V,W_1,\cdots,W_d)\) and \(\xi\) are independent. We assume, moreover that \(\xi\) satisfies ?? with \(\alpha = 1\).

The proof of this result is given in Section 8.3.3. We shall comment next on the subtle difference between the distribution of \(\Theta\) in 38 the marginal angular distribution of \(X\).

The following corollary shows that in the context of the Pareto-Breiman models, Theorem 2 yields optimal homogeneous predictors over the entire class of all \(\tau_X\)-support dominated and asymptotically calibrated predictors.

4 Statistical inference↩︎

In this section, the focus is on the statistical inference of the optimal predictors. We will work in the established peaks-over-threshold framework (Section 4.2) and obtain general results for the universal extremal consistency for homogeneous predictors constructed in this framework (Section 4.3).

4.1 Setup and main ideas↩︎

Throughout this section, we let \((Y,X)\) be a jointly regularly random vector in \(\mathbb{R}_+\times \mathbb{R}^{d}\) in the sense of Assumption 1, where for simplicity, we suppose that the radial function is \[\tau(y,x):= y + \|x\|_1 = y + \sum_{i=1}^d |x_i|, \;\;(y,x)\in \mathbb{R}_+\times \mathbb{R}^d,\] with \(\|\cdot\|_1\) denoting the \(\ell_1\)-norm. In this case, with the notation in 10 and 11 , we have: \[\label{e:wt-U-and-wt-Theta-setup} \widetilde{U} = \frac{Y}{ Y + \|X\|_1},\quadand \quad (1-\widetilde{U}) \widetilde{\Theta} = \frac{X}{ Y + \|X\|_1 },\quadwhere\widetilde{\Theta}:= \frac{X}{\|X\|_1}.\tag{44}\]

For the purpose of conditional quantile inference, we shall assume that the convergence of the angular distributions in 12 holds in a stronger, total variation norm sense, i.e., \[\label{e:wt-U-wt-Theta-convergence} (\widetilde{U}, (1-\widetilde{U})\widetilde{\Theta}) \mid \{\tau(Y,X)>t\}\overset{TV}{\longrightarrow}\Theta_{(Y,X)}=:(U,(1-U)\Theta),\quad as \;t\rightarrow\infty.\tag{45}\]

Finally, we shall suppose that the joint limit distribution \(P_{(U,\Theta)}\) of \(U\) and \(\Theta\) can be written as \[\label{e:P95U44Theta-TV} P_{(U,\Theta)}(du,d\theta) = f_{U|\Theta}(u|\theta) du p_\Theta(d\theta),\tag{46}\] where the Lebesgue density of the conditional distribution \(U|\Theta\), is assumed to be positive, i.e., \[\label{e:f95U124Theta620} f_{U|\Theta}(u|\theta) >0\quad modulo\;\;dup_\Theta(d\theta).\tag{47}\]

Recall that \(q_\alpha(\theta)\), defined in 21 , stands for the left-continuous \(\alpha\)-quantile of the weighted conditional density \(\propto (1-u)f_{U|\Theta}(u|\theta)\). Then, 47 implies that \(\alpha \mapsto q_\alpha(\theta) \in (0,1)\) is continuous and strictly increasing in \(\alpha\in (0,1)\) modulo \(p_\Theta(d\theta)\). In view of Theorem 2, we then have that for all \(0<\alpha<1\), \[\label{e:g95alpha-optimal} g_\alpha(x):= \|x\| \cdot \frac{q_\alpha(x/\|x\|)}{1-q_\alpha(x/\|x\|)},\;\;x\in \mathbb{R}^d,\tag{48}\] yields an optimal homogeneous predictor. Specifically, \[\label{e:Lambda-optimal} \Lambda(g_\alpha) = \sup_{g \in D(g_\alpha)} \Lambda(g),\tag{49}\] where \(\Lambda(\cdot)\) is as in 13 . Recall that the optimization domain \(D(g_\alpha)\) consists of all measurable positive homogeneous functions such that \[C(g) = \mathbb{E}[ (1-U) g(\Theta) ] = C(g_\alpha),\] where \(C(\cdot)\) stands for the constraint functional in 18 . Finally, notice that the optimal value is achieved by a continuous predictor \((\alpha,x) \mapsto g_\alpha(x)\).

Given the above discussion, obtaining an asymptotically optimal homogeneous predictor of \(\{Y>t\}\) amounts to estimating the quantile functions \(q_\alpha(\cdot)\), and finding the value of \(\alpha \in (0,1)\) such that the predictor in 48 is calibrated. This, in view of the continuity of \(g_\alpha\), amounts to having \[\label{e:C40g95alpha4241} C(g_{\alpha^*}) = \mathbb{E}[U],\tag{50}\] because \(\mathbb{P}[g_\alpha(X)>t]/\mathbb{P}[Y>t] \to C(g_\alpha)/\mathbb{E}[U],\) as \(t\to\infty\), by Proposition 2. This is the main idea behind the optimal predictor inference methodology described in Sections 5 and 10. Theorem 3 (below) shows that this methodology will yield predictors with asymptotically optimal precision. To make this precise, however, one needs to address several details.

If one had iid samples \((U_i,\Theta_i),\;i=1,\cdots,n\) from the limiting angular distribution \(P_\infty:= P_{(U,\Theta)}\), then \(q_\alpha(\cdot)\) can be consistently estimated (as \(n\to\infty\)) in many ways. For example, one can use classical quantile regression with a weighted loss or by using weighted quantile random forest estimators [16], [34] (cf Section 5).

In practice, however, one does not have iid samples from the limiting angular distribution, but instead a triangular array of samples \((\widetilde{U}_i,\widetilde{\Theta}_i)\)’s coming from the distribution \(P_{t_n}\) of \((\widetilde{U},\widetilde{\Theta})\, |\, \{\tau(Y,X)>t_n\}\) indexed by a threshold \(t_n\to\infty\). Thus, to show that the quantile estimates \(\widehat q_\alpha\) based on these data are asymptotically consistent, one needs an additional step that controls the rate at which \(P_{t_n}\) converges to \(P_\infty\). This will be done by appealing to a Le Cam-type contiguity argument in Section 4.2. This step will subsequently allow us to establish general results that can be thought of as extremal counterparts to universal consistency (Theorem 3 in Section 4.3).

4.2 Peaks-over-threshold: A contiguity perspective↩︎

Let \(\{(Y_i,X_i)\}_{i=1}^n\) be independent copies of \((Y,X)\) and define the radii \[R_i := \tau(Y_i,X_i) = Y_i + \|X_i\|_1, \qquad i=1,\ldots,n.\] Given a threshold \(t=t_n\to \infty\), consider the associated index set of exceedances: \[\mathcal{I}_n := \{ i \in \{1,\ldots,n\} : R_i > t_n \}.\]

Following 44 , the angular components of these exceedance observations are denoted as \[\label{e:wt-U95i-wt-Theta95i} \widetilde{U}_i := \frac{Y_i}{R_i},\quadand \;\;\widetilde{\Theta}_i := \frac{X_i}{\|X_i\|_1}, \quad i \in \mathcal{I}_n.\tag{51}\] The Peaks-over-Threshold (PoT) methodology stipulates that by Relation 45 , as \(t\to\infty\), one can regard the \((\widetilde{U}_i, \widetilde{\Theta}_i),\;i\in {\cal I}_n\) as independent samples coming from \(P_{\infty}:= P_{(U,\Theta)}\). Thus, if \(t=t_n\to\infty\) and at the same time the number of exceedances \[K(n):= |{\cal I}_n| \to\infty, \;\; asn\to\infty,\] one expects to be able to consistently estimate functionals of \(P_\infty\) with statistics based on the sample \((\widetilde{U}_i,\widetilde{\Theta}_i),\;i\in {\cal I}_n\). To formalize this idea, one needs to address two issues: (i) The sample size \(K(n)\) is random; (ii) The conditional distribution \(P_{t_n}\) of \((\widetilde{U},\widetilde{\Theta})\) given \(\tau(Y,X)>t_n\) may approach \(P_\infty\) at a slow rate. The first issue is elegantly addressed via the so-called dècoupage de Lévy, while the second by Proposition 8 using the ideas of contiguity.

This result is known see e.g. [35] or Remark 22 in Section 8.4.1 for a proof. Proposition 7 shows that the exceedance “angles” \((\widetilde{U}_i,\widetilde{\Theta}_i)\) can be regarded as independent draws from the distribution \(P_t\). Moreover, the random number of exceedances \(K(n)\) and their values are independent, i.e., can be decoupled, even though they come from the common sample \(\{(Y_i,X_i),\;i\in [n]\}\).
The following result makes the intuition behind Peaks-over-Threshold formally precise. It uses a contiguity-type argument in the spirit of Le Cam’s First Lemma [36].

The proof can be found in Section 8.4.1 (cf Proposition 16). This result shows that provided ?? holds, any statistic based on the sample \(\{(\widetilde{U}_i,\widetilde{\Theta}_i),\;i\in {\cal I}_n\}\) has the same consistency behavior as the same statistic computed from \(\{(U_i,\Theta_i)\}_{i=1}^{K(n)}\) – an iid sample from the asymptotic angular distribution \(P_\infty = P_{(U,\Theta)}\). This general contiguity argument allows us to seamlessly extend consistency results for iid samples from the unattainable asymptotic distribution \(P_\infty\) to the triangular array type setting of peaks-over-threshold.

4.3 On the universal consistency in optimal homogeneous prediction↩︎

In this section, we provide a general consistency result for the optimal homogeneous predictors under rather mild conditions. In Section 5, we demonstrate that these conditions are fulfilled for a class of random forest-based predictors. These general consistency results can be viewed as an extreme value counterpart to the classical universal consistency theory of Stone [15].

To be precise, the focus here is on showing convergence to the optimal precision (equivalently, loss) rather than on the consistency of the predictors themselves. Namely, consistency results are stated through the precision functional \(\Lambda\) defined in equation 13 for functions \(g\) as \[\Lambda(g) = \frac{1}{\mathbb{E}[U]} \mathbb{E}[U \wedge (1-U)g(\Theta)] = \lim_{t\to\infty} \frac{\mathbb{P}[Y>t, g(X)>t]}{\mathbb{P}[Y>t]};\] Observe that if \(g(X)\) is calibrated then we have \(\Lambda(g) = \lambda(Y,g(X))\). Hence the aim is to provide asymptotic calibration conditions and consistency of \(\Lambda(\hat{g}_n)\) for a suitable sequence of predictors \((\hat{g}_n)_{n\geq 1}\) based on a growing training data set.

Given the quantile estimators in Assumption 3, define \[\label{e:wh95g95alpha44k} \widehat g_{\alpha, k} (x) := \|x\| \cdot \frac{\widehat q_{\alpha,k}(x/\|x\|)}{1-\widehat q_{\alpha,k}(x/\|x\|)},\;\;x\in \mathbb{R}^d,\tag{52}\] where \(\widehat g_{\alpha,k} (x):= 0\).

Now, consider the Peaks-over-Threshold methodology of Section 4.2. That is, let \((Y_i,X_i),\;i\in [n]\) be independent copies of a jointly regularly varying random vector \((Y,X)\) in \(\mathbb{R}_+\times \mathbb{R}^d\) and focus on the random sample of exceedance-based angles \(\{(\widetilde{U}_i,\widetilde{\Theta}_i),\;i\in {\cal I}_n\}\) in 51 , for a sequence of thresholds \(t=t_n\to\infty\). Let also \(\alpha_n\) be a sequence of estimators for \(\alpha^*\) in 50 and define the predictors \[\label{e:wh95g95n} \widehat g_n(x):= \widehat g_{\alpha_n, K(n)}(x),\tag{53}\] computed with the \((U_i,\Theta_i)\)’s replaced by the the \((\widetilde{U}_i,\widetilde{\Theta}_i)\)’s, where \(K(n) = |{\cal I}_n|\).

The following theorem shows the universal extremal consistency of the predictors \(\{\widehat g_{K(n)}(X)>t\}\) for the events \(\{Y>t\}\), as \(t\to\infty\) under mild conditions. That is, we will show that these event predictors are asymptotically calibrated and that their precision is asymptotically optimal.

We will first comment on the consequences of this result and then present its proof.

5 Methodology and numerical examples↩︎

5.1 Weighted quantile inference↩︎

Consider a random variable \(\xi\) whose conditional density given \(\Theta=\theta\) is defined by \[f_{\xi \mid \Theta=\theta}(u) \propto (1-u)f_{U\mid \Theta}(u\mid \theta),\] and assume that Assumption 3 holds.

Let \[\rho_\alpha(u) = \bigl(\alpha-\mathbb{1}_{\{u\le 0\}}\bigr)u, \qquad \alpha\in(0,1),\] denote the standard check loss function.

Under Assumption 3, the conditional density \(f_{\xi\mid \Theta=\theta}\) is continuous and strictly positive in a neighborhood of its \(\alpha\)-quantile \[q_{\alpha,\xi}(\theta) := F^{-1}_{\xi\mid\Theta=\theta}(\alpha), \qquad \theta\in\mathbb{S}^{d-1}.\] It is well known (see, e.g., [38]) that, for every fixed \(\theta\in\mathbb{S}^{d-1}\), the conditional quantile \(q_{\alpha,\xi}(\theta)\) is the unique minimizer of \[a \longmapsto \mathbb{E}\!\left[ \rho_\alpha(\xi-a) \,\middle|\, \Theta=\theta \right].\]

Since \(q_\alpha\) is assumed to be continuous, the conditional quantile admits the representation \[\label{e:q-argmin} q_\alpha(\theta) = \underset{a\in(0,1)}{\operatorname{argmin}} \, \mathbb{E}\!\left[ \rho_\alpha(\xi-a) \,\middle|\, \Theta=\theta \right].\tag{71}\]

Moreover, integrating with respect to the distribution of \(\Theta\) and using the definition of \(\xi\), we obtain \[\int_{\mathbb{S}^{d-1}} \mathbb{E}\!\left[ \rho_\alpha(\xi-a) \,\middle|\, \Theta=\theta \right] p_\Theta(d\theta) = \mathbb{E}\!\left[ \rho_\alpha(U-a)(1-U) \right].\]

Consequently, introducing \[M_\alpha(a) := \mathbb{E}\!\left[ \rho_\alpha(U-a)(1-U) \right],\] we have \[q_\alpha(\theta) = \underset{a\in(0,1)}{\operatorname{argmin}} \, M_\alpha(a).\]

Since \(M_\alpha\) is convex and subdifferentiable, its minimizer is characterized by the first-order optimality condition \[\alpha = \frac{ \mathbb{E}\!\left[ \mathbb{1}_{\{U\le a\}}(1-U) \,\middle|\, \Theta=\theta \right] }{ \mathbb{E}\!\left[ 1-U \,\middle|\, \Theta=\theta \right] }.\] The right-hand side is precisely the conditional distribution function of \(\xi\) given \(\Theta=\theta\). This observation leads to the following weighted quantile estimator [39], [40]: \[\label{e:q-hat} \widehat q_\alpha(\theta) := \widehat F^{-1}_{k,\xi\mid\Theta=\theta}(\alpha),\tag{72}\] where \[\widehat F_{k,\xi\mid\Theta=\theta}(a) = \sum_{i=1}^{k} \frac{ \omega_h(\Theta_i,\theta)(1-U_i) }{ \sum_{j=1}^{k} \omega_h(\Theta_j,\theta)(1-U_j) } \, \mathbb{1}_{\{U_i\le a\}}.\]

Here \((U_i,\Theta_i)_{i=1}^{k}\) are independent copies of \((U,\Theta)\), and \[\omega_h(\theta,\theta') = \omega\!\left( \frac{d(\theta,\theta')}{h} \right),\] with \[d(\theta,\theta') = \arccos(\theta^\top\theta')\] denoting the geodesic distance on the sphere \(\mathbb{S}^{d-1}\) and \(h\) a bandwidth depending on \(k\).

For a fixed \(\theta\in\mathbb{S}^{d-1}\), consistency of the estimator \(\widehat q_\alpha(\theta)\) follows, as long as \(kh(k)^{d-1}\) goes to \(0\) as \(n\) goes to infinity, from standard arguments because the density is strictly positive in kernel estimation; see, for instance, Section 1.2 of tsybakov2009? for kernel estimator and [35] Proposition 0.1 for the convergence of inverse functions. Note that a wide range of weights functions \(\omega\) can be chosen in practice. As an example, implementation with quantile regression forests [16] are discussed in Supplements 10 and 11.

However, pointwise consistency alone is not sufficient to guarantee the uniform convergence required in Assumption 3. The estimator \(\widehat q_\alpha\) can nevertheless be modified to obtain a uniformly consistent version.

An explicit construction of the estimator \(\widehat q_{n,\alpha}(\cdot)\) is provided in the proof of Proposition 17. Note that the sequence \(m(n)\) therein is not very explicit since one only requires the pointwise consistency of \(q_{n,\alpha}(\theta)\). If one has more information on the rate of this consistency, then a conservative choice of \(m(n)\) can be obtained. For more details, see Remark 23.

5.2 Examples and simulations↩︎

In this section, we illustrate the performance of our optimal homogeneous predictors that were estimated from data in comparison with the oracle predictors in several simulated models. We implemented the peaks-over-threshold methodology described in Section 4.1. While the quantile \((\theta,\alpha)\mapsto q_\alpha(\theta)\) involved therein can be estimated using a variety of methods, we opted for the versatile nonparametric method of [16] known as quantile regression forest. This method works remarkably well in multiple dimensions and is easy to tune. Specifically, we implemented our estimator in R [18], [19] using a state-of-the-art computationally and memory efficient implementation of quantile regression forests in the R package ranger by [41]. For more details, see Section 10.

We present two types of models: spectrally discrete and spectrally continuous ones, that illustrate and the general case described in Section 4. For more details and proofs, see Section 9.
The spectrally discrete case. Consider the linear factor model \[\label{e:linear-model-preview} Y = \sum_{i=1}^{p} b_i \xi_i \quad \text{and} \quad X = \sum_{i=1}^{r} a_i \xi_i,\;\;(r\le p),\tag{73}\] where \(\xi_1,\ldots,\xi_p\) are independent standard \(1\)-Pareto random variables, viewed as latent factors; \(b_i \ge 0\), and \(a_i \in \mathbb{R}^d\). Then, vector \(Z=(Y,X)\) is regularly varying with a discrete angular measure concentrated on the \(p\) directions \(c_i/\|c_i\|_1\), where \(c_i=(b_i,a_i)\) (see Proposition 18 in Section 9.1).

Under mild non-proportionality assumptions on the \(a_i\)’s (Corollary 4), the optimal asymptotic precision takes the explicit form \[\label{e:lambda-opt-discrete-preview} \lambda^{(\mathrm{opt})}_{\mathcal{G}}(Y,X) = \frac{\sum_{i=1}^{r} b_i}{\sum_{i=1}^{p} b_i},\tag{74}\] where \(r \le p\) is the number of factors with \(a_i \ne 0\).

Given the \(b_i\)’s and \(a_i\)’s, one can construct a continuous, non-negative, \(1\)-homogeneous function \(h^{(\mathrm{opt})}\) achieving 74 , which we refer to as the oracle predictor (see 137 ). In practice the \(b_i\)’s and \(a_i\)’s are unknown. Nevertheless, the spectral measure of \((Y,X)\) in 136 can be estimated consistently from data by using the peaks-over-threshold methodology combined with a non-parametric estimator of the quantiles of \(q_\alpha(\theta)\) the tilted conditional angular distribution (cf Algorithm 4).

Figure 1 illustrates the finite-sample behavior of our nonparametric estimator across repeated experiments: The left panel corresponds to the under-determined case \(r < p\), where 74 is strictly less than one, while the right panel shows the perfect precision case \(r = p\), where \(\lambda^{(\mathrm{opt})}_{\mathcal{G}} = 1\). The estimators are based on a training sample \(\{(Y_i,X_i)\}\) of \(n_{\rm train}= 10^4\) using a random threshold based on the \(0.95\)th empirical quantile. The empirical precision is then computed over an independently generated testing sample \(\{(Y_i^*,X_i^*)\}_{i=1}^{n_{\rm test}}\) of size \(n_{\rm test}=10^4\), where \[\label{e:wh-lambda40p41} \widehat \lambda(p):= \frac{1}{1-p} \sum_{i=1}^{n_{\rm test}} I\{ \widehat g(X_i^*) > \widehat F_{\widehat g}(p) \} I\{ Y_i^* > \widehat F_{Y}(p) \},\;\;p\in (0,1),\tag{75}\] where \(\widehat F_{\widehat g}\) and \(\widehat F_Y\) are the empirical CDFs of \(\{\widehat g(X_i^*)\}\) and \(\{Y_i^*\}\), respectively. The boxplots are based on \(100\) independent replications of the \(\widehat \lambda(p)\)’s. The coefficients \(b_i\) and the components of the vectors \(a_{i}\) in 73 were fixed across different replications, but drawn independently and at random from the Uniform\((0,1)\) distribution, therefore ensuring that \(a_i/\|a_i\| \not = a_j/\|a_j\|\), for all \(i\not = j\), with probability one.

Figure 1 shows a close agreement of the estimated optimal homogeneous predictor and the oracle, and both are rather close to the optimal homogeneous extremal precision (dashed horizontal line). It is remarkable that the empirical tail dependence coefficient is nearly one for thresholds \(p\) as low as \(0.85\) in the asymptotically perfect precision scenario.
The Pareto-Dirichlet case. As an illustration of the spectrally continuous setting, consider \(Z = (Y,X) = \xi\cdot(V,W)\), where \(\xi\) is a standard \(1\)-Pareto random variable independent of \((V,W) \sim \mathrm{Dirichlet}(\beta)\) with parameter \(\beta = (\beta_0,\beta_1,\ldots,\beta_d) \in (0,\infty)^{d+1}\). We call this the Pareto-Dirichlet\((\beta)\) model (cf Proposition 20). It is \(\tau\)-regularly varying with respect to the radial function \(\tau(y,x) = |y| + \|x\|_1\), and its angular measure is the Dirichlet distribution on the unit simplex \(\Delta_d\).

The neutrality property of the Dirichlet distribution [42] implies that \(U := V\) and \(\Theta := W/\|W\|_1\) are independent, so the conditional quantile \(q_\alpha(\theta)\) in Theorem 2 is constant in \(\theta\). This yields a remarkably simple optimal predictor (Proposition 20 in Section 9): \[\label{e:h-opt-Dirichlet-preview} h^{(\mathrm{opt})}(X) = c\,\|X\|_1, \qquad c := \frac{\beta_0}{\sum_{i=1}^{d}\beta_i},\tag{76}\] with optimal precision \[\label{e:lambda-opt-Dirichlet-preview} \lambda^{(\mathrm{opt})}_{\mathcal{G}}(Y,X) = \mathbb{E}\!\left[\frac{U}{\mu_U} \wedge \frac{1-U}{1-\mu_U}\right], \qquad U \sim \mathrm{Beta}\!\left(\beta_0,\,\textstyle\sum_{i=1}^d \beta_i\right), \quad \mu_U = \mathbb{E}[U].\tag{77}\] Moreover, \(\lambda^{(\mathrm{opt})}_{\mathcal{G}}(Y,X) = \lambda_p(Y, h^{(\mathrm{opt})}(X))\) for all \(p > p_0 := \mu_U \vee (1-\mu_U)\), so the precision is exactly constant beyond a threshold. This is clearly observed in the right panel of Figure 2, where the empirical tail-dependence coefficient \(\hat{\lambda}_p\) boxplots for both the oracle and the nonparametric estimator plateau at \(\lambda^{(\mathrm{opt})}_{\mathcal{G}}\) across a wide range of \(p\) values (see also Figure 3).

Figure 2 (left panel) shows boxplots of the empirical tail-dependence coefficients as in 75 based on \(100\) replications of an estimated optimal predictor and the corresponding oracle. As in Figure 1, we have \(n_{\rm train} = n_{\rm test} = 10^4\) and a random peaks-over-threshold of the \(0.95\)-th percentile. The Pareto-Dirichlet model involves coefficients \(\beta_0 = 1\) and \(\beta_i = (i+1)/10,\;i=1,\cdots,9\). As in the spectrally discrete case, we have a remarkably close agreement in the precision of the estimated predictor and the oracle, which in this case is simply proportional to \(\|X\|_1\). Both agree closely with the asymptotically optimal precision as indicated above.

The right panel of Figure 2 shows the same type of experiment where \((Y,X)\) come from a Gumbel copula, where \(d=10\). By Corollary 2, in this case the optimal homogeneous predictor is in fact the ultimate optimal predictor for all values of \(p\in (0,1)\). In this scenario, the non-parametric estimator is remarkably close to the oracle.

These brief simulation results illustrate the success of the optimal homogeneous prediction methodology. The quantile regression forest based implementation is nearly tuning–free and quite adaptive – working well in both spectrally discrete and continuous settings. We have implemented this methodology in the context of the challenging open problem of Solar flare prediction, where we found it to work out-of-the-box and provide state-of-the-art prediction scores for the most extreme X-class flares (cf Figure 6). We plan to pursue more applied and methodological directions of this research in a future work, but the preliminary findings are available in Section 11. The code used to produce all simulation in the paper is freely available from the GitHub repository [18]. An R Shiny app illustrating the methodology over a challenging Solar flare prediction problem is deployed on [19].

6 Discussion and future work↩︎

In this paper, we considered the fundamental problem of optimal extreme event prediction in terms of covariates. Starting from a general the Neyman-Pearson-type characterization of the optimal predictors (Theorem 1), we formulated an optimality framework that is closely related to the notion of tail-dependence. This naturally led us to explore the problem of optimal extremal prediction in the context of jointly regularly varying response and covariates. For the general class where the predictors are positive homogeneous of the covariates, the problem was cast into a variational problem – a constrained optimization problem for a convex integral functional with respect to the angular measure. Our first main contribution was an explicit solution to this problem (Theorem 2), which can be represented in terms of the quantiles of a certain tilted conditional distribution arising from the angular measure. This abstract result let us to concrete and yet general inference methodology for the optimal homogeneous predictors. Our second contribution was to establish the consistency for estimators of these optimal predictors in the framework of peaks-over-threshold (Theorem 3). This result, modulo mild regularity conditions, may be viewed as an extremal counterpart to the celebrated Stone’s Universal Consistency Theorems for the class of homogeneous predictors. Last but not least, using established quantile regression forests tools, the estimators were implemented into a broadly applicable and nearly tuning free statistical software, which was shown to achieve the oracle performance over simulated examples.

The problem of extreme event prediction is fundamental and ubiquitous. It can be seen as a series of increasingly imbalanced classification problems, which is rather challenging. Notably, [9] have identified the shortcomings of the traditional empirical risk minimization approach and proposed a solution. While motivated by similar problems, our approach is different. Here, by conditioning on the regime of extremes, the loss function is directly related to tail-dependence, which allowed us to take a different look through the lens of calculus of variations.

We provide a rather complete solution to the optimal extremal prediction with homogeneous functions in the context of jointly regularly varying \(Y\) and \(X\). While in many cases the class of homogeneous predictors do contain the optimal extremal predictors (cf Example 1 and Corollary 2), this is not always the case. Indeed, as shown in Example 8, even under joint regular variation, the ‘extremes’ of \(X\) do not necessarily predict the extremes of \(Y\). To the best of our knowledge, the problem of finding universally consistent predictors in the extreme value sense, remains open both in terms of its precise formulation and solution.

Acknowledgments↩︎

We thank Victor Verma for providing curated the GOES flux time series and SHARP parameter data for the Solar flare prediction experiment in Section 11. SS was partially supported by the NSF grant CNS/CSE-2319592 “Collaborative Research: IMR: MM-1A: Scalable Statistical Methodology for Performance Monitoring, Anomaly Identification, and Mapping Network Accessibility from Active Measurements”.

7 Examples↩︎

7.1 Examples of optimal predictors↩︎

In this section, we elaborate on the closed-form optimal predictors discussed in Section 2. Specifically, we prove the stated characterizations of the optimal predictors in the bi-variate max-stable and Archimedean copula models with completely monotone generators.

7.1.1 Bivariate extreme value distributions↩︎

Suppose that \((Y,X)\) follow the a bivariate extreme value distribution \(G\), where without loss of generality the marginal distributions of \(X\) and \(Y\) are standard unit Fréchet. That is, \[F_X(u) = \mathbb{P}[X\le u] = F_Y(u) = \mathbb{P}[Y\le u] = \exp\{-1/u\},\;u>0.\] It is known that [35], for all \(x,y\ge 0,\) \[\label{e:G-EVD} G(x,y) = e^{-\mu(A_{x,y})},\tag{78}\] where \(A_{x,y} = [0,\infty)^2 \setminus ([0,x]\times[0,y])\) and \(\mu\) is a Borel measure on \([0,\infty)^2\setminus \{(0,0)\}\), such that \(\mu(A_{x,y})<\infty\) provided \(x>0\) or \(y>0\). It can be shown that \(\mu(t A) = t^{-1} \mu(A)\), for all \(t>0\) and Borel sets \(A\subset [0,\infty)^2\setminus \{(0,0)\}\) and hence for all \(x,y \ge 0\) \[\label{e:G-EVD-spectral} G(x,y) = \exp\Big\{ - V(1/x,1/y) \Big\},\;\; where V(a,b):= \int_0^1 \max\{ w a, (1-w) b)\} \sigma(dw)\tag{79}\] where \(\sigma\) is a finite measure on \([0,1]\) and by convention \(1/\infty = 0\) and \(1/0 = \infty\). The function \(V\) is referred to as the stable tail-dependence function of the bivariate extreme value distribution.

The following result is rather natural and it shows that the events \(\{X>x_0\}\) are optimal predictors of \(\{Y>y_0\}\).

7.1.2 Archimedean copula: Proofs for Section 2.2↩︎

We start with the following auxiliary result.

Proposition 1 implies the following result, which is of independent interest since it concerns an important class of multivariate max-stable models.

8 Proofs↩︎

8.1 Multivariate regular variation↩︎

In this section, we provide a brief review of the theory of multivariate regular variation of Borel measures on \(\mathbb{R}^k,\;k\ge 1\). We consider the general case of \(\tau\)-regular variation, where \[\label{e:tau} \tau:\mathbb{R}^k \to \mathbb{R}_+\tag{89}\] is a fixed non-negative, continuous \(1\)-homogeneous function [26], [27]. For an abstract and general theory of regular variation with respect to the notion of a bornology, see e.g., [28] and the references therein.

8.2 The space \(\mathbb{M}_\tau\) and the Tail Index Theorem↩︎

Let \(\mathbb{M}_\tau = \{ \mu\, :\, \mu(\tau >\epsilon)<\infty,\; for all \epsilon>0\}\) be the collection of Borel measures \(\mu\) on \(\mathbb{R}^k\setminus\{\tau=0\}\) such that \(\mu(\{\tau>\epsilon\})<\infty\), for all \(\epsilon>0\). Such measures assign finite masses to all Borel sets \(A\) that are bounded away from \(0\), relative to the \(\tau\)-loss, i.e., \(\tau(A):=\inf_{x\in A}\tau(x)\) is positive. We shall denote this class of Borel sets by \({\cal B}_0(\tau)\) and \({\cal B}(\tau) = \sigma({\cal B}_0(\tau))\) is the \(\sigma\)-field of all Borel sets in \(\mathbb{R}^k\setminus\{\tau=0\}\).

If \(\mu_n,\mu\in \mathbb{M}_\tau\), then we write \(\mu_n\stackrel{M_0}{\to} \mu\), as \(n\to\infty\), if \[\label{e:M0-convergence-via-f} \int f(x)\mu_n(dx) \to \int f(x)\mu(dx),\;\; as n\to\infty,\tag{90}\] for all bounded and continuous functions \(f\) such that \(\{|f|>0\} \subset \{\tau>\epsilon\}\) for some \(\epsilon>0\). It can be shown that \(\mu_n\stackrel{M_0}{\to} \mu\) if and only if \[\mu_n(A)\to \mu(A),\;\; as n\to\infty,\] for all \(A\in {\cal B}_0(\tau)\) such that \(\mu(\partial A) =0\), where \(\partial A:=\overline{A}\setminus A^\circ\) denotes the topological boundary of the set \(A\) [45]. The above-defined \(M_0\)-convergence of measures in \(\mathbb{M}_\tau\) is generated by a metric, which renders \(\mathbb{M}_\tau\) a complete separable metric space [45].

The proof of this result can be derived from many excellent treatments in the literature. Specifically, see, for example, Theorem 3.1 and Corollary 4.4 in [46], as well as [45], [47], the monographs [25], [29], the general theory in [28], and the references therein for more details.

8.3 Proofs for Section 3↩︎

8.3.1 Proofs for Section 3.1↩︎

For completeness, we will present the proof of Proposition 2 below although similar results can be found in Proposition 2.5 of [23] and Theorem 2.1 of [30]. The proof consequence of the following characterization of multivariate regular variation [25].

8.3.2 Proofs for Section 3.2↩︎

8.3.3 Proofs of Section 3.3↩︎

8.4 Proofs of Section 4↩︎

In this section we provide proofs and auxiliary results needed for Theorem 3. We start by providing a contiguity argument, which will be used to formally justify the Peaks-over-Threshold methodology. This argument is of independent interest since it allows one to effectively plug-in a sample from the limit angular distribution \(\sigma\) in place of the triangular array coming from the exceedance distributions for a sequence of growing thresholds.

8.4.1 More on contiguity and the dècoupage de Lévy↩︎

We begin by recalling established results and notation on Hellinger and total variation distances of measures. Let \(P\) and \(Q\) be two probability measures, defined on a common measure space \((S,{\cal S})\). Recall that the squared Hellinger distance of \(P\) and \(Q\) is defined as: \[H^2(P,Q):= \frac{1}{2} \int_S \Big(\sqrt{\frac{dP}{d\nu}} - \sqrt{\frac{dQ}{d\nu}}\Big)^2 d\nu,\] where \(\nu\) is any other measure dominating \(P\) and \(Q\), i.e., \(\nu = P + Q\).

We shall denote by \(P^{\otimes k}\) the product measure on the cartesian product space \((S^k,\sigma({\cal S}^k))\). Note that \(P^{\otimes k}\) is the probability distribution of a sample of \(k\) independent elements taking values in \(S\) and having probability distribution \(P\). The following inequalities are well known [36], [50], [51] \[\label{e:Hellinger-product-inequality} H^2( P^{\otimes k}, Q^{\otimes k}) \le k H^2(P,Q)\tag{123}\] and \[\label{e:Hellinger-TV-inequality} H^2(P,Q) \le TV(P,Q) \le \sqrt{2} H(P,Q),\tag{124}\] where \(TV(P,Q)\) stands for the total variation distance: \[\label{e:TV-def} TV(P,Q):= \sup_{A\in {\cal F}} |P(A)- Q(A)|.\tag{125}\]

Let now \(P_t\) be the conditional distribution of \(Z/\tau(Z)\, |\, \tau(Z)>t\) and let \(P_\infty\) denote its limit.

The next result is similar in spirit to Le Cam’s First Lemma (see also Remark 21 below).

The following remark shows why the assumption of independence of \(K(t)\) and the \(\Theta_i\)’s in Proposition 16 is not restrictive in the setting of peaks over threshold.

8.4.2 On the construction of uniformly consistent estimators↩︎

In this section we demonstrate that a point-wise consistent estimator of a continuous function can be modified to obtain a uniformly consistent estimator over compacta. Although surprising, at a first glance, this fact is akin to the method of sieves and is a direct consequence of the principle of the diagonal. The caveat is that this modification may ultimately converge at a potentially slow rate. We were unable to find a reference for this simple fact and therefore we give it here for the sake of completeness.

We begin with the following

9 Examples of asymptotically optimal homogeneous predictors↩︎

In this section, we will illustrate the optimal homogeneous predictors of Theorem 2 for two classes of spectral measures – purely discrete and absolutely continuous. The spectrally discrete case can be understood in a simpler way and yet, it is interesting since a curious perfect asymptotic precision phenomenon emerges in this context. On the other hand, the spectrally continuous case is illustrated with a simple and yet elegant Dirichlet-type model, which may be of independent interest.

9.1 The spectrally discrete case: A perfect precision phenomenon↩︎

Consider the linear independent factor model: \[\label{e:linear-model} Y = \sum_{i=1}^p b_i \xi_i\;\;\; and \;\;\; X = \sum_{i=1}^p a_i \xi_i,\tag{135}\] where the \(\xi_i\)’s are independent standard \(1\)-Pareto random variables, i.e.,\(\mathbb{P}[\xi_i> x] = 1/x,\;x\ge 1\).

We suppose that \(b_i\ge 0\) and \(a_i\in \mathbb{R}^d,\;i=1,\cdots,p,\) where not all \(b_i\)’s and not all \(a_i\)’s are zero. Letting \(\tau(y,x):= y_+ + \|x\|\) for some (any) fixed norm \(\|\cdot\|\) on \(\mathbb{R}^d\), it follows by Proposition 7.3 in [52] that \(Z:= (Y,X) \in RV_1(\mathbb{R}_+\times \mathbb{R}^d,\{a_n^{(Z)}:=n\},c_Z,\tau,\sigma)\), where with \(c_i:= (b_i,a_i)\), we have \[c_Z=\sum_{i=1}^p \tau(c_i) = \sum_{i=1}^p b_i+\|a_i\|.\] In this case, the angular measure \(\sigma\) is discrete and takes the form: \[\label{e:p:linear-model} \sigma(d\theta) = \frac{1}{c_Z} \sum_{i=1}^p \tau(c_i) \delta_{c_i/\tau(c_i)}(d\theta).\tag{136}\]

We begin with a counterpart to Proposition 3 and for the sake of completeness, we provide a proof.

The following result as an immediate but rather curious consequence of the formula in ?? .

9.2 The spectrally continuous case↩︎

In this section, we illustrate Theorem 2 in the so-called spectrally continuous case where the conditional distribution \(p(du |\theta)\) has a density. Inference methodology for the optimal predictors in this rather general regime will be developed in Section 4.

We consider here non-negative random vectors \(Z=(Y,X)\) that are \(\tau\)-regularly varying with \[\tau(y,x) := y_+ + \sum_{i=1}^d (x_i)_+.\] In this case the unit sphere \(S_\tau\) becomes unit simplex: \[S_\tau = \Delta_d:= \Big\{ w = (w_i)_{i=0}^d\, :\, \sum_{i=0}^d w_i= 1,\;w_i\ge 0,\;i=0,\cdots,d\Big\}.\] A natural family of probability distributions on \(\Delta_d\) are the Siri holet models. Recall that a random vector \(W = (W_0,W_1,\cdots,W_d)\) has the Dirichlet distribution with parameter vector \(\beta = (\beta)_{i=0}^d \in (0,\infty)^{d+1}\) if \(W\) takes values in \(\Delta_d\) and the joint density of \(W_1,\cdots,W_d\) is given by: \[f_{W_1,\cdots,W_{d}}(x_1,\cdots,x_{d}) = \frac{\Gamma(\sum_{i=0}^d\beta_i)}{\prod_{i=0}^d \Gamma(\beta_i)}\cdot \Big(1-\sum_{i=1}^d x_i\Big)^{\beta_0-1}\times \prod_{i=1}^{d} x_i^{\beta_i-1},\;\;x_i \ge 0,\;\sum_{i=1}^d x_i\le 1,\] In this case, we shall write \(W\sim {\rm Dirichlet}(\beta)\).

Thus, in this context, it is natural to consider \(\tau\)-regularly varying vectors \(Z\) with Dirichlet angular distribution [26], [53]–[55].

The optimal homogeneous predictors in the class of spectrally Dirichlet model are particularly elegant.

Recall that a non-negative random vector \[(V,W) = (V,W_1,\cdots,W_d),\;d\ge 1\] has the multivariate Dirichlet distribution with parameter vector \(\beta = (\beta_0,\beta_1,\cdots, \beta_d)\in (0,\infty)^{d+1}\) if: (i) \(V + \sum_{i=1}^d W_i = 1\) and (ii) the joint probability density of \(W_1,\cdots,W_d\) is given by: \[f_{W_1,\cdots,W_{d}}(x_1,\cdots,x_{d}) = \frac{\Gamma(\sum_{i=0}^d\beta_i)}{\prod_{i=0}^d \Gamma(\beta_i)}\cdot \Big(1-\sum_{i=1}^d x_i\Big)^{\beta_0-1}\times \prod_{i=1}^{d} x_i^{\beta_i-1},\;\;x_i \ge 0,\;\sum_{i=1}^d x_i\le 1,\] In this case, we shall write \((V,W)\sim {\rm Dirichlet}(\beta)\).

Since the Dirichlet distribution is a probability distribution on the positive unit simplex, it is natural to consider polar coordinates with respect to the \(L^1\)-norm: \[\tau(y,x) := |y| + \|x\|_1,\] and define a version of the Pareto-Breiman models in Section 3.3.

The Pareto-Dirichlet models and their spectral mixtures have been considered due to the tractable analytical form of its spectral distribution [26], [53]–[55]. The following result shows that the optimal homogeneous predictor takes a particularly simple form for the Pareto-Dirichlet model.

Figure 3 (left) illustrates this optimal extremal precision as a function of \(a_0=\beta_0\) and \(a_1=\|\widetilde{\beta}\|_1\), where the expectation in 139 is expressed via incomplete Beta functions in terms of \(\beta_0\) and \(\|\beta\|_1\).

The right panel of Figure 3 illustrates the finite-sample performance of the optimal (oracle) predictor and a general estimator for the optimal homogeneous predictor implemented in Section 10.

9.3 On the sub-optimality of homogeneous predictors↩︎

In some cases, the optimal predictors are in fact homogeneous functions of the covariates (cf Example 4). We will demonstrate next, however, that the class of homogeneous functions does not always contain neither optimal nor asymptotically optimal predictors. The spirit of this example is to show that the extremes of \(X\) do not necessarily predict well the extremes of \(Y\).

10 Practical Implementation↩︎

Many of the existing implementations of quantile regression can accommodate importance weights and thus achieve estimators for the quantiles of the tilted distribution in 21 . Here, we opted for using the versatile quantile regression forest of [16] and its efficient implementation in the package ranger of [41]. Specifically, this can be achieved with the R commands: \[\begin{align} &\tt weights <- 1 - u\\ &{\tt data<-data.frame(u=u, theta)}\\ &\tt rf<-ranger(u\sim.,data=data, case.weights=weights, quantreg=TRUE) \end{align}\] Here \({\tt u}\) and \({\tt \theta}\) are \((n_r\times 1)\) and \((n_r\times d)\) arrays, whose \(i\)th rows correspond to \(U_i\) and the vector \(\Theta_i\), respectively, \(i=1,\cdots,n_r\). The ranger function returns a random forest S3 object, which can be used to generate predictions of \(\widehat q_\alpha(\theta_{\rm new})\) using the commands \[\begin{align} &{\tt theta\_new = x\_new/L1\_norm(x\_new)}\\ &{\tt pred <- predict(rf, data=theta\_new, type="quantiles", quantiles=alpha)}\\ &{\tt q\_alpha <- pred\$predictions} \end{align}\]

The following Algorithm 4 summarizes the inference. The R code implementing the estimator using the ranger package of [41] is given in [18]. An R Shiny App demonstrating the estimator is deployed on [19], where further covariate selection heuristics are implemented.

Note that the sequence \((\alpha_n)_{n\in \mathbb{N}}\) provided by Algorithm 4 fulfill the assumption of theorem 3 by application of lemma 14.

11 An application to Solar flare prediction↩︎

In this section, we apply the general optimal homogeneous prediction methodology described in the previous section to an important open problem from space physics [4], [20].

In the following, the responses \(\{Y(t),\;t=1,2,\cdots\}\) come from a curated X-ray flux time series derived from a collection up to 15 different GOES satellite measurements (Figure 5). The covariates \(X_1(t),\cdots,X_d(t),\;t=1,2,\ldots\) are time series of \(d=22\) statistics referred to as SHARP time series. They are derived from solar disk images and capture fundamental physical characteristics of the active regions (sun spots) which are believed to drive the formation of solar flares [56], [57].

The objective is to predict extreme flux events \(\{Y(t+h)>y_0\}\) at some future time period \(t+h\), based on \(\{X_i(t),\;i=1,\cdots,d\}\). The events when the X-ray flux exceeds high thresholds are known as Solar flares and they are classified in several severity classes [58]. We will focus on predicting the most extreme M- and X-class flares, which correspond to X-ray fluxes exceeding \(10^{-5}\;W/m^2\) and \(10^{-4}\;W/m^2\), respectively (Figure 5).

Standardization. For both the response and the predictors, we consider non-overlapping block maxima over a window \(w\) corresponding to \(24\) hours: \[Y^{(w)} (t):= \max_{ w(t-1)+1 < s\le wt} Y(s)\;\; and \;\;X_{j}^{(w)}(t):= \max_{w(t-1)+1 < s\le wt} X_j(s).\] We then construct estimates \(\widehat F_{Y^{(w)}}\) and \(\widehat F_{X_j^{(w)}}\) of \(Y^{(w)}(t)\) and \(X_j^{(w)}(t)\) based on a training sample \(t\in {\cal D}_{\rm train}\). Tail extrapolation is used to ensure that the CDFs have unbounded supports and define the training sample with standardized approximately \(1\)-Pareto marginals: \[\widetilde{Y}(t) := \frac{1}{1-\widehat F_{Y^{(w)}}(Y^{(w)}(t))} \;\; and \;\; \widetilde{X}_j(t) := \frac{1}{1-\widehat F_{X_j^{(w)}}(X_j^{(w)}(t))},\] for all \(t\in {\cal D}_{\rm train}\), and \(j=1,\cdots,d\).

Then, under the assumption that the so-standardized vectors \((\widetilde{Y}(t+h), \widetilde{X}_1(t),\cdots, \widetilde{X}_d(t)), t, t+h\in {\cal D}_{\rm train}\) are realizations from a jointly regularly varying vector \((Y,X)\), we apply the inference methodology of Section 4. Specifically, we select the threshold \(r\) as a high quantile of the sample of norms, and apply the quantile random forest to estimate the optimal homogeneous predictor \(\widehat h\).
Prediction. Given a test sample \({\cal D}_{\rm test}\), we compute \[\widehat Y(t+h) := \widehat h( \widetilde{X}(t)),\;t, t+h\in {\cal D}_{\rm test},\] where \(\widetilde{X}(t) = (\widetilde{X}_j(t),\;j=1,\cdots,d)\) are standardized as: \[\widetilde{X}(t): =\frac{1}{\widehat F_{X_j^{(w)}} (X_j^{(w)}(t))},\;t\in {\cal D}_{\rm test}.\] Recall that \(\widehat F_{X_j^{(w)}}\) is an estimate of the CDF of \(X_j^{(w)}\) obtained from the training sample.

We predict a certain type of flare whenever \(\widehat Y(t+h)\) exceeds a threshold. Specifically, we let \(\hat{p} := \widehat F_{Y^{(w)}}(u)\) and produce the binary predictors: \[I_t(h):= I\{\widehat Y(t+h) > 1/(1-\widehat p) \}.\] This choice, under the assumption of stationarity, ensures the approximate calibration \(\mathbb{E}[ I_t(h)] \approx \mathbb{P}[ Y^{(w)}(t+h)> u]\). More precisely, for the M- and X-class flares, we use the thresholds \(u=10^{-5}W/m^2\) and \(u = 10^{-4} W/m^2\), respectively. The X-class flares are among the most extreme and rare events, which can lead to major magnetic storms. While less severe, the M-class flares can also cause major disruption of Earth’s ionosphere and lead to satellite and electric outages.

Results. Figure 6 shows a simple, nearly tuning-free application of the optimal homogeneous prediction methodology to Solar flare prediction. Specifically, the homogeneous predictor was trained over the period (2010-2017) and applied to forecasting M- and X-class flares over the period (2017-2025). Confusion matrices for the binary forecasts were used to obtain various metrics displayed on the plots. Specifically, letting TP, FP, TN, FN denote true-positive, false-negative, true-negative, and false-negative counts of the predictions, we obtain the standard metrics: \[{\rm prec} = \frac{TP}{TP+FN}, \;\;\text{TSS} = \frac{\mathrm{TP}}{\mathrm{TP} + \mathrm{FN}} - \frac{\mathrm{FP}}{\mathrm{FP} + \mathrm{TN}}, \; and \; {\rm missed} = \frac{FN}{TP+FN}.\] These metrics depend on the chosen quantile \(\tau_p\) for the variable \(\widehat Y=h(\widetilde{X})\), which ultimately determines the alarm rate: \[{\rm alarm} = \frac{{\rm TP} + {\rm FP }}{{ n_{\rm test}}},\] where \(n_{\rm test}\) is the size of the test sample.

As the alarm rate decays, generally, the precision increases but so does the proportion of missed events. Observe that the popular TSS statistic can be non-informative or at best misleading. Indeed, very high values of TSS may be associated with quite poor precision and hence quite high false-alarm rate. On the other hand, well-calibrated predictors with alarm rate equal to the event rate may have quite low TSS values. Nevertheless, the peak of the TSS values seem to be aligned with an elbow in the missed rate, so it happens to be a good heuristic in this case.

In principle, however, in addition to TSS, we advocate for examining the precision and missed metrics, as a function of the alarm rate. The practitioner may choose to settle for an acceptable alarm rate that achieves good precision and low missed rate that may or may not correspond to high TSS. Optimizing a flare prediction method for achieving high TSS without considering calibration or alarm rates could be misleading. This said, the TSS values achieved by the generic optimal homogeneous predictor methodology are surprisingly high. While the best M-class flare forecasting machines claim to achieve TSS rates as high as \(0.7\) (or sometimes \(0.8\) in carefully designed studies), our method reaches \(0.66\) with hardly any tuning at all, should one look at a suitable alarm rate of about \(30\%\). At the same time the existing X-class flare forecasting methods barely reach TSS as high as \(0.5\) in operational settings (i.e., without off-line post-processing of the data, which can often induce hard-to-control leakage of information from the testing data into the training data). The maximum TSS values for our method of predicting X-class flares is \(0.64\), comfortably surpassing the coveted state-of-the-art threshold, should one contend with alarm rate of \(40\%\). While we do not claim to have achieved close to optimal forecasts, this brief analysis demonstrates the great potential of the optimal homogeneous prediction methodology in the solar flare forecasting challenge.

References↩︎