Bayesian Non-Parametric Inference for Multivariate Peaks-over-Threshold Models

We consider a constructive definition of the multivariate Pareto that factorizes the random vector into a radial component and an independent angular component. The former follows a univariate Pareto distribution, and the latter is defined on the surface of the positive orthant of the infinity norm unit hypercube. We propose a method for inferring the distribution of the angular component by identifying its support as the limit of the positive orthant of the unit p-norm spheres and introduce a projected gamma family of distributions defined through the normalization of a vector of independent random gammas to the space. This serves to construct a flexible family of distributions obtained as a Dirichlet process mixture of projected gammas. For model assessment, we discuss scoring methods appropriate to distributions on the unit hypercube. In particular, working with the energy score criterion, we develop a kernel metric that produces a proper scoring rule and presents a simulation study to compare different modeling choices using the proposed metric. Using our approach, we describe the dependence structure of extreme values in the integrated vapor transport (IVT), data describing the flow of atmospheric moisture along the coast of California. We find clear but heterogeneous geographical dependence.


Introduction
The statistical analysis of extreme values focuses on inferences for rare events that correspond to the tails of probability distributions.As such, it is a key ingredient in the risk assessment of phenomena that can have strong societal impacts like floods, heat waves, high concentration of pollutants, crashes in the financial markets, among others.The fundamental challenge of extreme value theory (EVT) is to use information, collected over limited periods of time, to extrapolate to long time horizons.This sets EVT apart from most of statistical inference, where the focus is on the bulk of the distribution.Extrapolation to the tails of the distributions is possible thanks to theoretical results that give asymptotic descriptions of the probability distributions of extreme values.
Inferential methods for the extreme values of univariate observations are well established, and software is widely available; see, for example, [1].For variables in one dimension, the application of EVT methods considers the asymptotic distribution of either the maxima calculated for regular blocks of data, or the values that exceed a certain threshold.The former leads to a generalized extreme value (GEV) distribution that depends on three parameters.The latter leads to a generalized Pareto (GP) distribution, which depends on a shape and a scale parameter.Likelihood-based approaches to inference can be readily implemented in both cases.In the multivariate case, the GEV theory is well developed; see, for example, [2], but the inferential problem is complicated by the fact that there is no parametric representation of the GEV.This problem is inherited by the peaks over threshold (PoT) approach and compounded by the fact that there is no unique definition of an exceedance of a multivariate threshold, as there is an obvious dependence on the norm that is used to measure the size of a vector.
During the last decade or so, much work has been done in the exploration of the definition and properties of an appropriate generalization of the univariate GP distribution to a multivariate setting.To mention some of the papers on the topic, the work of [3] defines the generalized Pareto distribution, with further analysis of these classes of distributions presented in [4,5].A recent review of the state of the art in multivariate peaks over threshold modeling using generalized Pareto is provided in [6] while [7] provides insight on the theoretical properties of possible parametrizations.These are used in [8] for likelihoodbased models for PoT estimation.A frequently used method for describing dependence in multivariate distributions is to use a copula.Refs.[9,10] provide successful examples of this approach in an EVT framework.
Ref. [11] presents a constructive definition of the Pareto process, which generalizes the GP to an infinite-dimensional setting.It consists of decomposing the process into independent radial and angular components.Such an approach can be used in the finitedimensional case, where the angular component contains information pertaining to the dependence structure of the random vector.Based on this definition, we present a novel approach for modeling the angular component with families of distributions that provide flexibility and can be applied in a moderately large dimensional setting.Our focus on the angular measure is similar to that in [12][13][14], which consider Bayesian non-parametric approaches.Yet, our approach differs in that it is established in the peaks-over-threshold regime and uses a constructive definition of the multivariate GP based on the infinity norm.The approach proposed in this paper adds to the growing literature on Bayesian models for multivariate extreme value analysis (see, for example, [12][13][14][15]), providing a model that has strong computational advantages due its structural simplicitly, achieves flexibility using a mixture model, and scales well to moderately large dimensions.
The remainder of this paper is outlined as follows.Section 2 comprises a brief review of multivariate PoT, detailing the separation of the radial measure from the angular measure.Section 3 details our approach for estimating the angular measure, based on transforming an arbitrary distribution supported in R d + onto unit hyper-spheres defined by p-norms.Section 4 develops criteria for model selection in the support of the angular measure.Section 5 explores the efficacy of the proposed approach on a set of simulated data, and, acknowledging the relevance of extreme value theory to climatological events [16][17][18], estimates the extremal dependence structure for a measure of water vapor flow in the atmosphere, used for identifying atmospheric rivers.Finally, Section 6 presents our conclusions and discussion.
Throughout the paper, we adopt the operators ∧ to denote minima and the ∨ to denote maxima.Thus ∧ i s i = min i s i , and ∨ i s i = max i s i .These operators can also be applied component-wise between vectors such as a ∧ b = (a 1 ∧ b 1 , a 2 ∧ b 2 , . ..).Similarly, we apply inequality and arithmetic operators operators to vectors.For example, a ≤ b, and interpret them component-wise.We use uppercase to indicate random variables, lowercase to indicate points, and bold-face to indicate vectors or matrices thereof.

A Multivariate PoT Model
To develop a multivariate PoT model for extreme values, consider a d-dimensional random vector W = (W 1 , . . ., W d ) with cumulative distribution F. A common assumption on W is that it is in the so-called domain of attraction of a multivariate max-stable distribution, G. Thus, following [7], assume that there exists sequences of vectors a n and b n , such that lim n→∞ F n (a n w + b n ) = G(w).G is a d-variate generalized extreme value distribution.Notice that, even though the univariate marginals are obtained from a three-parameter family, there is no parametric form to represent G. Taking logarithms and expanding, we have that lim where a −1 n indicates element-wise inversion, and {W ̸ ≤ b n } denotes the set where at least one coordinate is above the corresponding component of b n .H is a multivariate Pareto distribution.It corresponds to a joint distribution conditional on exceeding a multivariate threshold.H is defined by a non-parametric function governing the multivariate dependence and two d-dimensional vectors of parameters that control the shapes and scales of the marginals.We denote these as ξ for the shapes and σ for the scales.Ref. [7] provides a number of stochastic representations for H.In this paper, we focus on a particular one that is proposed in [11].To this end, we denote as Z a random variable with distribution H where ξ = 1 and σ = 0.Then, Z = RV where R and V are independent.R = ∥Z∥ ∞ = ∨ d i=1 Z i is distributed as a standard Pareto random variable, and ∞ , the positive orthant of the unit sphere under L ∞ norm, with distribution Φ.This representation is central to the methods proposed in this paper.R and V are referred to, respectively, as the radial and angular components of H.The angular measure controls the dependence structure of Z in the tails.In view of this, to obtain a PoT model, we seek a flexible model for the distribution of V ∈ S d−1 ∞ , based on a Bayesian non-parametric model.The approach considered in [6] focuses on the limiting conditional distribution H.An alternative approach to obtaining a limiting PoT distribution consists of assuming that regular variation (see, for example, [19]) holds for the limiting distribution of W, implying that for some measure µ that is referred to as the exponent measure.µ has the homogeneity property µ(tA) = t −1 µ(A).Letting ρ = ∥W∥ p , p > 0 and θ = W /ρ, define Ψ(B) = µ({w : ρ > 1, θ ∈ B}), which is referred to as the angular measure.After some manipulations, we obtain that Thus, a model for the exponent measure induces a model for the limiting distribution conditional on the observations being above a threshold defined with respect to their p-norm.The constraint that all marginals of µ correspond to a standard Pareto distribution leads to the so-called moment constraints on Ψ, consisting of Inference for the limiting distribution of the exceedances needs to account for the normalizing constant in Equation (1) as well as the moment constraints.Because of these issues, in this paper, we prefer to follow the limiting conditional distribution approach.An example of the application of the regular variation approach using p = 1 is developed in [13].

Estimation of the Angular Measure
To infer the PoT distribution, we consider two steps: First we estimate the shape and scale parameters for the multivariate Pareto distribution, using the univariate marginals; then we focus on the dependence structure in extreme regions by proposing a flexible model for the distribution of V .Consider w i , i = 1, . . ., n a collection of realizations of W. We start by setting a large threshold b t,ℓ for the ℓ-th marginal, ℓ = 1, . . ., d.Then, the distribution of W ℓ , conditional on exceeding the threshold, can be approximated with a generalized univariate Pareto.Thus, where (•) + indicates the positive part function.We set b t,l = F−1 ℓ (1 − 1/t), the empirical (1 − 1/t)-quantile.We then estimate ξ ℓ and σ ℓ , for each ℓ, using likelihood-based methods.To estimate the angular distribution, we standardize each of the marginals.The standardization yields Note that z iℓ > 1 implies that w iℓ > b t,ℓ , meaning that the observation w i is extreme in the ℓ-th dimension.Thus, r i = ∥z i ∥ ∞ > 1 implies that at least one dimension has an extreme observation and corresponds to a very extreme observation when t is large.We focus on the observations that are such that r i > 1.These provide a sub-sample of the standardized original sample.We define ∞ .These vectors are used for the estimation of Φ.

Projected Gamma Family
At the core of our PoT method is the development of a distribution on S d−1 p = {s : s ∈ R d + , ∥s∥ p = 1}, where, for p > 0, ∥ • ∥ p is the L p -norm of a vector x ∈ R d , defined as The absolute and Euclidean norms are obtained for p = 1 and p = 2 respectively, and the L ∞ norm can be obtained as a limit: To obtain a distribution on S d−1 p , we start with a vector in x ∈ R d + and normalize it to obtain y = x/∥x∥ p ∈ S d−1 p . Figure 1a shows the progression of S 1 p asymptotically towards S 1 ∞ as p increases; Figure 1b shows the relative positions of data points normalized to S 1 p for selected p.A natural distribution to consider in R d + is given by a product of independent univariate Gamma distributions.Let X ∼ ∏ d ℓ=1 Ga(X ℓ | α ℓ , β ℓ ).α ℓ and β ℓ are the shape and scale parameters, respectively.For any finite p > 0, letting y d = (1 − ∑ d−1 ℓ=1 y p ℓ ) 1/p , the transformation is invertible with The Jacobian of the transformation takes the form The normalization provided by T maps a vector in R d + onto S d−1 p .With a slight abuse of terminology, we refer to it as a projection.Using Equations (3)-( 5), we have the joint density Integrating out r yields the resulting Projected Gamma density defined for y ∈ S d−1 p and for any finite p > 0. To avoid identifiability problems when estimating the shape and scale parameters, we set β 1 = 1.Ref. [20] obtain the density in Equation ( 7) for p = 2 as a multivariate distribution for directional data using spherical coordinates.For y ∈ S d−1 1 and β ℓ = β for all ℓ, the density in Equation ( 7) corresponds to that of a Dirichlet distribution.
The projected gamma family is simple to specify and has very tractable computational properties.Thus, we use it as a building block for the angular measure Φ models.To build a flexible family of distributions in S d−1 p , we consider mixtures of projected gamma densities defined as where θ = (α, β).Following a Bayesian non-parametric approach [21][22][23], we assume that G is drawn from a random measure.In particular, assuming a Dirichlet process prior for G, we have a hierarchical formulation of the mixture model that, for a vector of observations y i , is given by where DP denotes a Dirichlet process, η is the precision parameter, and G 0 is the centering distribution.
Unfortunately, in the limit when p → ∞, the normalizing transformation is not differentiable.Thus, a closed-form expression like Equation (7) for the projected gamma density on S d−1 ∞ is not available.Instead, we observe that for a sufficiently large p, S d−1 p will approach S d−1 ∞ .With that in mind, our strategy consists of describing the angular distribution Φ using a sample-based approach with the following steps: (i) Apply the transformation in Equation ( 2) to the original data; (ii) Obtain the subsample of the standardized observations that satisfy R > 1; (iii) Take a finite p and project the observations onto S d−1 p ; (iv) Fit the model in Equation (8) to the resulting data and obtain samples from the fitted model; (v) project the resulting samples onto S d−1 ∞ .For step (iv), we use a Bayesian approach that is implemented using a purposely developed Markov chain Monte Carlo sampler described in the next section.

Tail Probabilities for the PoT Model
A measure that is used to characterize the strength of the dependence in the tail for two random variables, Z 1 and Z 2 , with marginal distributions F 1 and F 2 is given by [1] χ 12 provides information about the distribution of extremes for the variable Z 1 given that Z 2 is very large.When χ 12 > 0, Z 1 and Z 2 are said to be asymptotically dependent; otherwise, they are asymptotically independent.The following result provides the asymptotic dependence coefficient between two components of Z for our proposed PoT model.
Proof.Denote as F ℓ the marginal distribution of Z ℓ .To obtain χ ȷℓ , we need Using the fact that V ℓ > 0, ∀ℓ almost surely, we have that the former is equal to where the second identity is justified by the fact that V i is bounded and Equation (10) implies that χ ȷℓ > 0, and so, no asymptotic independence is possible under our proposed model.For the analysis of extreme values, it is of interest to calculate the multivariate conditional survival function.The following result provides the relevant expression as a function of the angular measure.Proposition 2. Assume the same conditions as Proposition 1.Let α ⊂ {1, . . ., d} be a collections of indexes.Then The proof uses a similar approach to the proof of Proposition 1. Equations ( 10) and (11) provide relevant tools for inference on the tail behavior of the joint distribution of the observations.The expressions can be readily calculated within a sample-based inferential approach like the one considered in the following section.

Inference for the Projected Gamma Mixture Model
To perform inference for our proposed PoT model, we develop an iterative samplebased approach.We implement a Markov chain Monte Carlo method that, for a given iteration, groups observations into stochastically assigned clusters, where members of a cluster share distributional parameters [23,24].Building out the methods of inference for Equation (9), let n be the number of observations in cluster j, not including observation i.Let J (−i) be the number of extant clusters, not including any singleton containing obser-vation i.Under this model, the probability of cluster membership for a given observation is proportional to where the top case is iterating over extant clusters j = 1, . . ., J (−i) , and the bottom case is for a new cluster.If G 0 is not a conjugate prior for the kernel density, the integral in the above formula may not be available in closed form.We sidestep this using Algorithm 8 from [25]: by Monte Carlo integration, we draw m candidate clusters, α j , β j for j = J (−i) + 1, . . ., J (−i) + m from G 0 .Then, we sample the cluster indicator γ i from extant or candidate clusters, where the probability of cluster membership is proportional to Again, the top case is iterating over extant clusters, and now the bottom case is iterating over new candidate clusters.If a candidate cluster is selected, then γ i = J = J (−i) + 1, and the associated cluster parameters are saved.
A key feature of the the projected Gamma distribution is its computational properties.We augment P G(y i | α i , β i ) by introducing a latent radial component r i , for each observation.Using Equation ( 6) we observe that the full conditional of r i is easy to sample from, as it is given as Moreover, the full conditional for α j , β j is then proportional to Note that the ordering of the products can be reversed in Equation ( 14), indicating that with appropriate choice of centering distribution, the full conditionals for α j , β j can become separable by dimension, and thus inference on α jℓ , β jℓ can be done in parallel for all j, ℓ.
We first consider a centering distribution given by a product of independent Gammas: This model is completed with independent Gamma priors on ξ ℓ , τ ℓ , ζ ℓ , and σ ℓ .We also assume a Gamma prior on η, which is updated via the procedure outlined in [26].We refer to this model as the projected gamma-gamma (PG-G) model.An advantage of the PG-G model is that, thanks to conjugacy, the rate parameters β jℓ can easily be integrated out for inference on α j .Then, the full conditional for α jℓ takes the form for ℓ = 2, . . ., d.For ℓ = 1, as β 1 := 1, the full conditional takes the simpler form Samples of α jℓ can thus be obtained using a Metropolis step.In our analysis, we first transform α jℓ to the log scale and use a normal proposal density.The full conditional for β is for ℓ = 2, . . ., d. Updating β jℓ is done via a Gibbs step.The hyper-parameters ξ ℓ , τ ℓ , ζ ℓ , σ ℓ follow similar gamma-gamma update relationships.We also explore a restricted form of this model, where β ℓ = 1 for all ℓ.Under this model, we use the full conditional in Equation ( 17) for all ℓ and omit inference on ζ, σ.We refer to this model as the projected restricted gamma-gamma (PRG-G) model.
The second form of centering distribution we explore is a multivariate log-normal distribution on the shape parameters α j with independent gamma β jℓ rate parameters.
This model is completed with a normal prior on µ, an inverse Wishart prior on Σ, and Gamma priors on ζ ℓ , σ ℓ , and η.This model is denoted as the projected gamma-log-normal (PG-LN) model.We also explore a restricted Gamma form of this model as above, where β ℓ = 1 for all ℓ.This is denoted as the projected restricted gamma-log-normal (PRG-LN) model.
Updates for α can be accomplished using a joint Metropolis step, where β jℓ for ℓ = 2, . . ., d have been integrated out of the log-density.That is, The inferential forms for β jℓ and its priors are the same as for PG-G.The normal prior for µ is conjugate for the log-normal α j and can be sampled via a Gibbs step.Finally, the inverse Wishart prior for Σ is again conjugate to the log-normal α j , implying that it can also be sampled via a Gibbs step.
To effectively explore the sample space with a joint Metropolis step, as well as to speed convergence, we implement a parallel tempering algorithm [27] for the log-normal models.This technique runs parallel MCMC chains at ascending temperatures.That is, for chain s, the posterior density is exponentiated by the reciprocal of temperature t s .For the cold chain, t 1 := 1.Let E s be the log-posterior density under the current parameter state for chain s.Then states between chains r, s are exchanged via a Metropolis step with probability Higher temperatures serve to flatten the posterior distribution, meaning hotter chains have a higher probability of making a given transition or will make larger transitions.As such, they will more quickly explore the parameter space and share information gained through state exchange.

Scoring Criteria for Distributions on the Infinity-Norm Sphere
In order to assess and compare the estimation of a distribution on S d−1 ∞ , we consider the theory of proper scoring rules developed in [28].As mentioned in Section 3.1, our approach does not provide a density on S d−1 ∞ , restricting our ability to construct model selection criteria to sample-based approaches.To this end, we employ the energy score criterion introduced therein.
The energy score criterion, defined for a general probability distribution P with a finite expectation, is developed as where g is a kernel function.The score defined in Equation ( 20) can be evaluated using samples from P with the help of the law of large numbers.Moreover, Theorem 4 in [28] states that if g(•, •) is a negative definite kernel, then S(P, x) is a proper scoring rule.Recall that a real valued function g is a negative definite kernel if it is symmetric in its arguments, and ∑ n i=1 ∑ n j=1 a i a j g(x i , x j ) ≤ 0 for all positive integers n and any collection a 1 , . . ., a n ∈ R such that ∑ n i=1 a i = 0.In a Euclidean space, these conditions are satisfied by the Euclidean distance [29].However, for observations on different faces of S d−1 ∞ , the Euclidean distance will underestimate the geodesic distance, the actual distance required to travel between the two points.Let ∞ , x ℓ = 1} comprise the ℓth face.For points on the same face, the Euclidean distance corresponds to the length of the shortest possible path in S d−1 ∞ .For points on different faces, the Euclidean distance is a lower bound for that length.
For a finite p, the shortest connecting path between two points in S d−1 p is the minimum geodesic; its length satisfying the definition of a distance.Thus its length can be used as a negative definite kernel for the purpose of defining an energy score.Unfortunately, as p → ∞, the resulting surface S d−1 ∞ is not differentiable, implying that routines to calculate geodesics are not readily available.However, as S d−1 ∞ is a portion of a d-cube, we can borrow a result from geometry [30] stating that the length of the shortest path between two points on a geometric figure corresponds to the length of a straight line drawn between the points on an appropriate unfolding, rotation, or net of the figure from a d-dimensional to a d − 1-dimensional space.The optimal net will have the shortest straight line between the points, as long as that line is fully contained within such a net.As S d−1 ∞ has d faces-each face pairwise adjacent, there are d! possible nets.However, we are only interested in nets that begin and end on the source and destination faces, respectively, reducing the number of nets under consideration to ∑ d−2 k=0 ( d−2 k ).This is still computationally burdensome for a large number of dimensions.However, we can efficiently establish an upper bound on the geodesic length.We use this upper bound on geodesic distance as the kernel function for the energy score.
To calculate the energy score, we define the kernel , for ℓ, ȷ ∈ {1, . . ., d}.For ℓ ̸ = ȷ, the transformation P ȷℓ (•) required to rotate the ȷth face along the ℓth axis produces the vector b ′ , where We then have that The last equality is due to the fact that a and b ′ belong to the same hyperplane.
Using the rotation in Proposition 3, we obtain the following result.
Proposition 4. g is a negative definite kernel.
Proof.For a given n, consider an arbitrary set of points a 1 , . . ., a n ∈ S d−1 ∞ , and α 1 , . . ., where a ′ ȷ is defined as in Proposition 3. The last equality holds as ∥x − x ′ ∥ 2 , x, x ′ ∈ R d is negative definite [28] Proposition 3 provides a computationally efficient way to evaluate the proper scoring rule S ES defined on S d−1 ∞ for each observation.For the purpose of model assessment and comparison, we report the average S ES taken across all observed data and notice that the smaller the score, the better.

Data Illustrations
We apply the aforementioned models to simulated angular data.We then consider the analysis of atmospheric data.To tackle the difficult problem of assessing the convergence an MCMC chain for a large-dimensional model, we monitor the log-posterior density.In all the examples considered, MCMC samples produced stable traces of the log-posterior in less than 40,000 iterations.We use that as a burn-in and thereafter sample 10,000 additional iterations.We then thin the chain by retaining one every ten samples, to obtain 1000 total samples.These are used to generate samples from the posterior predictive densities.We used two different strategies to implement the MCMC samplers.For the models whose DP prior is centered around a log-normal distribution, we used parallel tempering.This serves to overcome the very slow mixing that we observed in these cases.The temperature ladder was set as t s = 1.3 s for s ∈ {0, 1, . . ., 5}.This was set empirically in order to produce acceptable swap probabilities both for the simulated data and real data.Parallel tempering produces chains with good mixing properties but has a computational cost that grows linearly with the number of temperatures.Thus, for the gamma-centered models, we used a single chain.We leverage the fast speed of each iteration to obtain a large number of samples, which are appropriately thinned to deal with a mild autocorrelation.In summary, the strategy for log-normal centered models is based on a costly sampler with good mixing properties.The strategy for the gamma-centered models is based on a cheap sampler that can be run for a large number of iterations.
Our hyperprior parameters are set as follows: For the gamma-centered models (PG-G, PRG-G), the shape parameter for the centering distribution ξ ℓ ∼ G(1, 1), and rate parameter τ ℓ ∼ G(2, 2).For the log-normal centered models (PG-LN, PRG-LN), the centering distribution's log-mean µ ∼ N d (0, I d ), and covariance matrix Σ ∼ IW (d + 10, (d + 10)I d ).These values are intended such that draws from the prior for Σ will weakly tend towards the identity matrix.For models learning rate parameters β jℓ (PG-G, PG-LN), the centering distribution's shape parameter ζ ∼ G(1, 1) and rate parameter σ ∼ G(2, 2).The choice of the G(2, 2) for rate parameters places little mass near 0 in order to draw estimates for the value away from 0 for numerical stability.

Simulation Study
The challenging problem in multivariate EVT is to capture the dependence structure of the limiting distribution.To this end, we focus our simulation study specifically on the angular component.To evaluate our proposed approach for angular measure estimation, we consider simulated datasets on S d−1 ∞ for values of d ranging from 2 to 32.We generated each dataset as a mixture of multivariate log-normal distributions projected onto S d−1 ∞ .The generation procedure is detailed in Algorithm 1.We produced ten replicates of each configuration.We consider two gamma-centered and two log-normal-centered DP mixture models, with and without restrictions in each case.To perform a comparative analysis, we fitted the pairwise betas model proposed in [31].We chose this model for comparison, as it similarly works to capture a complex dependence structure on an S d−1 p sphere, albeit with p = 1, and is implemented in the readily available package BMAmevt in R [32], which can provide samples from the posterior predictive distribution.These samples are needed for the calculation of the energy scores that are at the basis of our comparison.In addition, BMAmevt can be fitted to moderately large multivariate observations.For the DP mixture models, the data are projected onto S d−1 10 .For the other two models, they are projected onto S d−1 1 .We sampled each model for 50,000 iterations, dropping the first 40,000 as burn-in and thinning to keep every 10th iteration after.These settings were intended to provide a consistent sampling strategy that would work with every model, even if inefficient for some.
Algorithm 1 Simulated Angular Dataset Generation Routine.µ j , Σ j are the parameters of the mixture component distribution; π is the probability vector assigning weight mixture components; δ i is the mixture component identifier for each simulated observation.
for n iter in [1, . .[2,4,8,16,24,32] and save.end for end for end for Figure 2 shows the average rise over baseline in energy score, as calculated on S d−1 ∞ using the kernel metric described in Proposition 3, for models trained on simulated data.After training a model, a posterior predictive dataset is generated, and the energy score is calculated as a Monte Carlo approximation of Equation (20).In our analysis, after burnin and thinning, we had 1000 replicates from the posterior distribution and generated 10 posterior predictive replicates per iteration.The baseline value is the energy score of a new dataset from the same generating distribution as the training dataset evaluated against the training dataset.For the simulated data, we observe that the projected gamma models dominate the other two options considered, regardless of the choice of centering distribution.The projected restricted gamma models with a multivariate log-normal centering distribution appear to be dominated by the models based on the alternative centering distributions.Moreover, the performance deteriorates with the increase in dimensionality.Additionally, models centered around the log-normal distribution incur the computational cost of multivariate normal evaluation and parallel tempering, taking approximately six times longer to sample relative to the gamma models.We also note that the computational cost of the pairwise betas model grows combinatorially, with a sample space of dimension ( d 2 ) + 1.By comparison, the sample space for PG-G and PRG-G are 2(J + 1)d and (J + 1)d, respectively, where J is the number of extant clusters, with much of that inference able to be done in parallel.In our testing, for low-dimensional problems, BMAmevt was substantially faster than any of our proposed DP mixture models.However, for examples with high numbers of dimensions, the computational time for BMAmevt was greater than that for PG-G.We compare computing times in our data analysis in Table 1b.) for various models fitted to simulated data, with ascending count of mixture components (indicated by plot heading) and number of dimensions (indicated by horizontal axis).Note that pairwise betas is a moment-restricted model.

Integrated Vapor Transport
The integrated vapor transport (IVT) is a two-component vector that tracks the flow of the total water volume in a column of air over a given area [33].IVT is increasingly used in the study of atmospheric rivers because of its direct relationship with orographically induced precipitation [34].Atmospheric rivers (AR) are elongated areas of high local concentration of water vapor in the atmosphere that transport water from the tropics around the world.AR can cause extreme precipitation, something that is usually associated with very large values of the IVT magnitude over a whole geographical area.In spite of this, AR are fundamental for the water supply of areas like California.Thus, the importance of understanding the extreme behavior of IVT includes extreme tail dependence.We consider datasets that correspond to IVT estimated at two different spatial resolutions.The coarse-resolution dataset is obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) Interim reanalysis (ERA-Interim) [35,36].The high-resolution dataset corresponds to the latest ECMWF observational product, ERA5 [37].
Our data correspond to daily average values for the IVT magnitude along the coast of California.The ERA-Interim data used cover the time period 1979 through 2014 (37 years), omitting leap days, and eight grid cells that correspond to the coast of California.The ERA5 data cover the time period 1979 through 2019 (42 years) with the same restriction and 47 grid cells for the coast of California.This gives us the opportunity to illustrate the performance of our method in multivariate settings of very different dimensions.Figure 3 provides a visual representation of the area these grid cells cover.Fitting our models to the IVT data requires some pre-processing.First, we subset the data to the rainy season, which in California runs roughly from November to March.Following the approach described in Section 3, we estimate the shape and scale parameters of a univariate GP in each dimension, using maximum likelihood.We set the threshold in each dimension ℓ as b t,ℓ = F−1 ℓ (1 − t −1 ), where F is the empirical CDF and t = 20, which corresponds to the 95 percentile.We then use the transformation in Equation ( 2) to standardize the observations.Dividing each standardized observation by its L ∞ norm, we obtain a projection onto S d−1 ∞ .As the data correspond to a daily time series, the observations are temporally correlated.For each group of consecutive standardized vectors z i such that ∥z i ∥ ∞ > 1, we retain only the vector with the largest L ∞ norm.The complete procedure is outlined in Algorithm 2.
After subsetting the ERA-Interim data to the rainy season, we have 5587 observations.After the processing and declustering described in Algorithm 2, this number reduces to 511 observations.A pairwise plot of the transformed data after processing and declustering is presented in Figure 4. From this, we note that the marginal densities display strong similarities, with a large spike near 0 and a small spike near 1.A value of 1 in a particular axis indicates that the standardized threshold exceedance was the largest in that dimension.The off-diagonal plots correspond to pairwise density plots.We observe that some site pairs, such as (1, 2), (7,8), and especially (4,5), have the bulk of their data concentrated in a small arc along the 45 • ; in other site combinations, such as (3, 6), (2, 7), or (1,8), the data is split, favoring one side or the other of the 45 • line.For the ERA5 data, after subsetting, we have 6342 observations, which reduces to 532 observations after processing and declustering.We fit the PG-G, PRG-G, PG-LN, and PRG-LN models to both datasets.Table 1a shows the values of the estimated energy scores for the different models considered.We observe that, contrary to the results in the simulation study in Figure 2, the preferred model is the projected restricted gamma models, though for the lower-dimensional ERA-Interim data, all models perform comparably.Table 1b shows the computing times needed to fit the different models to the two datasets.We see the effect of dimensionality on the various models; for gamma-centered models, it grows linearly; for the log-normal centered model, it will grow superlinearly as matrix inversion becomes the most costly operation.For BMAmevt, its parameter space grows combinatorially with the number of dimensions, and thus so does computational complexity and sampling time.
We consider an exploration of the pairwise extremal dependence using Monte Carlo estimates of the coefficients in Equation (10).For this, we use samples obtained from the PRG-G model.Figure 5 provides a graphical analysis of the results.The coefficients achieve values between 0.286 and 0.759 for the ERA-Interim data and between 0.181 and 0.840 for the ERA5 data.The greater range in dependence scores observed in the ERA5 data versus ERA-Interim speaks to the greater granularity of the ERA5 data, indicating that distance between locations is a strong contributor to the strength of the pairwise asymptotic dependence.The highest coefficients are 0.759 for locations 4 and 5 in the ERA-Interim data and 0.840 for locations 1 and 2 in the ERA5 data.Clearly, pairwise asymptotic dependence coefficients tell a limited story, as a particular dependence may include more than two locations.We can, however, glean some information from the patterns that emerge in two dimensions.For the ERA-Interim data, we observe a possible cluster between cells 5-8, indicating a strong dependence among these cells.Analogously, for the ERA5 data, we observe three possible groups of locations.Figure 6 shows, for the ERA-Interim data under the PRG-G model, the conditional survival curve defined in Equation (11) for one dimension conditioned on all other dimensions being greater than their (fitted) 90th percentile.Figure 7 presents the bi-variate conditional survival function, conditioning on all other dimensions.These results illustrate quantitatively how extremal dependence affects the shape of the conditional survival curves.The two top panels represent the joint survival function between grid locations 4 and 5, which are shown in Figure 5 to exhibit strong extremal dependence.We observe that the joint survival surface is strongly convex.The bottom panels represent the joint survival surface between grid locations 1 and 5, which exhibited low extremal dependence.In this case, the shape of the contours tends to be concave, quite different from the shapes observed in the top panels.Using our proposed scoring criteria, we explored the effect of the choice of p on the final results.Using the simulated data, generated from a mixture of projected gammas, we were unable to observe sizeable differences in the scores for p ranging between 1 and 15.However, for the IVT data, we observed a drop in the energy score associated with higher p, with a diminishing effect as p increased.We observed no significant differences in the performance of the model that uses p = 10, which corresponds to the analysis presented, relative to the one that uses p = 15.

Conclusions
In this paper, we have built upon the definition of the multivariate Pareto described in [11] to establish a useful representation of its dependence structure through the distribution of its angular component, which is supported on the positive orthant of the unit hypersphere under the L ∞ norm, S d−1 ∞ .Due to the inherent difficulty of obtaining the likelihood of distributions with support on S d−1 ∞ , our method transforms the data to S d−1 p , fits then using mixtures of products of independent gammas, then transforms the predictions back to S d−1 ∞ .As S d−1 p converges to S d−1 ∞ as p → ∞, we expect the proposed resampling to be efficient for large enough p.In fact, our exploration of the simulated and real data indicates that the procedure is robust to the choice of moderately large values of p.Our method includes two inferential steps.The first consists of the estimation of the marginal Pareto distributions; the second consists of the estimation of the angular density.Parameter uncertainty incurred in the former is not propagated to the latter.Conceptually, an integrated approach that accounts for all the estimation uncertainty is conceivable.Unfortunately, this leads to posterior distributions with complex data-dependent restrictions that are very difficult to explore, especially in large dimensional settings.In fact, our attempts to

Figure 2 .
Figure 2. Average energy score rise over baseline (on S d−1∞ ) for various models fitted to simulated data, with ascending count of mixture components (indicated by plot heading) and number of dimensions (indicated by horizontal axis).Note that pairwise betas is a moment-restricted model.

Algorithm 2 1 ℓ 1
Data preprocessing to isolate and transform data exhibiting extreme behavior.r i represents the radial component, and v i the angular component.The declustering portion is relevant for data correlated in time.for ℓ = 1, . . ., d do Set b t,ℓ = F−− 1 t .With x ℓ > b t,ℓ , fit σ ℓ , ξ ℓ via MLE according to generalized Pareto likelihood.end for for i = 1, . . ., n do Define

Figure 4 .
Figure 4. Pairwise plots from ERA-Interim data after transformation and projection to S 7∞ .Down the diagonal are marginal kernel densities with two-dimensional histograms on the off-diagonal.In those plots, red indicates a higher density.All data are between 0 and 1.

χFigure 5 .
Figure 5. Pairwise extremal dependence coefficients for IVT data using the PRG-G model.

Figure 6 .
Figure 6.Conditional survival curves for selected locations using ERA-Interim and PRG-G model conditioning on all other dimensions at greater than 90th percentile (fitted).The left panel uses original units.Right panel uses standardized units.

Figure 7 .
Figure 7. Pairwise conditional survival curves for selected locations, using ERA-Interim and PRG-G models, conditioning on all other dimensions at greater than the 90th percentile (fitted).

Table 1 .
Model fit assessment and computation time on ERA-Interim and ERA5 data.(a) Energy score criterion from fitted models against the IVT data.Lower is better.(b) Time to sample (in minutes) 50,000 iterations for various models.
such that r i ≥ 1 if declustering then for i = 1, . . ., n do If r i ≥ 1 and r i−1 ≥ 1, drop the lesser (and associated v i ) from dataset.