Estimating the parameters of max-stable parametric models poses significant challenges, particularly when some parameters lie on the boundary of the parameter space. This situation arises when a subset of variables exhibits extreme values simultaneously, while the remaining variables do not—a phenomenon commonly referred to as an extreme direction. A novel estimator is proposed for the parameters of a general parametric mixture model, incorporating a threshold exceedances approach based on a pseudo-norm penalization. The latter plays a crucial role in accurately identifying parameters at the boundary of the parameter space. Additionally, the estimator comes with a data-driven algorithm to detect groups of variables corresponding to extreme directions. The performance of the estimator is assessed in terms of both parameter estimation and the identification of extreme directions through extensive simulation studies. Finally, the method is applied to two real-world datasets: discharge measurements at stations along the Danube river, and financial portfolio losses from stocks listed on the NYSE, AMEX, and NASDAQ. In both applications, the sets of variables that can become large simultaneously are identified.

Multivariate generalized Pareto distributions arise as limits of threshold exceedances and form a central model class for multivariate extremes. Existing inference methods based on the extremal variogram condition on the value of a single component, which can be statistically suboptimal. We generalize this approach by conditioning the multivariate generalized Pareto random vector Y to lie on arbitrary half-spaces. Specifically, for a direction vector v, we introduce the random vector Y^v = (Y | v^⊤ Y > 0) and define the associated v-variogram Γ_ij^v=Var(Y_i^v-Y_j^v). We establish the decomposition Y^v d= W^v+E1 into the so-called v-extremal function W^v and an independent exponential random variable E, and derive several results relating these random variables to each other. For logistic, Dirichlet, and Hüsler-Reiss multivariate generalized Pareto models, we derive closed-form expressions for Γ^v. In the Hüsler-Reiss case, we further derive new density representations and identify a distinguished resistance-curvature vector v_0 that uniquely centers the Gaussian law of W^v_0 while characterizing the least-mass half-space. On the statistical side, we introduce empirical v-variograms and show in a simulation study that the choice of v induces a pronounced bias-variance trade-off that is strongly related to the mass of the conditioning half-space. Moreover, combining information across multiple directions v can substantially reduce estimation variance relative to methods based on a single vector.

Economically responsible mitigation of multivariate extreme risks-such as extreme rainfall over large areas, large simultaneous variations in many stock prices, or widespread breakdowns in transportation systems-requires assessing the resilience of the systems under plausible stress scenarios. This paper uses Extreme Value Theory (EVT) to develop a new approach to simulating such multivariate extreme events. Specifically, we assume that after transformation to a standard scale the distribution of the random phenomenon of interest is multivariate regular varying and use this to provide a sampling procedure for extremes on the original scale. Our procedure combines a Wasserstein-Aitchison Generative Adversarial Network (WA-GAN) to simulate the tail dependence structure on the standard scale with joint modeling of the univariate marginal tails on the original scale. The WA-GAN procedure relies on the angular measure-encoding the distribution on the unit simplex of the angles of extreme observations-after transformation to Aitchison coordinates, which allows the Wasserstein-GAN algorithm to be run in a linear space. Our method is applied both to simulated data under various tail dependence scenarios and to a financial data set from the Kenneth French Data Library. The proposed algorithm demonstrates strong performance compared to existing alternatives in the literature, both in capturing tail dependence structures and in generating accurate new extreme observations.

In multivariate extreme value analysis, the tail dependence between some of the risk variables at hand may be weak, even when other variables do tend to become large simultaneously. Weak tail dependence may induce a substantial bias in estimation procedures based on the limiting multivariate (generalized) Pareto distribution of excesses over high thresholds. We consider a Hüsler–Reiss multivariate generalized Pareto model and, motivated by this issue, propose first- and second-order moment estimators of its variogram matrix constructed from a lower-tail-clipped version of the underlying random vector. The asymptotic normality of the proposed estimators is established. We demonstrate by simulation studies that they have lower bias than the empirical variogram estimator in certain cases, particularly when the dependence between components is weak. The estimators are applied to flood discharge data from the Danube river basin and the US flight delay data, showing that the tail dependence structure implied by the fitted model based on the first-order clipped moment estimator aligns more closely with the empirical tail dependence of the data than that based on the empirical variogram estimator.

We develop a unified statistical framework for attributing heatwaves as spatio-temporal phenomena under climate change. We quantify the impact of anthropogenic forcing on the probability and persistence of heatwaves not captured by standard marginal extreme-value approaches. Our methodology constructs a generative model for daily temperature fields that separates marginal nonstationarity from spatio-temporal dependence. We combine three components: a Bayesian spatial quantile regression model for the bulk of the data; a nonstationary spatial generalized extreme value model for tail behavior; and a copula-based model capturing both asymptotic dependence and independence in the extremes. The framework is applied to the CMIP6 MRI-ESM2 climate model, contrasting factual and counterfactual scenarios for probabilistic attribution. Our results show that the approach captures key heatwave characteristics inaccessible to traditional methods, enabling direct estimation of event-level attribution metrics. Overall, it provides a flexible basis for analyzing and attributing complex climate extremes as space-time objects.

We study cyclically monotone transport plans between measures in M_0(ℝ^d), the class of Borel measures on ℝ^d ∖{0} that are finite on sets bounded away from the origin but may have infinite total mass. We avoid moment assumptions and allow the transport cost to be infinite. This framework naturally arises for exponent measures in multivariate regular variation and includes other examples such as Lévy measures. We introduce the notion of a zero-coupling and establish existence of cyclically monotone zero-couplings for arbitrary pairs of measures in M_0(ℝ^d). Under a Hausdorff-dimension condition on the first measure and when at least one of the two measures has infinite mass, we prove uniqueness of the cyclically monotone zero-coupling, yielding an analogue of the Brenier–McCann theorem in this infinite-measure setting. We further derive a representation of such couplings through gradients of closed convex functions and identify conditions under which the zero-coupling is proper in the sense that the second measure is equal to the restriction to the punctured space of the push-forward of the first measure by a cyclically monotone transport map. Finally, we apply these results to regularly varying probability measures. We show that a cyclically monotone coupling between two such distributions admits a tail limit that coincides with the unique proper cyclically monotone zero-coupling between the corresponding exponent measures.

Probabilistic forecasts comprehensively describe the uncertainty in the unknown future outcome, making them essential for decision making and risk management. While several methods have been introduced to evaluate probabilistic forecasts, existing evaluation techniques are ill-suited to the evaluation of tail properties of such forecasts. However, these tail properties are often of particular interest to forecast users due to the severe impacts caused by extreme outcomes. In this work, we introduce a general notion of tail calibration for probabilistic forecasts, which allows forecasters to assess the reliability of their predictions for extreme outcomes. We study the relationships between tail calibration and standard notions of forecast calibration, and discuss connections to peaks-over-threshold models in extreme value theory. Diagnostic tools are introduced and applied in a case study on European precipitation forecasts

Estimating the parameters of max-stable parametric models poses significant challenges, particularly when some parameters lie on the boundary of the parameter space. This situation arises when a subset of variables exhibits extreme values simultaneously, while the remaining variables do not – a phenomenon referred to as an extreme direction in the literature. In this paper, we propose a novel estimator for the parameters of a general parametric mixture model, incorporating a penalization approach based on a pseudo-norm. This penalization plays a crucial role in accurately identifying parameters at the boundary of the parameter space. Additionally, our estimator comes with a data-driven algorithm to detect groups of variables corresponding to extreme directions. We assess the performance of our estimator in terms of both parameter estimation and the identification of extreme directions through extensive simulation studies. Finally, we apply our methods to data on river discharges and financial portfolio losses.

This paper explores strong and weak consistency of M-estimators for non-identically distributed data, extending prior work. Emphasis is given to scenarios where data is viewed as a triangular array, which encompasses distributional regression models with non-random covariates. Primitive conditions are established for specific applications, such as estimation based on minimizing empirical proper scoring rules or conditional maximum likelihood. A key motivation is addressing challenges in extreme value statistics, where parameter-dependent supports can cause criterion functions to attain the value $-\infty$, hindering the application of existing theorems.

A novel linear integration rule called $\textit{control neighbors}$ is proposed in which nearest neighbor estimates act as control variates to speed up the convergence rate of the Monte Carlo procedure on metric spaces. The main result is the $\mathcal{O}(n^{-1/2} n^{-s/d})$ convergence rate -- where $n$ stands for the number of evaluations of the integrand and $d$ for the dimension of the domain -- of this estimate for H\"older functions with regularity $s \in (0,1]$, a rate which, in some sense, is optimal. Several numerical experiments validate the complexity bound and highlight the good performance of the proposed estimator.

Motivated by examples from extreme value theory, but without using the theory of regularly varying time series or any assumptions about the marginal distribution, we introduce the general notion of a cluster process as a limiting point process of returns of a certain event in a time series. We explore general invariance properties of cluster processes that are implied by stationarity of the underlying time series. Of particular interest in applications are the cluster size distributions, and we derive general properties and interconnections between the size of an inspected and a typical cluster. While the extremal index commonly used in extreme value theory is often interpreted as the inverse of a “mean cluster size”, we point out that this only holds true for the expected value of the typical cluster size, caused by an effect very similar to the inspection paradox in renewal theory.

Regular vine sequences permit the organisation of variables in a random vector along a sequence of trees. Regular vine models have become greatly popular in dependence modelling as a way to combine arbitrary bivariate copulas into higher-dimensional ones, offering flexibility, parsimony, and tractability. In this project, we use regular vine structures to decompose and construct the exponent measure density of a multivariate extreme value distribution, or, equivalently, the tail copula density. Although these densities pose theoretical challenges due to their infinite mass, their homogeneity property offers simplifications. The theory sheds new light on existing parametric families and facilitates the construction of new ones, called X-vines. Computations proceed via recursive formulas in terms of bivariate model components. We develop simulation algorithms for X-vine multivariate Pareto distributions as well as methods for parameter estimation and model selection on the basis of threshold exceedances. The methods are illustrated by Monte Carlo experiments and a case study on US flight delay data.

The severity of multivariate extreme events is driven by the dependence between the largest marginal observations. The H\"usler-Reiss distribution is a versatile model for this extremal dependence, and it is usually parameterized by a variogram matrix. In order to represent conditional independence relations and obtain sparse parameterizations, we introduce the novel H\"usler-Reiss precision matrix. Similarly to the Gaussian case, this matrix appears naturally in density representations of the H\"usler-Reiss Pareto distribution and encodes the extremal graphical structure through its zero pattern. For a given, arbitrary graph we prove the existence and uniqueness of the completion of a partially specified H\"usler-Reiss variogram matrix so that its precision matrix has zeros on non-edges in the graph. Using suitable estimators for the parameters on the edges, our theory provides the first consistent estimator of graph structured H\"usler-Reiss distributions. If the graph is unknown, our method can be combined with recent structure learning algorithms to jointly infer the graph and the corresponding parameter matrix. Based on our methodology, we propose new tools for statistical inference of sparse H\"usler-Reiss models and illustrate them on large flight delay data in the U.S., as well as Danube river flow data.

Multivariate extreme value distributions are a common choice for modeling multivariate extremes. In high dimensions, however, the construction of flexible and parsimonious models is challenging. We propose to combine bivariate max-stable distributions into a Markov random field with respect to a tree. Although in general not max-stable itself, this Markov tree is attracted by a multivariate max-stable distribution. The latter serves as a tree-based approximation to an unknown max-stable distribution with the given bivariate distributions as margins. Given data, we learn an appropriate tree structure by Prim's algorithm with estimated pairwise upper tail dependence coefficients as edge weights. The distributions of pairs of connected variables can be fitted in various ways. The resulting tree-structured max-stable distribution allows for inference on rare event probabilities, as illustrated on river discharge data from the upper Danube basin.

The angular measure on the unit sphere characterizes the first-order dependence structure of the components of a random vector in extreme regions and is defined in terms of standardized margins. Its statistical recovery is an important step in learning problems involving observations far away from the center. In this paper, we test the goodness-of-fit of a given parametric model to the extremal dependence structure of a bivariate random sample. The proposed test statistic consists of a weighted $L_1$-Wasserstein distance between a nonparametric, rank-based estimator of the true angular measure obtained by maximizing a Euclidean likelihood on the one hand, and a parametric estimator of the angular measure on the other hand. The asymptotic distribution of the test statistic under the null hypothesis is derived and is used to obtain critical values for the proposed testing procedure via a parametric bootstrap. Consistency of the bootstrap algorithm is proved. A simulation study illustrates the finite-sample performance of the test for the logistic and H\"usler--Reiss models. We apply the method to test for the H\"usler--Reiss model in the context of river discharge data.

The Sliced-Wasserstein (SW) distance between probability measures is defined as the average of the Wasserstein distances resulting for the associated one-dimensional projections. As a consequence, the SW distance can be written as an integral with respect to the uniform measure on the sphere and the Monte Carlo framework can be employed for calculating the SW distance. Spherical harmonics are polynomials on the sphere that form an orthonormal basis of the set of square-integrable functions on the sphere. Putting these two facts together, a new Monte Carlo method, hereby referred to as Spherical Harmonics Control Variates (SHCV), is proposed for approximating the SW distance using spherical harmonics as control variates. The resulting approach is shown to have good theoretical properties, e.g., a no-error property for Gaussian measures under a certain form of linear dependency between the variables. Moreover, an improved rate of convergence, compared to Monte Carlo, is established for general measures. The convergence analysis relies on the Lipschitz property associated to the SW integrand. Several numerical experiments demonstrate the superior performance of SHCV against state-of-the-art methods for SW distance computation.

Graphical models with heavy-tailed factors can be used to model extremal dependence or causality between extreme events. In a Bayesian network, variables are recursively defined in terms of their parents according to a directed acyclic graph (DAG). We focus on max-linear graphical models with respect to a special type of graphs, which we call a tree of transitive tournaments. The latter are block graphs combining in a tree-like structure a finite number of transitive tournaments, each of which is a DAG in which every two nodes are connected. We study the limit of the joint tails of the max-linear model conditionally on the event that a given variable exceeds a high threshold. Under a suitable condition, the limiting distribution involves the factorization into independent increments along the shortest trail between two variables, thereby imitating the behavior of a Markov random field. We are also interested in the identifiability of the model parameters in case some variables are latent and only a subvector is observed. It turns out that the parameters are identifiable under a criterion on the nodes carrying the latent variables which is easy and quick to check.

When modeling a vector of risk variables, extreme scenarios are often of special interest. The peaks-over-thresholds method hinges on the notion that, asymptotically, the excesses over a vector of high thresholds follow a multivariate generalized Pareto distribution. However, existing literature has primarily concentrated on the setting when all risk variables are always large simultaneously. In reality, this assumption is often not met, especially in high dimensions. In response to this limitation, we study scenarios where distinct groups of risk variables may exhibit joint extremes while others do not. These discernible groups are derived from the angular measure inherent in the corresponding max-stable distribution, whence the term extreme direction. We explore such extreme directions within the framework of multivariate generalized Pareto distributions, with a focus on their probability density functions in relation to an appropriate dominating measure. Furthermore, we provide a stochastic construction that allows any prespecified set of risk groups to constitute the distribution’s extreme directions. This construction takes the form of a smoothed max-linear model and accommodates the full spectrum of conceivable max-stable dependence structures. Additionally, we introduce a generic simulation algorithm tailored for multivariate generalized Pareto distributions, offering specific implementations for extensions of the logistic and Hüsler-Reiss families capable of carrying arbitrary extreme directions.

When passing from the univariate to the multivariate setting, modelling extremes becomes much more intricate. In this introductory exposition, classical multivariate extreme value theory is presented from the point of view of multivariate excesses over high thresholds as modelled by the family of multivariate generalized Pareto distributions. The formulation in terms of failure sets in the sample space intersecting the sample cloud leads to the over-arching perspective of point processes. Max-stable or generalized extreme value distributions are finally obtained as limits of vectors of componentwise maxima by considering the event that a certain region of the sample space does not contain any observation.

The angular measure on the unit sphere characterizes the first-order dependence structure of the components of a random vector in extreme regions and is defined in terms of standardized margins. Its statistical recovery is an important step in learning problems involving observations far away from the center. In the common situation that the components of the vector have different distributions, the rank transformation offers a convenient and robust way of standardizing data in order to build an empirical version of the angular measure based on the most extreme observations. We provide a functional asymptotic expansion for the empirical angular measure in the bivariate case based on the theory of weak convergence in the space of bounded functions. From the expansion, not only can the known asymptotic distribution of the empirical angular measure be recovered, it also enables to find expansions and weak limits for other statistics based on the associated empirical process or its quantile version.