Léo R Belzile
email: leo.belzile@hec.ca

Christophe Dutang

Paul J Northrop

Thomas Opitz

A modeler's guide to extreme value software

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

Extreme value analysis has seen strong development over the years. While software development typically lags behind methodological developments due in part to lack of recognition of the effort needed to provide reliable software, reproducibility requirements and individual efforts have led to a growth in the coverage of statistical methods. Many procedures developed in the last decades are now available, but the diversity of numerical implementations complicates somewhat the choice of routine to adopt.

Our intention, rather than to solely provide a catalog of existing software, is to discuss and compare existing implementations of statistical methods and to highlight numerical issues that are of practical importance yet are not typically discussed in theoretical or methodological papers. Our work also provides an update including the most recent software development since the reviews of [START_REF] Stephenson | Software for the analysis of extreme events: The current state and future directions[END_REF]; [START_REF] Gilleland | A software review for extreme value analysis[END_REF]; [START_REF] Gilleland | Computing software[END_REF].

Given its ongoing popularity, we focus on implementations using the R programming language, unless stated otherwise. The Comprehensive R Archive Network (CRAN) Task View on Extreme values provides an extensive list of package functionalities organized by topics; we follow this approach and broadly separate implementations into univariate, multivariate and functional extremes rather than present functionalities package by package. Using the RWsearch package [START_REF] Kiener | RWsearch: Lazy search in R packages, task views, CRAN, the web[END_REF], we automated the process of searching for extreme-related packages on the CRAN and inspected all of the packages that have "extreme value" or "peak over threshold" as keywords in the package description. Additional searches were done for unpublished packages.

Univariate extremes

Asymptotic theory for univariate extremes

The starting point for univariate extreme value analysis is the extremal types theorem. For convenience, we state the result in terms of point process representations, i.e., of the behavior of random point clouds: let X i be independent and identically distributed random variables with distribution F , where F has lower endpoint x * = inf{x : F (x) > 0} and upper endpoint x * = sup{x : F (x) < 1}. If there exist normalizing sequences a n > 0 and b n ∈ R such that the distribution of the bidimensional point process

P n = i n + 1 , X i -b n a n , i = 1, . . . , n
converges to an inhomogeneous Poisson point process on sets of the form (a, b) × (z, ∞) for 0 ≤ a ≤ b ≤ 1 and z > z * = lim n→∞ {(x *b n )/a n }, the intensity measure of the limiting process, which gives the expected number of points falling in a set, is (cf., [START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF])

Λ{(a, b) × (z, ∞)} = (b -a) 1 + ξ z-µ σ -1/ξ + , ξ = 0 (b -a) exp z-µ σ , ξ = 0 , µ ∈ R, σ > 0, (1)
where x + = max{0, x}. Equipped with this convergence statement, we can consider various tail events. For example, the limiting distribution of the maximum of n observations is

lim n→∞ Pr max n i =1 X i -b n a n ≤ x = exp [-Λ{(0, 1) × (x, ∞)}] . (2)
The right hand side of Equation ( 2) is the distribution function of the generalized extreme value distribution, say G(x), with location parameter µ ∈ R, scale parameter σ ∈ R + and shape parameter ξ ∈ R, defined on {x ∈ R : ξ(xµ)/σ > -1}. For historical reasons, the distribution is categorized based on the sign of ξ in so-called "domains of attraction". If ξ < 0, the distribution has a bounded upper tail, ξ = 0 leads to an exponential "light" tail and ξ > 0 to a "heavy tail" with polynomial decay with finite moments only of order r < 1/ξ.

We can also consider conditional exceedances of a threshold u as the latter increases to the upper endpoint x * ; for threshold exceedances over u > z * , there exists a u > 0 such that

lim u→x * Pr(a -1 u X > x + u) Pr(a -1 u X > u) = Λ{(0, 1) × (x + u, ∞)} Λ{(0, 1) × (u, ∞)} = 1 -H (x), (3)
where

1 -H (x) =    1 + ξ x σ u -1/ξ + , ξ = 0, exp(-x/σ u ) ξ = 0;
(4) and σ u = σ+ξ(u-µ). The right-hand side of Equation ( 4) is the survival function of the generalized Pareto distribution with scale σ u and shape ξ ∈ R.

Maximum likelihood estimation

We can approximate the log likelihood by taking the limiting relations of Equations ( 1) to (3) as exact for the maximum of a finite block of m observations or for exceedances above a large quantile u; the unknown normalizing constants a n , b n (respectively a u ) are absorbed by the location and scale parameters. For example, the log likelihood obtained through the limiting inhomogeneous Poisson point process for a sample of N independent observations with exceedances y 1 , . . . , y n u is (µ, σ, ξ; y) = -n u log(cσ) -

n u i =1 1 + 1 ξ log 1 + ξ y i -µ σ + -c 1 + ξ u -µ σ -1/ξ + , µ, ξ ∈ R, σ > 0. ( 5
)
The quantity c, which does not appear in Equation (1), is a tuning parameter (Coles, 2001, § 7.5): if we take c = N /m, the parameters of the point process likelihood correspond to those of the generalized extreme value distribution fitted to the maximum of blocks of m observations. This parametrization however induces strong correlation between the parameters (µ, σ, ξ).

Likelihood-based inference for extreme value distributions is in principle straightforward, even if there is no closed-form solution for the maximum likelihood estimators (MLE). The latter must solve the score equations S(θ) = 0 (i.e., the log likelihood must have gradient zero) unless ξ = -1, but with nonlinear inequality constraints depending on the sign of ξ (e.g., max(x) ≤ µσ/ξ if ξ < 0 for the generalized extreme value distribution). Constrained gradient-based optimization algorithms are thus logical choices for finding the MLE. Many numerical implementations of the log likelihood simply return very large finite values for parameter combinations outside of the support, which can impact the convergence of gradient-based optimization routines: the user is invited to check convergence of whichever software is employed by evaluating the score function at the MLE configuration to make sure it is indeed a root of the score vector S(θ). Even then, the solution returned may not be a global maximum: Figure 1 shows the conditional log likelihood surface for an inhomogeneous Poisson process model, obtained by fixing the scale. The feasible region is defined by a hyperbola and features two local maxima; depending on the starting value, a gradient algorithm would converge to different values.

Particular attention must be paid to numerical overflow when implementing the likelihood, score and information matrix of the generalized extreme value distribution, especially for terms of the form log(1+ξx) when ξ → 0 for the information and cumulants. For example, the expected information matrix satisfies I ξξ = f (ξ)/ξ -4 [START_REF] Prescott | Maximum likelihood estimation of the parameters of the generalized extreme-value distribution[END_REF] and the limit as ξ → 0 is well-defined, but this expression is numerically unstable when ξ ≈ 0. High precision functions such as log1p can be used to alleviate this somewhat, and interpolation of the cumulants when ξ ≈ 0 is recommended.

The extreme value distributions do not satisfy regularity conditions because their endpoint depends on the parameter values. Maximum likelihood estimators for both the generalized extreme value and generalized Pareto distributions fitted to respectively block maxima and threshold exceedances are nevertheless asymptotically Gaussian and the models are regular whenever ξ > -1/2 [START_REF] Smith | Maximum likelihood estimation in a class of nonregular cases[END_REF][START_REF] De Haan | On maximum likelihood estimation of the extreme value index[END_REF][START_REF] Bücher | On the maximum likelihood estimator for the generalized extreme-value distribution[END_REF]de Haan and Ferreira, 2006, § 3.4).

Moments of some of the kth order derivatives of the log likelihood exist only if the shape ξ > -1/k. Thus, when ξ ≤ -1, the MLE does not solve the score equation. The likelihood functions for the inhomogeneous Poisson point process of exceedances and the r -largest observations, the generalized Pareto, the generalized extreme value, are unbounded if ξ < -1, as there exists a combination of parameters that lead to infinite log likelihood values. This means one should restrict the parameter space to {ξ : ξ ≥ -1} and check that the solution does not lie on the boundary of the parameter space: for the generalized extreme value distribution, the conditional maximum likelihood estimator for ξ = -1 is µ = y and σ = max(y)y, whereas for the generalized Pareto distribution, σ = max(y), for the likelihood of the r -largest order statistics, σ = {max(y)y (r ) }/r, µ = max(y) -σ, etc.

The (lack of) existence of cumulants also impacts the calculation of standard errors, as elements of the Fisher information matrix are defined only if ξ > -1/2. Most software implementation compute standard errors based on the inverse Hessian matrix computed via finite-difference, but these are misleading if ξ ∈ [-1/2, -1).

We can sometimes deploy dimension reduction strategies to facilitate numerical optimization. For the generalized Pareto distribution, [START_REF] Grimshaw | Computing maximum likelihood estimates for the generalized Pareto distribution[END_REF] uses a profile likelihood to reduce the problem to a one-dimensional optimization. This is arguably one of the safest maximum likelihood estimation procedures and the exponential sub-case, for which the profile likelihood is unbounded, can be easily handled separately. The left panel of Figure 1 shows profile log likelihoods for two simulated datasets, including one for which ξ lies on the boundary of the parameter space.

As noted before, the MLE of the parameters of the inhomogeneous Poisson point process are hard to obtain because of the strong dependence between parameters: the optimization in packages such as ismev [START_REF] Heffernan | ismev: An introduction to statistical modeling of extreme values[END_REF][START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF] or evd [START_REF] Stephenson | evd: Extreme value distributions[END_REF] often fails to converge, mostly because of poor starting values and existence of several local maxima. [START_REF] Sharkey | A Poisson process reparameterisation for Bayesian inference for extremes[END_REF] propose a reparametrization that ensures approximate orthogonality of the parameters, but the following trick can also help with this: if the estimated probability of exceedance is small, the Poisson approximation implies

c 1 + ξ u -µ σ -1/ξ ≈ n u .
We can thus fit a generalized Pareto distribution to threshold exceedances, whose maximum likelihood estimates we denote ( σ u , ξ), and then use as starting values for the point-process optimization routine The black line shows one data set for which the conditional maximum ξ η = -1 lies on the boundary of the parameter space (black) and one with ξ η > -1 (grey). Right: conditional log likelihood surface for the inhomogeneous Poisson process at σ for simulated data (larger values have darker grey-scale shade). The white dot indicates the maximum likelihood estimate, while the hyperbola defines the feasible region of the parameter space given by the support constraints.

µ 0 = u -σ 0 {(n u /c) -ξ -1}/ ξ, σ 0 = σ u × (n u /c) ξ , ξ 0 = ξ.

Case study

We performed some sanity checks for various implementations of maximum likelihood estimation routines and parametric models. Specifically, we verified that density functions are non-negative and evaluate to zero outside of the domain of the distribution, and that distribution functions are non-decreasing and map to the unit interval. Certain packages have incorrect implementations of density and distribution functions.

To assess the quality of the optimization routines for extreme value distribution, we simulated exceedances and block maxima from parametric distributions with varying tail behaviors. We compared the maximum likelihood estimates returned by default estimation procedures for different packages for simulated data, checking that the value returned is a global optimum and the gradient evaluated at the value is approximately zero whenever ξ > -1. The purpose of the exercise was to check the reliability for a range of sample sizes; when systematic differences in maximum log likelihood and/or parameter estimates arose compared to other packages, then they are attributable to poor starting values, incorrect implementation of the density function, lack of handling of boundary constraints or the choice of optimization algorithm.

Figure 2 shows results for 1000 simulations, each based on n = 400 observations from a gamma distribution with shape 3 and scale 2; we set the threshold to the 0.95 theoretical quantile of the distribution, which leads to an average of 20 threshold exceedances. The left panel shows the distribution of the gradient of the log likelihood of the generalized Pareto distribution evaluated at the maximum likelihood estimate over all replicates for the shape parameter. Most instances of non-zero gradient are attributable to boundary cases, with ξ = -1 not accounted for. Other discrepancies are due to differences in numerical tolerance, but the differences in log likelihood relative to the maximum over all routines are negligible in most non-boundary cases investigated.

The right panel of Figure 2 shows the source of some oddities: the QRM package has unexpected small spread and a positive bias for estimation of ξ, different from other packages because it fails more often Figure 2: Diagnostics for maximum likelihood estimation of the generalized Pareto distribution based on 1000 samples simulated from a gamma distribution with shape 3 and scale 2 and a random number of exceedances (average 20) above the theoretical 0.95 quantile. Left: density plots of the absolute value of log gradient ∂ /∂ξ evaluated at the maximum likelihood estimator ( σ, ξ), split by simulations yielding a boundary case ( ξ = -1, blue) and regular case ( ξ > -1, orange); the y-axis scale for each package is different to ease visualization. Right: dot plots of the shape parameter estimates. Results for samples for which the numerical routines failed to converge are omitted.

when ξ is negative due to poor starting values. Both ercv and extRemes fit have noticeable point masses at ξ = 0, suggesting something is amiss. Many packages do not enforce boundary constraints and may return ξ < -1 with near-infinite gradient estimates, responsible for the large gradients in the left panel of Figure 2: Renext returns a hard-coded lower bound that is larger than -1, while SpatialExtremes and mev correctly return -1.

By contrast, the optimization routines for the generalized extreme value distribution yielded similar behavior and nearly all packages give identical results. However, we noticed that some packages fare poorly when location or scale parameters are orders of magnitude larger than scaled components: since the generalized extreme value distribution is a location-scale family, scaling the data before passing them to the routine and back-transforming after the optimization may solve such issues.

Regression modeling

Most data encountered display various forms of nonstationarity, including trends, seasonality and covariate effects, which the extreme value distributions cannot capture without modification. One can thus consider regression models in which the parameters of the extreme value distributions are functions of covariates or vary smoothly in space or time. These parameters may be suitably transformed via a link function to ensure that the functions satisfy the usual range or positivity constraints. If we assume independent observations, then maximum likelihood estimates, standard errors, etc. are obtained as before by maximizing the log likelihood function, which is now a function of the regression coefficients and of other parameters arising in the nonstationary formulation of the extreme value distribution. In models with a relatively large number of parameters, it becomes useful to include an additive penalty term in the log likelihood: for example, generalized additive models for the parameters include smooth functions (smooths in short) via basis function representations (e.g., B -splines), with a penalty that controls the wiggliness of the estimated predictor functions, which is typically evaluated using the second-order derivative of the basis functions. The obvious difficulty for numerical maximization of the log likelihood is again the presence of support constraints, since there are now potentially as many inequality constraints as there are observations. A general advice for models with covariates is that inputs should be centred and scaled to facilitate the optimization.

The ismev, texmex, eva and extRemes packages allow users to provide a model matrix (containing one covariate in each of its columns) for each parameter of the generalized Pareto and generalized extreme value distributions, thus enabling generalized linear modeling of the parameters, while the evd package only allows for linear modeling of the location parameter of the generalized extreme value distribution and bivariate counterparts; in both cases, no penalty terms are added to the log likelihood, while texmex allows for L 1 and L 2 penalties for the coefficients. The lax package supplements the functionality of these, and other, packages by providing robust sandwich estimation of parameter covariance matrix and log likelihood adjustment [START_REF] Chandler | Inference for clustered data using the independence loglikelihood[END_REF] for their fitted model objects. The GEVcdn package uses a neural network to relate the parameters of the generalized extreme value distribution with covariates [START_REF] Cannon | A flexible nonlinear modelling framework for nonstationary generalized extreme value analysis in hydroclimatology[END_REF].

The scale parameter of the generalized Pareto distribution, σ u , varies with the threshold u: it is recommendable to pay special attention to the parametrization of the scale and shape functions with covariates to ensure that the threshold stability property, which is used for extrapolation, is not lost [START_REF] Eastoe | Modelling non-stationary extremes with application to surface level ozone[END_REF]. It may be tempting to use directly the likelihood of eq. ( 5) instead (see Northrop et al., 2016). Chavez-Demoulin and [START_REF] Chavez-Demoulin | Generalized additive modelling of sample extremes[END_REF] use an orthogonal reparametrization (η, ξ), where η = σ(1 + ξ) along with bootstrap routines for uncertainty quantification; their generalized additive modeling framework is available via ismev.

Penalty terms are often used in regression models with many coefficients to control the local variability of the estimated predictors, and they allow obtaining useful estimates for coefficients whose value would otherwise not directly depend on any observations in the sample (e.g., in the case where the regression function is intended to capture spatial variability of extreme value parameters but some regions of space do not contain any observation locations). Many general packages implement generalized additive modeling with some support for extreme value distributions, including VGAM [START_REF] Yee | Vector generalized linear and additive extreme value models[END_REF]. The recent evgam package [START_REF] Youngman | evgam: An R package for generalized additive extreme value models Zhang J, Stephens MA[END_REF], dedicated to extreme value models, uses the methodology of [START_REF] Wood | Smoothing parameter and model selection for general smooth models[END_REF] for general distributions to marginalize out the regression coefficients using Laplace's method to obtain estimates of the hyperparameters (e.g., variance and autocorrelation of regression coefficients) controlling the penalty strength and shape -these are estimated simultaneously with all of the other parameters through maximum likelihood.

The evgam package builds on generic model building tools available in packages such as mgcv [START_REF] Wood | Generalized Additive Models: An Introduction with R, 2nd edn[END_REF] and provides state-of-the-art methodology tailored for extremes, including generalized additive models for extreme value distributions, quantile regression and in addition functionalities for obtaining return levels for the nonstationary models.

There usually is a Bayesian interpretation of such quadratic penalty terms as corresponding to some multivariate Gaussian prior distribution on the regression coefficients. Fitting regression (or multilevel) models is natural in the Bayesian setting, and many of the packages discussed in the next section have capabilities for fitting multilevel models.

is proportional as a function of θ to the product of the likelihood and the priors in the numerator, but the integral appearing in the denominator of Equation ( 6) is untractable in general. In such cases, the posterior density p(θ | y) usually does not correspond to any well-known distribution family, and posterior inferences about the components of θ further involve marginalizing out the other components.

For instance, to obtain the posterior density p(θ 1 | y) of the first parameter in θ, we have to evaluate the (m -1)-dimensional integral p(θ | y) d(θ 2 , . . . , θ m ). Most of the field of Bayesian statistics revolves around the creation of algorithms that circumvent the calculation of the normalizing constant (or else provide accurate numerical approximation of the latter) or that allow for marginalizing out all parameters except for one. Rather than a point estimate θ of the parameter vector, the target of inference is the whole posterior distribution. The majority of Bayesian estimation algorithms are simulation-based, and their typical output is an (approximate) sample drawn from the posterior distribution p(θ | y), from which any functional of interest can be estimated by Monte Carlo methods. Of particular interest is the posterior predictive distribution, which is obtained by simulating new observations from the response model by forwardsampling from p(y | θ (b) ) one new observation for each draw of θ (b) from the posterior.

In simple problems, exact sampling algorithms can provide independent and identical samples from the posterior. This is the exception rather than the norm; most of the time, users resort to Markov chain Monte Carlo (MCMC) algorithms for more complex settings: these algorithms admit the posterior distribution as the stationary distribution of a Markov chain with appropriately designed transition probabilities and provide autocorrelated samples from it. Another popular solution is through Laplace approximation for regression models when multivariate Gaussian priors are put on the vector of regression coefficients arising in the latent layer of the model, from which observations are conditionally independent; see the discussion in Section 2.3. In this setting, Laplace approximations give fast deterministic approximation of highdimensional integrals, which avoids resorting to simulation-based MCMC estimation. Laplace approximations are particularly accurate when they are applied twice in a certain nested way, which is known as the integrated nested Laplace approximation (INLA, Rue et al., 2009), implemented in the general INLA package offering extreme value functionality for generalized Pareto and generalized extreme value distributions.

Despite the computational overhead associated, the Bayesian paradigm has many benefits, including the capacity to incorporate physical constraints and expert opinion through the prior distributions [START_REF] Coles | A Bayesian analysis of extreme rainfall data[END_REF]. It is easier and more natural to define hierarchical structures for parameters to pool information. For multivariate and functional extremes, priors can be used for regularization purposes to pool information, for instance across time and space.

Specificity of extremes

Readers wishing to learn more about Bayesian modeling for extreme values are referred to the extensive overview in [START_REF] Stephenson | Bayesian inference for extreme value modelling[END_REF]. While Bayesian inference for extreme value models does not differ much from that of general models, additional care is required with prior specification. For example, in order to get a well-defined posterior distribution, improper reference priors such as the maximal data information (MDI) and Jeffreys priors for ξ may need to be truncated [START_REF] Northrop | Posterior propriety in Bayesian extreme value analyses using reference priors[END_REF] to result in proper (i.e., integrable) posterior distributions. [START_REF] Martins | Generalized maximum-likelihood generalized extreme-value quantile estimators for hydrologic data[END_REF] proposed using a shifted Beta distribution for ξ to constrain the support of the latter to [-0.5, 0.5]. Other popular choices are vague normal priors for location, log-scale and shape parameters, or else penalized complexity priors [START_REF] Simpson | Penalising model component complexity: A principled, practical approach to constructing priors[END_REF][START_REF] Opitz | INLA goes extreme: Bayesian tail regression for the estimation of high spatio-temporal quantiles[END_REF]. To avoid issues related to the finite and parameter-dependent lower endpoint in the generalized extreme value distribution for ξ > 0, the INLA package implements so-called blended generalized extreme value distribution that replaces the bounded lower distribution tail with the unbounded one of a Gumbel distribution through a mixture representation [START_REF] Castro-Camilo | Practical strategies for GEV-based regression models for extremes[END_REF].

Three packages, evdbayes, extRemes and MCMC4Extremes, provide MCMC algorithms for extreme value distributions, which implement so-called random walk Metropolis-Hastings steps. MCMC4Extremes (do Nascimento and Moura e Silva, 2016) is superseded by the implementation of evdbayes [START_REF] Stephenson | evdbayes: Bayesian analysis in extreme value theory[END_REF] which is more efficient, thanks to default tuning of proposal standard deviations and more flexible choices of priors. The user guide of evdbayes gives more details.

The underlying implementation of the MCMC algorithm for the function posterior in evdbayes allows for a linear trend in the location parameter and supports inclusion of positive point mass at zero for the shape parameter (Gumbel distribution for generalized extreme value distribution/exponential distribution for the generalized Pareto distributions) by leveraging an MCMC implementation using reversible jumps, following [START_REF] Stephenson | Bayesian inference for extremes: Accounting for the three extremal types[END_REF]. The gamma priors for quantile differences, used for expert prior elicitation, are provided.

Contrary to most implementations, evdbayes returns a list of posterior samples and relies on methods implemented in coda [START_REF] Plummer | CODA: Convergence diagnosis and output analysis for MCMC[END_REF] for diagnostic, summary and plots. The extRemes package [START_REF] Gilleland | extRemes 2.0: An extreme value analysis package in R[END_REF]) also has functionalities for computing posterior summaries for univariate extremes through the fevd function, which allows users to specify their own priors and proposal distributions, but the sampling is notably slower than in other packages and more cumbersome to set up, as the default values are not adequate in most cases. Linear modeling of the parameters with covariates is also possible, and Bayes factors for comparisons between models are also supported even if the methods used to compute them are not recommended. Package texmex also includes maximum a posteriori estimation and simulation from the posterior for extreme value distributions (with linear modeling of covariates) via the function evm, but only with normal priors. Behind the scenes, the texmex implementation uses an independent Metropolis-Hastings step with multivariate Cauchy/normal proposals with location vector and scale matrix based on a normal approximation to the posterior, using maximum a posteriori estimates. This translates into smaller autocorrelation (and thus larger effective sample size) than other package implementation, and it is the fastest of all MCMC implementations. The data-driven prior proposed by Zhang and Stephens (2009), reputed to give better results than maximum likelihood, is implemented in mev and is the default method for Pareto-smoothed importance sampling [START_REF] Vehtari | Practical Bayesian model evaluation using leave-one-out crossvalidation and WAIC[END_REF] from the loo package [START_REF] Vehtari | loo: Efficient leaveone-out cross-validation and WAIC for Bayesian models[END_REF]. However, because it uses the data to construct the prior, performance benchmarks alleging superior performances are misleading because of double dipping.

The current state-of-the-art method for sampling from the posterior of univariate models in simple analyses without covariates is the revdbayes package, which relies on the ratio-of-uniforms method to generate independent samples from the posterior distribution of the models. Use of advanced techniques such as mode relocation, marginal Box-Cox transformations and rotation can drastically improve the efficiency of this accept-reject scheme and make it very competitive. The ratio-of-uniforms method generates independent draws, thus avoiding the need to monitor convergence to the stationary distribution of the Markov chain and removing tuning parameters. The package uses the Rcpp interface (Eddelbuettel and Balamuta, 2018) to speed up sampling, and the sampling is an order of magnitude faster than other implementations.

While the aforementioned packages are dedicated to extreme value distributions, other popular programming languages could be used even if they would require users to implement likelihood functions themselves. Notably, the Stan programming language [START_REF] Team | RStan: the R interface to Stan[END_REF] uses Hamiltonian Monte Carlo, a state-of-the-art MCMC method, for simulating samples from the posterior distribution. The latter can easily be combined with multilevel models, but requires implementation of bespoke code for likelihood and priors that are specific to extreme value analysis. The Hamiltonian Monte Carlo sampling algorithm leads to rejection due to boundary constraints and leads to incorrect posterior draws for, e.g., the generalized extreme value distribution when ξ ≈ 0; this can be corrected by using a Taylor series approximation. The Matlab package NEVA uses a differential evolution Markov chain algorithm for estimating univariate nonstationary models [START_REF] Cheng | Non-stationary extreme value analysis in a changing climate[END_REF].

Semiparametric inference for univariate extremes

In the semiparametric approach to extremes, some components of the probability structure are handled through a relatively general (and nonparametric) asymptotic structure, which can be extrapolated towards higher yet unobserved quantile levels, for instance for the purpose of extreme-quantile estimation. The parametric form includes the shape parameter ξ and potentially second-order regular variation indices, ρ. [START_REF] Caeiro | Threshold selection in extreme value analysis[END_REF] provides a review of many estimators discussed next with an emphasis on the choice of the number of order statistics to keep for inference, which has close ties to threshold selection methods discussed in Section 2.6. Consider a sample of independent and identically distributed variables X 1 , . . . , X n with quantile function Q(x). Assuming heavy-tailed distributions with limiting shape parameter ξ > 0, the survival function

S(x) = x -1/ξ L F (x) if and only if Q(1 -1/x) = x ξ L U (x)
, with ξ the shape parameter, and L F and L U slowly varying functions, meaning lim x→∞ L(t x)/L(x) = 1 for any t > 0 (Ledford and Tawn, 1996, § 5).

Nonparametric estimators of the extreme value index are widespread, most of them variants of the [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF] estimator. Let X 1,n ≤ • • • ≤ X n,n be the order statistics of the sample of size n: the Hill estimator is the mean excess value of log-transformed data of the k largest values,

H k,n = 1 k k j =1 log X n-j +1,n X n-k,n . (7)
The Hill estimator is generally computed for a wide range of values of k, which leads to so-called Hill plots (k, H k,n ), k = 1, . . . , n. Under certain regularity conditions, the Hill estimator ( 7) is asymptotically normal (e.g., Beirlant et al., 2004, Section 4.4). A large number of R packages provides functions to estimate ( 7)

and to make Hill plots such as evir [START_REF] Pfaff | evir: Extreme values in R. R package version 1[END_REF], evmix [START_REF] Hu | evmix: An R package for extreme value mixture modeling, threshold estimation and boundary corrected kernel density estimation[END_REF], extremefit [START_REF] Durrieu | extremefit: A package for extreme quantiles[END_REF], ExtremeRisks [START_REF] Padoan | ExtremeRisks: Extreme risk measures[END_REF], ptsuite [START_REF] Munasinghe | ptsuite: Tail index estimation for power law distributions[END_REF], QRM [START_REF] Pfaff | QRM: Provides R-language code to examine quantitative risk management concepts[END_REF], ReIns [START_REF] Reynkens | ReIns: Functions from[END_REF] and tea [START_REF] Ossberger | tea: Threshold estimation approaches[END_REF].

The performance of the Hill estimator strongly depends on the number of observations kept to estimate the tail index: H k,n has a large variance if k is too small, whereas the Pareto-type tail behavior might not be verified for the selected k largest values if k is too large. The choice of k is typically based either on an empirical rule to find the area where H k,n is "stable" or by minimizing the asymptotic mean squared error (AMSE). A large number of those algorithms to minimize the latter are provided in tea along with the bootstrap methods of [START_REF] Hall | Adaptive estimates of parameters of regular variation[END_REF], [START_REF] Hall | Using the bootstrap to estimate mean squared error and select smoothing parameter in nonparametric problems[END_REF], [START_REF] Danielsson | Using a bootstrap method to choose the sample fraction in tail index estimation[END_REF], [START_REF] Caeiro | On the bootstrap methodology for the estimation of the tail sample fraction[END_REF] and [START_REF] Caeiro | Threshold selection in extreme value analysis[END_REF].

The Hill estimator (7) has a number of drawbacks: it supposes ξ > 0, it is not consistent, is not translation invariant and it behaves erratically for large k. Many alternative estimators try to palliate these lacks.

Since the Hill estimator has nondifferentiable sample paths with respect to the threshold value, the choice of threshold is notoriously difficult. [START_REF] Resnick | Smoothing the Hill estimator[END_REF] proposed a smoothed version of the Hill estimator based on averaging consecutive estimates via a moving window; these plots are provided in evmix and tea. The random block maximum estimator [START_REF] Wager | Subsampling extremes: From block maxima to smooth tail estimation[END_REF], constructed as a U statistic, has infinitely differentiable sample paths and is thus much less sensitive to the choice of k than most Hill-type estimators. The unpublished package rbm is available from Github.

Packages evt0 [START_REF] Manjunath | evt0: Mean of order p, peaks over random threshold Hill and high quantile estimates[END_REF] and ReIns implement the generalized Hill estimator based on a uniform kernel estimation [START_REF] Beirlant | Excess functions and estimation of the extreme-value index[END_REF].

evt0 also provides functions for the location-scale invariant version of the Hill estimator introduced by [START_REF] Santos | Peaks over random threshold methodology for tail index and quantile estimation[END_REF] and the biased-reduced version of [START_REF] Figueiredo | A computational study of a quasi-PORT methodology for VaR based on second-order reduced-bias estimation[END_REF], as well as a mixed moment estimator and location invariant alternative.

The package extremefit implements the kernel-weighted version of the Hill estimator of [START_REF] Grama | Statistics of extremes by oracle estimation[END_REF]; the authors provide an automatic selection procedure for the threshold u, with functions to handle these weighted estimators either for user-supplied weights or for weights automatically selected using an adaptative selection.

Moment estimators and other alternative estimators

While maximum likelihood estimation and Hill-type estimators are most commonly used for the shape parameter, other estimators are available and may be more robust in small samples. One such was proposed by [START_REF] Dekkers | A moment estimator for the index of an extreme-value distribution[END_REF] and evt0 provides a generalization of the latter by [START_REF] Santos | Peaks over random threshold methodology for tail index and quantile estimation[END_REF]. Since moments of extreme value distributions may not exist if ξ > 0, we can consider instead a bijection between the parameter vector θ and probability weighted moments of the form E[X p F (X ) q {1 -F (X )} r ] for integers p, q, r [START_REF] Hosking | Parameter and quantile estimation for the generalized Pareto distribution[END_REF]. Another avenue is to match sample linear combinations of order statistics with their theoretical counterparts using (trimmed) L-moments [START_REF] Hosking | l -moments: Analysis and estimation of distributions using linear combinations of order statistics[END_REF]. A group of R packages, including lmom [START_REF] Hosking | L-moments[END_REF], lmomco [START_REF] Asquith | lmomco -L-moments, censored $L$-moments, trimmed $L$-moments, $L$comoments, and many distributions[END_REF], TLMoments [START_REF] Lilienthal | TLMoments: Calculate TL-moments and convert them to distribution parameters[END_REF] implement these approaches for a variety of common distributions (as does the Python package lmoments), but some also allow custom distribution functions. extRemes also implements L-moments, while RobExtremes [START_REF] Ruckdeschel | RobExtremes: Optimally robust estimation for extreme value distributions[END_REF] provides robust estimators of the extreme value parameters and laeken [START_REF] Alfons | Estimation of social exclusion indicators from complex surveys: The R package laeken[END_REF] proposes robust modelling of Pareto data. Package extremeStat [START_REF] Boessenkool | extremeStat: Extreme value statistics and quantile estimation[END_REF] includes functionalities to compute extreme quantiles based on L-moments estimator.

Quantile, expectile and extremiles

In the heavy-tailed setting with ξ > 0, [START_REF] Weissman | Estimation of parameters and large quantiles based on the k largest observations[END_REF] proposed estimating the tail quantile Q(1p) for given small tail probability p using

q W k,n (p) = X n-k,n k + 1 p(n + 1) H k,n
, where H k,n is the Hill estimator (7). ReIns implements the Weissman estimator either specified by the probability level p or by the return period 1/p. The Weissman-type estimator for the class of estimators proposed by [START_REF] Santos | Peaks over random threshold methodology for tail index and quantile estimation[END_REF] are provided by evt0, whereas R package extremefit gives the quantile corresponding to weighted Hill estimator. Bias-corrected versions of the Weissman estimator also exist, yet are seemingly not implemented in software.

Quantiles can be formulated as the solution of an asymmetric piecewise linear loss function. Taking instead an asymmetric quadratic loss function yields expectiles [START_REF] Newey | Asymmetric least squares estimation and testing[END_REF], another risk measure gaining popularity in risk management [START_REF] Bellini | Risk management with expectiles[END_REF]. Many recent work studies their extremal property: on the software side, ExtremeRisks implements the methodology of [START_REF] Davison | Tail risk inference via expectiles in heavy-tailed time series[END_REF]; [START_REF] Padoan | Joint inference on extreme expectiles for multivariate heavy-tailed distributions[END_REF], including estimation of expectiles using Hill-type estimators, test of equality of tail expectiles and confidence regions for extreme expectiles. An alternative risk measures, the so-called extremiles, has been developed recently (e.g., [START_REF] Daouia | Extremile regression[END_REF] in the previous section. The underlying theory considers limiting behavior as the threshold increases.

In practice, a suitably high threshold is set empirically, balancing the bias from using a low threshold that violates the theory with statistical imprecision from using a threshold that is unnecessarily high.

For information about many of the following methods, see the review of [START_REF] Scarrott | A review of extreme value threshold estimation and uncertainty quantification[END_REF]. Methods for semiparametric estimators based on variants of Hill's estimator for the shape were presented in Section 2.5.

Visual threshold selection diagnostics

In a threshold stability plot, point and interval estimates of parameters are plotted against a range of threshold values. A particular example is the Hill plot featured in Section 2.5 (see Table 3 for an overview of available implementations). In the univariate case, the focus is often on the shape parameter ξ: we choose the lowest threshold above which we judge the point estimates of ξ to be approximately constant in threshold, bearing in mind statistical uncertainty quantified by the interval estimates. These inferences may be based on the generalized Pareto distribution (4) for threshold excesses or the inhomogeneous Poisson process model (1), using a frequentist or Bayesian analysis. In the generalized Pareto case, the threshold-independent scale parameter σ * u = σ u -ξu is used. In the frequentist case, it is useful to have the option to calculate the intervals using profile likelihoods, because they tend to have better coverage properties than Wald intervals, especially for high thresholds.

If a generalized Pareto distribution with ξ < 1 applies at threshold u then the mean excess E(Xv | X > v) is a linear function of v for all v > u. This motivates the mean residual life (MRL) plot, in which the sample mean of excesses of a range of thresholds are plotted against the threshold, with pointwise confidence intervals superimposed. We choose the lowest threshold above which the plot appears linear.

Table 4 summarises the functionality of R packages in terms of these plots.

The lmomplot function in the POT package can help to identify for which thresholds the sample L-skewness and L-kurtosis of excesses are related as expected under a generalized Pareto distribution. These plots require the use of subjective judgement to select a threshold. More formal methods seek to reduce subjectivity and perhaps introduce a greater degree of automation.

More formal methods

Penultimate models. Formal testing procedures compare the null hypothesis of having a generalized Pareto distribution above a threshold u against an alternative model. Theoretically-justified alternative models can be derived from the penultimate approximation to extremes, either by selecting piecewise constant shape [START_REF] Northrop | Improved threshold diagnostic plots for extreme value analyses[END_REF] or by using tilting function to provide more general models that should have faster convergence. The models proposed in [START_REF] Papastathopoulos | Extended generalised Pareto models for tail estimation[END_REF] lead to a threshold stability plot for an additional parameters. These approaches are implemented in mev.

Goodness-of-fit diagnostics. One drawback of the threshold stability plot and tests is that they do not entirely indicate whether the tail model fits the data well. Goodness-of-fit diagnostics can thus complement other diagnostics. The eva package [START_REF] Bader | eva: Extreme value analysis with goodness-of-fit testing[END_REF] provides multiple testing methods with the Cramér-von Mises and Anderson-Darling criteria and Moran's tests, all with control for the false discovery rate [START_REF] Bader | Automated threshold selection for extreme value analysis via ordered goodness-of-fit tests with adjustment for false discovery rate[END_REF]. The benefit of this approach, compared to visual diagnostics, is that it does not require user input and is more readily implementable with large multivariate or spatial data sets.

The approach of [START_REF] Dupuis | Exceedances over high thresholds: A guide to threshold selection[END_REF], based on examination of the weights attached to the largest observations from the sample and obtained using a robust fitting procedure, can be obtained via mev.

Sequential analysis and changepoints. Parameter estimates obtained by fitting a tail model at multiple consecutive thresholds are dependent because of the non-negligible sample overlap. The mev package provides the method of [START_REF] Wadsworth | Exploiting structure of maximum likelihood estimators for extreme value threshold selection[END_REF], which exploits a technique from sequential analysis by fitting a point process over a range of thresholds and building an approximate white noise sequence from the differences between consecutive estimates using their asymptotic covariance matrix, suitably rescaled to be standard normal. The tea package provides the Pearson χ 2 test of normality applied to sequences of differences of scale estimates, following [START_REF] Thompson | Automated threshold selection methods for extreme wave analysis[END_REF], while threshold stability plots based on estimates of the coefficient of variation and sequential testing of del Castillo and Padilla (2016) are included in ercv.

Gerstengarbe-Werner plots (cf., Cebrián et al., 2003, Appendix B) are graphical diagnostic plots derived from consecutive differences between order statistics. Sequential Mann-Kendall tests are performed from smallest to largest order statistics (and vice-versa), with the intersection point used as a changepoint candidate; such plots can be created with tea.

Predictive performance. The threshr package [START_REF] Northrop | Cross-validatory extreme value threshold selection and uncertainty with application to ocean storm severity[END_REF] looks at the predictive performance of the generalized Pareto for a binomial-generalized Pareto model fitted using the Bayesian approach. The scheme uses a leave-one-out cross validation scheme for values at a fixed validation threshold v at or above the range of potential thresholds considered.

Mixture models. The generalized Pareto specifies a distribution only for exceedances above a threshold u, but having a model below this threshold may be desirable, with some options enabling automatic threshold selection. The evmix package [START_REF] Hu | evmix: An R package for extreme value mixture modeling, threshold estimation and boundary corrected kernel density estimation[END_REF] provides implementations of most of the mixture models listed in [START_REF] Scarrott | A review of extreme value threshold estimation and uncertainty quantification[END_REF]: this includes parametric models for the bulk of the data (for which users can inform threshold selection by looking at the profile likelihood for u), nonparametric and kernel-based approaches for the data below the threshold. Many such models are discontinuous at the threshold and require choosing a fixed threshold. The extremefit package [START_REF] Durrieu | extremefit: A package for extreme quantiles[END_REF] provides a mixture model implementation with a kernel-based bulk model and adaptive selection rules for the bandwidth parameter. The mev package provides the extended generalized Pareto model of [START_REF] Naveau | Modeling jointly low, moderate, and heavy rainfall intensities without a threshold selection[END_REF] for modeling rainfall.

Univariate extremes implementations in other programming languages

While R is arguably the programming language boasting the most software implementations used for extreme value analyses, some basic routines are available elsewhere for estimation of univariate models using maximum likelihood or probability weighted moments: these include the Julia package Extremes, the Matlab EVIM package and the Python packages wafo, pyextremes and scikit-extremes.

Multivariate extremes

The lack of ordering of R D leads to multiple definitions of extremes [START_REF] Barnett | The ordering of multivariate data (with discussion)[END_REF]. We focus on componentwise maxima and concomitant exceedances, which lead to the multivariate analog of block maximum and peaks over threshold methods. Another option, structure variables, reduces the data to univariate summaries and can be dealt with using tools presented before.

Multivariate maxima

Consider independent and identically distributed sequence of D-variate random vectors {Y i } i ≥1 , where each vector Y i has marginal distribution functions F j ( j = 1, . . . , D). Paralleling the univariate case, we consider the vector of componentwise maxima M n = (M n,1 , . . . , M n,D ), where M n, j = max{Y 1, j , . . . , Y n, j }. If there exists sequences of location and scale vectors

a n ∈ R D + and b n ∈ R D such that lim n→∞ Pr a -1 n (M n -b n ) = G(y),
for G non-degenerate, then the limiting distribution is of the form

G(y) = exp -D S D max w z dH (w ) , (8)
where the so-called spectral measure H is a probability measure on the D-simplex S D = {ω ∈ R D + : ω 1 = 1} and the margins are generalized extreme value distributed. Distributions of the form eq. ( 8) are termed multivariate extreme value, or max-stable. The simple max-stable distribution is obtained upon setting µ = 1, σ = 1, ξ = 1 for each margin, corresponding to the unit Fréchet.

The lack of specification for H , other than the D moment constraints E(S j ) = 1/D ( j = 1, . . . , D) for S ∼ H , means that the set of probability measures satisfying the moment constraint is infinite, unlike in the univariate case. The copula package includes three tests of the max-stability assumption; see [START_REF] Kojadinovic | Large-sample tests of extreme-value dependence for multivariate copulas[END_REF]; [START_REF] Kojadinovic | Nonparametric rank-based tests of bivariate extreme-value dependence[END_REF]; Ben [START_REF] Ghorbal | On the Ghoudi, Khoudraji, and Rivest test for extreme-value dependence[END_REF], while the graphical diagnostic proposed by [START_REF] Gabda | Discussion of "Statistical modeling of spatial extremes[END_REF] is part of mev. Likelihood-based estimation. The likelihood of a simple max-stable random vector Z with a parametric model for V is obtained by differentiating the distribution function exp{-V (z)} with respect to each z 1 , . . . , z D .

The number of terms in the likelihood is the Dth Bell number, which is the total number of partitions of D elements into k (k = 1, . . . , D) elements. Even in moderate dimensions, the number of distinct likelihood contributions is huge and the calculations become prohibitive. One way to circumvent this problem is to add the information about the partition if occurrence times are recorded [START_REF] Stephenson | Exploiting occurrence times in likelihood inference for componentwise maxima[END_REF]; this is implemented in the evd, but only for bivariate models. In practice, we replace the limiting partition with the empirical one. The likelihood is biased unless n D since the empirical partition also needs to converge to the limiting hitting scenario; for weakly dependent processes, use of the observed partition may induce bias [START_REF] Wadsworth | On the occurrence times of componentwise maxima and bias in likelihood inference for multivariate max-stable distributions[END_REF]. Instead, [START_REF] Thibaud | Bayesian inference for the Brown-Resnick process, with an application to extreme low temperatures[END_REF] propose to impute the partition using a Gibbs sampler, while Huser et al. (2019) use a stochastic expectation-maximisation algorithm; the E -step for the missing partition uses a Monte-Carlo estimator, where approximate draws are obtained from the Gibbs sampler of [START_REF] Dombry | Conditional simulation of max-stable processes[END_REF]. None of these extensions have been implemented in publicly available software packages. Parametric models. While max-stable models have been around for a while, there are few software implementations for estimating such models. The evd and copula packages provide functionalities that are restricted to the bivariate setting. The SpatialExtremes and CompRandFld packages have methods for fitting max-stable processes using pairwise composite likelihood for spatial models; see Section 4.

There are only handful of useful parametric models that generalize to dimension D > 2. The prime example is the logistic multivariate extreme value model, which is overly simplistic and lacks flexibility since the distribution is exchangeable. Many existing models are special cases of a Dirichlet family of distributions [START_REF] Belzile | Extremal attractors of Liouville copulas[END_REF] and obtained through tilting [START_REF] Coles | Modelling extreme multivariate events[END_REF] to satisfy the moment constraint. These all have the drawback that the number of parameters is constant or grows linearly with the number of dimensions O(D) and this typically isn't enough for characterizing complex data. Two models derived from elliptical distributions, the Hüsler-Reiss model [START_REF] Hüsler | Maxima of normal random vectors: between independence and complete dependence[END_REF] and the extremal Student-t [START_REF] Nikoloulopoulos | Extreme value properties of multivariate t copulas[END_REF], are more useful in large dimensions because their scale matrix can be used to parametrize the pairwise dependence individually with O(D 2 ) entries, and they can be more readily adapted to the functional setting, with extensions for skew-symmetric families [START_REF] Beranger | Models for extremal dependence derived from skew-symmetric families[END_REF]. The last parametric family, of which the most prominent example is the asymmetric logistic distribution, are max-mixtures [START_REF] Stephenson | Exploiting occurrence times in likelihood inference for componentwise maxima[END_REF] that assign different weights to multiple simultaneous combinations of extremes. This allows for some degree of asymmetry and asymptotic independence, but such models are overparametrized with O(2 D ) coefficients.

Joint estimation of all marginal and dependence parameters is complicated because of the potential high-dimensionality of the optimization problem, but also because of potential model misspecification that leads to unplausible parameter estimates. It is therefore common to use a two-stage approach, whereby data are first transformed to standardized margins and then dependence parameters are fitted separately. The function fbvevd in evd allows the user to pass fixed values for some parameters. The tailDepFun package contains routines for fitting the continuous updating weighted least squares estimator, along with goodness-of-fit tests, for multivariate and functional models including max-linear models [START_REF] Einmahl | A continuous updating weighted least squares estimator of tail dependence in high dimensions[END_REF].

A different avenue is to estimate an equivalent form of H termed the Pickands dependence function (cf., Falk et al., 2011, p. 150), which is equivalent to eq. ( 8).

Enforcing properties of the Pickands dependence function, notably convexity and known values on the corners of the simplex, can improve estimation.

The evd package allows users to estimate nonparametrically the bivariate dependence function based on the estimators of [START_REF] Pickands | Multivariate extreme value distributions[END_REF] and [START_REF] Caperaa | A non-parametric estimation procedure for bivariate extreme value copulas[END_REF] for block maxima; additional options correct for boundary and convexity constraints.

While most estimators are restricted to the bivariate case, copula with its function An provides generalization to higher dimensions for multiple estimators [START_REF] Gudendorf | Nonparametric estimation of multivariate extreme-value copulas[END_REF]. Multivariate estimators based on Bernstein polynomials that guarantee convexity [START_REF] Marcon | Multivariate nonparametric estimation of the Pickands dependence function using Bernstein polynomials[END_REF] are provided by the beet procedure in ExtremalDep, along with the madogram estimator. Bayesian estimation is also available in the bivariate case, imposing a prior on the order of the Bernstein polynomials. The package also includes a procedure for computing pointwise confidence intervals using a nonparametric bootstrap. The UniExtQ from ExtremalDep provides credible intervals for bivariate extreme quantile regions [START_REF] Beranger | Estimation and uncertainty quantification for extreme quantile regions[END_REF], estimated using an extension of this approach. Lastly, fCopulae provides parametric dependence function, correlation coefficient and tail dependence measures for bivariate extreme value copulas. Unconditional simulation algorithms. For a long time, exact unconditional simulation algorithms for max-stable processes were elusive outside of special cases [START_REF] Schlather | Models for stationary max-stable random fields[END_REF]. Both mev and graphicalExtremes implement the algorithm of [START_REF] Dombry | Exact simulation of max-stable processes[END_REF] for selected multivariate models (including for the latter extremal graphical models on trees) ensuring exact simulation, whereas evd uses dedicated algorithms for logistic and asymmetric logistic models in arbitrary dimensions [START_REF] Stephenson | Simulating multivariate extreme value distributions of logistic type[END_REF]. The copula (evCopula objects) [START_REF] Yan | Enjoy the joy of copulas: With a package copula[END_REF] and SimCop packages [START_REF] Tajvidi | A general approach to generate random variates for multivariate copulae[END_REF] have functionalities for simulation of some bivariate extreme value distributions and the multivariate logistic model, or Gumbel copula. Packages mev and BMAmevt provide random number generators for selected parametric angular density models.

Threshold models

Multivariate regular variation, which underlies the max-stable distribution of Equation ( 8), can also be used for threshold exceedances by considering the associated Poisson point process of extremes with intensity measure Λ on a risk region R ⊂ R D + \ {0 D }, i.e., the positive orthant excluding the origin [START_REF] Resnick | Extreme values, regular variation, and point processes, Applied Probability[END_REF]. Assuming the intensity measure is absolutely continuous, the intensity function λ

(x) = ∂ D Λ(x)/(∂x 1 • • • ∂x D )
exists and we can define a density over R by renormalizing λ(x) by the measure of the risk region Λ(R). The resulting likelihoods of the point process, multivariate generalized Pareto distributions and more general threshold models are much simpler than their max-stable counterpart, but there are typically two numerical bottlenecks associated to fitting these models. The first arises from the calculation of the measure of the risk region, which is often intractable and must thus be estimated using Monte Carlo methods. There are closed-form expressions for the special case

R = {x ∈ R D + : x i > u}; if ξ = 1 D , then R = {x ∈ R D + : X 1 > u} has risk measure Λ(R) = Du -1 irrespective of the model for Λ.
The second bottleneck is due to censoring: not all components of a random vector may be extreme and the limiting model may be a poor approximation at finite levels for weakly dependent vectors [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF]. To reduce the bias arising from consideration of the asymptotic distribution, it is customary to left-censor observations falling below marginal thresholds. Most multivariate peaks over threshold models are based on the multivariate generalized Pareto [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF], whereby R = {x ∈ R D + : max D j =1 x j > u} -other constructions are described in [START_REF] Rootzén | Multivariate peaks over thresholds models[END_REF]. [START_REF] Kiriliouk | Peaks over thresholds modeling with multivariate generalized Pareto distributions[END_REF] provide expressions for the likelihood of many parametric models with strategies for diagnostics; these are not currently implemented in software. The point process likelihood can also be used in place of the multivariate generalized Pareto: in the bivariate case, the evd package proposes the point process likelihood [START_REF] Smith | Markov chain models for threshold exceedances[END_REF], but the bivariate censored likelihood implemented therein actually uses the max-stable copula instead [START_REF] Ledford | Statistics for near independence in multivariate extreme values[END_REF].

Most implementations are restricted to the bivariate setting or are reserved for spatial data. The graphicalExtremes package (Engelke and Hitz, 2020) is a notable exception: it implements the multivariate Hüsler-Reiss generalized Pareto distribution for graphical models. Exploiting the relation between the model and conditional extremal dependence, the parameters of the Hüsler-Reiss or Brown-Resnick process are directly related to the variogram matrix, whose entries are estimated empirically using pairwise empirical estimators of χ. The full likelihood can be used (including censoring), but the factorization of the likelihood over cliques allows for higher-dimensional models to be fitted through maximum likelihood at reasonable cost, since each component is low dimensional. The bvtcplot function in the evd package provides threshold stability plots in the bivariate case based on the spectral measure. mev provides composition sampling algorithms for threshold models for various risk functionals R in the multivariate setting [START_REF] Ho | Simple models for multivariate regular variation and the Hüsler-Reiss Pareto distribution[END_REF].

Rather than condition on the maximum component exceeding a threshold, we can focus instead on exceedances of the j th component, i.e., consider a limiting model for Y -j | Y j > u. [START_REF] Heffernan | A conditional approach for multivariate extreme values (with discussion)[END_REF] showed that a particular choice of normalizing sequences allows for the existence of nondegenerate limiting measure, including for asymptotically independent models. Inference for the conditional extremes model is usually performed in two stages. In the first, the marginal distributions are estimated semiparametrically and data are transformed to Laplace margins [START_REF] Keef | Estimating the probability of widespread flood events[END_REF]. In the second step, the dependence parameter vectors are estimated using a nonlinear regression model under the assumption of Gaussian residuals.

Inference for the conditional extremes model as implemented in the texmex package relies on simulation: the probability of extreme events is obtained by calculating the fraction of simulated points falling in the risk region and uncertainty quantification is done using the bootstrap scheme described in [START_REF] Heffernan | A conditional approach for multivariate extreme values (with discussion)[END_REF].

The multivariate regular variation representation provides another modeling approach for peaks over threshold using radial exceedances. For this, the data are first transformed to standardized scale with unit shape and the resulting components Y are mapped to pseudo-angles (R, Θ), with, e.g., R = Y 1 and Θ = Y -D /R. Since R and Θ become stochastically independent as R tends to infinity, one can focus on modelling the spectral measure H (θ) appearing in eq. ( 8). ExtremalDep supports composite likelihood maximum estimation with pseudo-angles for D-dimensional distributions, with composite likelihood information criterion estimates to compare models, with density functions, plots and associated functions for each parametric model. Nonparametric estimation of the spectral measure only requires the user to impose mean contraints. Starting from a sample of pseudo-angles, these can be enforced through empirical likelihood method [START_REF] Einmahl | Maximum empirical likelihood estimation of the spectral measure of an extreme-value distribution[END_REF] or Euclidean likelihood [START_REF] De Carvalho | A Euclidean likelihood estimator for bivariate tail dependence[END_REF]. The extremis package [START_REF] De Carvalho | extremis: Statistics of extremes[END_REF] implements these functionalities in the bivariate setting, and mev in higher dimensions. The unpublished EVcopula package implements the bivariate model of Wadsworth ( 2016) along with likelihood-based estimation methods and can be used to estimate probabilities of large bivariate quantiles for both asymptotic (in)dependence scenarios. The BMAmevt package is dedicated to the implementation of a Bayesian nonparametric model that uses a trans-dimensional Metropolis algorithm for fitting a Dirichlet mixture to the spectral measure based on pseudo-angles in moderate dimensions [START_REF] Sabourin | Bayesian Dirichlet mixture model for multivariate extremes: A reparametrization[END_REF].

Coefficients of tail dependence and structural variables

In multivariate settings, knowing the speed of decay of the dependence between pairs of random variable is useful for risk assessment. This also helps validate empirically if asymptotic multivariate extreme value models that assume multivariate regular variation are warranted or not. The tail correlation coefficient [START_REF] Coles | Dependence measures for extreme value analyses[END_REF],

χ(v) = Pr[min i {F i (Y i ) > v}] 1 -v , ( 9
)
with χ = lim v→1 χ(v), is used to assess whether extremes are asymptotically independent (χ = 0) or dependent (χ > 0). Equation ( 9) suggests replacing the unknown distribution functions by their empirical counterpart. The bivariate empirical estimator is usually defined as 2 -log{ C (v1 2 )/ log(v) rather than

{1 -C (v, v)}/(1 -v) by making the approximation 1 -v ∼ -log(v), 1 -C (v, v) ∼ -log{C (v, v)} for v 1.
A related coefficient measuring dependence is the coefficient of tail dependence, often denoted η. With data transformed to unit Pareto margins, say Y p , the structural variable T p = min D j =1 Y p j is such that, for large u (Ledford and Tawn, 1996, eq. 5.6),

Pr T p > u + t | T p > u ∼ L(u + t ) L(u) (1 + t /u) -1/η , ( 10
)
with L(x) a slowly varying function. The η coefficient can be estimated by fitting a generalized Pareto distribution with shape η and scale ηu to exceedances of T p above u. If data are transformed to the exponential scale instead, the scale parameter of the structural variable is η and the maximum likelihood estimator of the latter coincides with Hill's estimator (Section 2.5).

The coefficient of tail dependence takes values in (0, D -1 ) if the variables are negatively associated, η = D -1 for independent variables, and η ∈ (D -1 , 1] if the variables exhibit positive association. In the multivariate setting, the coefficients η C for subsets C ⊂ {1, . . . , D} satisfy ordering constraints with η C ≤ η [START_REF] Schlather | A dependence measure for multivariate and spatial extreme values: Properties and inference[END_REF].

In the bivariate setting, it is customary to consider χ = 2η-1 instead of η, which gives χ ∈ (-1, 1] [START_REF] Coles | Dependence measures for extreme value analyses[END_REF]. The evd package function chiplot provides plots of χ and χ based on the empirical distribution of the minimum, with approximate pointwise standard errors through the delta-method. package provides various estimators of η and χ, while graphicalExtremes includes empirical estimators emp_chi that can be used to obtain empirical estimates of the dependence matrix of the Hüsler-Reiss distribution.

Extensions that consider different tail decays have emerged in the last decade: for example, [START_REF] Beirlant | Bias-reduced estimators for bivariate tail modelling[END_REF] and [START_REF] Dutang | Robust and bias-corrected estimation of the coefficient of tail dependence[END_REF] consider projections of the form Z ω = min{Y p 1 , Y p 2 ω/(1 -ω)} for ω ∈ (0, 1) a fixed angle. Under a regular variation assumption, the distribution of Z ω can be approximated by the so-called extended Pareto distribution. The parameters of the latter can be estimated using the minimum density power divergence (MDPD) criterion [START_REF] Dutang | Robust and bias-corrected estimation of the coefficient of tail dependence[END_REF], which includes the maximum likelihood estimator as a special case.

The RTDE package [START_REF] Dutang | RTDE: Robust tail dependence estimation[END_REF] provides various functions to estimate the parameters of this model, and the returned objects allow users to summarize/plot fitted outputs, to compute the bivariate tail probability as well as to perform a simulation analysis. A similar approach is considered in [START_REF] Wadsworth | A new representation for multivariate tail probabilities[END_REF] and implemented in lambdadep function of the mev package; the authors look at different extrapolation paths by replacing the multivariate regular variation by a collection of univariate regular variation assumptions. Mhalla et al. ( 2019) also use such ideas to implement generalized additive regression for extremal dependence parameters.

The drawback of these approaches, termed structural variables since they use univariate projections, is that estimation is carried independently for every angle ω.

Time series and graphical models

Data on a single variable collected over time often exhibit short-term temporal dependence, which can lead to extremes occurring in clusters. As a minimum, statistical methods for time series extremes need to account for dependence in the data and to estimate the extent to which extremes cluster, either directly or using a dependence model. For reviews of this area see Chavez-Demoulin and Davison (2012) and [START_REF] Reich | Bayesian inference for extreme value modelling[END_REF].

Extremal index estimation

For stationary processes satisfying the D(u n ) condition, which limits long-range dependence at extreme levels, the strength of local serial extremal dependence is commonly measured by the extremal index. The latter can be interpreted as the reciprocal of the limiting mean cluster size in a Poisson cluster process of exceedances of increasingly high thresholds. Table 6 gives basic information about the direct estimators of the extremal index that feature in this section, while Table 7 summarises implementations of these estimators, including information about diagnostics for the choice of tuning parameters. When a threshold is involved these diagnostics can be used for threshold selection. The diagnostics in the evd, evir, exdex, fExtremes and texmex packages are threshold stability plots for the extremal index. The information matrix test of [START_REF] Süveges | Model misspecification in peaks over threshold analysis[END_REF], which is based on a model for truncated inter-exceedance times called K -gaps, is provided by the mev and exdex packages. The packages evd (function clusters), ex-tRemes (decluster), fExtremes (deCluster), POT (clust) and texmex (declust) use an estimate of the extremal index to decluster exceedances of a threshold to form a series of sample cluster maxima.

Marginal modeling

Suppose that interest is limited to marginal extremes. The limiting distributions of cluster maxima and a randomly chosen threshold exceedance are identical, so inferences can be made using a marginal generalized Pareto model for sample cluster maxima or for all exceedances. The texmex [START_REF] Southworth | texmex: Statistical modelling of extreme values. R package version 2[END_REF] package is the most complete implementation of the analysis of cluster maxima: it uses a semiparametric bootstrap procedure to account for uncertainty in declustering and in marginal inference and can also accommodate covariate effects. The declustering approach is wasteful of data and [START_REF] Fawcett | Estimating return levels from serially dependent extremes[END_REF] show that the difficulty of identifying clusters reliably can lead to substantial bias. When using all exceedances appropriate adjustment must be made for dependence in the data and for the value of the extremal index [START_REF] Fawcett | Estimating return levels from serially dependent extremes[END_REF]: the lite package uses the methodology of estimator reference tuning parameter(s) runs [START_REF] Smith | Estimating the extremal index[END_REF] run length, threshold blocks (blocks 1) [START_REF] Smith | Estimating the extremal index[END_REF] block size, threshold modified blocks (blocks 2) [START_REF] Smith | Estimating the extremal index[END_REF] block size, threshold intervals (FS) [START_REF] Ferro | Inference for clusters of extreme values[END_REF] threshold iterative least squares (ILS) Süveges ( 2007) threshold K -gaps [START_REF] Süveges | Model misspecification in peaks over threshold analysis[END_REF] run length K , threshold semiparametric maxima (SPM) [START_REF] Northrop | An efficient semiparametric maxima estimator of the extremal index[END_REF] block size 8: Overview of packages and main functions for modeling time series extremes by area: marginal modelling (m); exploratory analysis (e); dependence modelling (d). [START_REF] Chandler | Inference for clustered data using the independence loglikelihood[END_REF] to estimate a marginal log likelihood that has been adjusted for clustering using a sandwich estimator of the covariance matrix of the marginal parameters and combines this with a log likelihood for the extremal index under the K -gaps model. The extremefit package provides a semiparametric procedure for time series extremes, as described in Section 2.5. Table 8 gives summaries of these packages and the packages that enable the estimation of time series dependence.

Models for dependence

In some applications it is important to infer more about the behavior of an extreme event than the size of a cluster of extreme values. For example, the duration of an extreme event or an accumulation of the extreme values may be of interest. This requires the nature of serial extremal dependence to be modeled.

The extremogram package implements the extremogram [START_REF] Davis | The extremogram: A correlogram for extreme events[END_REF][START_REF] Davis | Estimating extremal dependence in univariate and multivariate time series via the extremogram[END_REF][START_REF] Davis | Towards estimating extremal serial dependence via the bootstrapped extremogram[END_REF] to inform modeling by exploring quantitatively serial extremal dependence within stationary time series and between different time series. In the univariate case, it gives estimates of the conditional probabilities that a variable exceeds a user-supplied high threshold at time t + l given that it exceeded this threshold at time t . The stationary bootstrap is used to provide confidence intervals.

The fitmcgpd function in the POT package performs maximum likelihood inference using a firstorder Markov chain model, in which one of several bivariate extreme value distributions is used as a model for successive threshold exceedances [START_REF] Smith | Markov chain models for threshold exceedances[END_REF]. The function simmc simulates from this type of model, as does the evmc function in the evd package. The tsxtreme package models time series dependence using the conditional extremes approach of [START_REF] Heffernan | A conditional approach for multivariate extreme values (with discussion)[END_REF], which enables a greater range of dependence structures to be modeled. Inferences are performed using two-step maximum likelihood fitting and a Bayesian approach in which inferences are made about a more flexible model in which all inferences are performed simultaneously [START_REF] Lugrin | Bayesian uncertainty management in temporal dependence of extremes[END_REF]. The functions theta2fit (MLE) and thetafit (Bayesian) provide inferences for the sub-asymptotic extremal index of [START_REF] Ledford | Diagnostics for dependence within time series extremes[END_REF].

The ev.trawl package implements the modeling approach described in [START_REF] Noven | A latent trawl process model for extreme values[END_REF], which is based on the representation of a generalized Pareto distribution as a mixture of exponential distributions in which the exponential rate has a gamma distribution. An exponential trawl process introduces time series dependence in a latent gamma process, while a marginal probability integral transform allows both negative and positive shape parameter values. The CTRE package deals with processes for which interexceedance times have a heavy-tailed distribution and therefore a Poisson cluster representation is not appropriate [START_REF] Hees | Statistical inference for inter-arrival times of extreme events in bursty time series[END_REF]. Parameter stability plots are provided as a means of selecting a suitable threshold.

Graphical extremes

Under the first-order Markov chain model for time series extremes of [START_REF] Smith | Markov chain models for threshold exceedances[END_REF], the value of a variable at time t + 1 is assumed to be conditionally independent of its value prior to time t given the value at time t . This simple dependence structure could be represented as a graphical model in which nodes representing the value of the variable are only connected by an edge if they correspond to adjacent time points.

The packages graphicalExtremes [START_REF] Engelke | Graphical models for extremes (with discussion)[END_REF] and gremes [START_REF] Asenova | Inference on extremal dependence in the domain of attraction of a structured Hüsler-Reiss distribution motivated by a Markov tree with latent variables[END_REF] provide more general graphical modeling frameworks for extremes, based on a multivariate Hüsler-Reiss generalized Pareto model for peaks over thresholds; see also Section 3.2. A graph represents conditional independences between variables. If the graph is sparse then the joint distribution decomposes into the product of lower-dimensional distributions, which results in a more parsimonious and tractable model. If the graph is a tree, that is, there is exactly one path along edges between any pair of nodes, then this decomposition is particularly simple. The graphicalExtremes and gremes packages provide functions to fit a multivariate Hüsler-Reiss generalized Pareto model given a user-supplied graph and functions to simulate from this model. The specifics of the theory underlying these packages differ but the resulting model structures coincide when based on a tree.

In some applications, such as the analysis of extreme river flows, there is a physical network from which the graph can be constructed. In other cases the graph is conceptual: graphicalExtremes also provides a means to infer the structure of a graph from data.

Functional extremes (including spatial extremes)

Functional extremes designates a relatively recent branch of extreme value analysis concerned with stochastic processes over infinite-dimensional spaces, especially spatial and spatio-temporal extremes in geographic space [START_REF] Davison | Statistical modeling of spatial extremes[END_REF][START_REF] Huser | Advances in statistical modeling of spatial extremes[END_REF]. We here use the term space for R d with d ≥ 1, including the combination of geographic space and time (d = 3), and we explicitly refer to time only where necessary. In practice, we usually work with finite discretizations of the study domain, such that many multivariate results and techniques carry over to the functional setting, although usually in relatively high dimension.

Common exploratory tools for extremal dependence are coefficients for bivariate distributions assessed as a function of spatial distance or temporal lag (e.g., extremal coefficient function based on bivariate extremal coefficients θ 2 , tail correlation function based on the χ measure, F -madogram, concurrence probability for maxima).

The asymptotic mechanisms for functional maxima and threshold exceedances are similar to the multivariate setting. Available statistical implementations are summarized in Section 4.1. Marginal and dependence modeling is discussed in Section 4.3. Aspects that we consider as still underdeveloped in existing implementations are listed in Section 4.4.

We use notation X = {X (s)} for stochastic processes indexed by s ∈ S ⊂ R d , representing the process of the original event data. Usually we have a dataset of observations X i (s j ) for j = 1, . . . , m locations observed at i = 1, . . . , n time points.

Max-stable processes are the natural class of models for locationwise maxima taken over temporal blocks of the same length, such as annual maxima observed at fixed spatial locations. A max-stable process possesses finite-dimensional max-stable distributions, and convergence to a max-stable process can be defined through the convergence of all finite-dimensional distributions, such that strong links arise with the univariate and multivariate setting. If there exist sequences of normalizing functions a n (s) > 0 and b n (s) such that

a n (s){M n (s) -b n (s)} → Z (s), s ∈ S , n → ∞, (11)
with Z (s) a nondegenerate limit process, then Z (s) is max-stable. Generalized r -Pareto processes [START_REF] Ferreira | The generalized Pareto process; with a view towards application and simulation[END_REF][START_REF] Thibaud | Efficient inference and simulation for elliptical Pareto processes[END_REF][START_REF] Dombry | Functional regular variations, Pareto processes and peaks over threshold[END_REF][START_REF] De Fondeville | High-dimensional peaks-over-threshold inference[END_REF][START_REF] Engelke | Extremal behaviour of aggregated data with an application to downscaling[END_REF] arise asymptotically when a sum-mary functional r (X ) of the process X exceeds a threshold that tends towards the upper endpoint of the probability distribution of the functional r .

Limit theory and statistical methodology in this peaks over threshold setting was formulated under the assumption that realizations of X correspond to continuous functions over a compact domain S . The most widely used setting for functional peaks over threshold follows the multivariate setting by assuming that data have been standardized, i.e., marginally transformed with a transformation T that is strictly monotonic (i.e., T (x 2 ) > T (x 1 ) if x 2 > x 1 ), and that ensures positivity (i.e., T (x) ≥ 0) with standardized tails of transformed random variables for which lim x→∞ x Pr[T {X (s)} > x] = 1. We can choose T s {X (s)} = 1/[1 -F s {X (s)}], where F s is the marginal distribution of X (s). Typically, summary functionals r are homogeneous (i.e., r (t x) = t r (x) for t > 0) ; examples include the average r (x) = |S | -1 S x(s) ds, the minimum r (x) = min s∈S x(s), the maximum r (x) = max s∈S x(s), or the median. Convergence is assumed in the space of continuous functions over compact S , such that the distribution of the functional r [T s {X (s)}] is well defined. Functional convergence of maxima in (11) implies functional convergence to r -Pareto processes Y r = {Y r (s)}:

T s {X (s)} | r [T s {X (s)}] ≥ u → {Y r (s)} u → ∞. (12)
Max-stable and generalized Pareto processes have different probabilistic structures, but there always is a one-to-one correspondence between their dependence structures. Estimation of the marginal distributions and of the dependence structure is often conducted in two separate steps. The space-time dependence between sites is normally captured by correlation functions or variograms, which leads to much fewer parameters to infer than in the unstructured multivariate setting.

These asymptotic models can accommodate either asymptotic dependence or full independence among the variables X (s 1 ) and X (s 2 ) at locations s 1 , s 2 ∈ S . The coefficient of tail dependence introduced in (10), if considered for D sites s 1 , . . . , s d , is therefore restricted to values η ∈ {1/D, 1}. More flexible dependence structures can be achieved within the conditional extremes framework with conditioning on a fixed location [START_REF] Wadsworth | Higher-dimensional spatial extremes via single-site conditioning[END_REF][START_REF] Simpson | High-dimensional modeling of spatial and spatio-temporal conditional extremes using INLA and the SPDE approach[END_REF]. Finally, so-called subasymptotic models do not arise as classical extreme value limits but focus on flexibly capturing dependence remaining at subasymptotic levels, for instance with asymptotic independence where 1/D < η(s 1 , s 2 ) < 1 is possible;

for example, the class of max-infinitely divisible processes [START_REF] Huser | Max-infinitely divisible models and inference for spatial extremes[END_REF], which is useful for flexible modeling of location-wise maxima. Most such proposals do not come with packaged and generic software implementations so far.

Max-stable processes for maxima data

Suppose that data consist of locationwise block maxima M i (s j ), where i = 1, . . . , n block indexes the blocks, e.g., the observation year in case of annual block maxima. The SpatialExtremes package provides the most comprehensive collection of functions for exploration and statistical inference with max-stable processes for spatial maxima data in geographical space (d = 2). While standard full likelihoods are not tractable even for moderately many locations with the common models, pairwise likelihood has become the standard approach for fitting max-stable processes, with implementations in SpatialExtremes and CompRandFld. Global dependence measures such as concurrence maps [START_REF] Dombry | Probabilities of concurrent extremes[END_REF], available from concurrencemap in SpatialExtremes, can be constructed from bivariate summaries. The intractability of the multivariate max-stable distribution function G, described in eq. ( 8), has led to pairwise likelihood becoming the standard estimation method for spatial max-stable processes. In SpatialExtremes, joint frequentist estimation of marginal and dependence parameters is possible, where auxiliary variables can be flexibly included in the three parameters of the marginal generalized extreme value distribution. Similar to generalized additive models, smoothness penalties can be imposed on nonlinear effects modeled through spline functions. In contrast to the aforementioned generalized additive model approach without dependence, the numerical optimization becomes more involved here, such that only a moderate number of marginal parameters can be reasonably estimated.

RandomFields [START_REF] Schlather | Analysis, simulation and prediction of multivariate random fields with package RandomFields[END_REF] provides a large variety of max-stable models and particularly of tail correlation functions, with a focus on implementing simulation from such models. Moreover, the package encapsulates vast functionality, especially simulation, for Gaussian random fields, which are often the building blocks for the more sophisticated extreme value models. The package provides multiple state-of-the-art algorithms for simulating Brown-Resnick max-stable processes. Exact unconditional simulation of max-stable processes is available in RandomFields, mev and SpatialExtremes, but only the latter offers conditional simulation of max-stable random fields (conditional on observed values at given locations) using Gibbs sampling [START_REF] Dombry | Conditional simulation of max-stable processes[END_REF]. CompRandFld's simulation routine for maxstable processes uses an interface to RandomFields. For a particular hierarchically structured max-stable dependence model, known as the Reich-Shaby model [START_REF] Reich | A hierarchical max-stable spatial model for extreme precipitation[END_REF]) that is constructed using spatial kernel functions and is derived from the spectral representation of a max-stable process based on a l p -norm [START_REF] Oesting | Equivalent representations of max-stable processes via p -norms[END_REF], estimation tools are available in the hkevp package and the experimental abba function in extRemes. It is difficult to fit because of the dual role of its nugget parameter α > 0. The hkevp package provides a Metropolis-within-Gibbs algorithm for Bayesian estimation of the model and for simulation, whereas the abba is flagged as experimental and is less comprehensive.

Peaks-over-threshold modeling

For functional peaks over threshold, mvPot (de Fondeville and Belzile, 2021) provides parametric simulation and estimation tools for various r -Pareto processes using Brown-Resnick and extremal Student-t dependence structures. Parameter estimates are computed using optimization of either full likelihood or gradient score functions; the latter remains computationally tractable for settings where full likelihood does not. Estimation of the marginal transformation T is not implemented and has to be performed prior to estimating the extremal dependence parameters using mvPot. A competitive estimation procedure is the gradient score estimating equation of de Fondeville and [START_REF] De Fondeville | High-dimensional peaks-over-threshold inference[END_REF], which does not require calculation of the normalizing constant of the model and also replaces censoring with downweighting. While statistically less efficient than full likelihood estimation, the procedure is more robust and can be applied in very high-dimensional settings. For estimation, numerical implementations are currently restricted to the Brown-Resnick model but tools for simulation and calculation of likelihoods are available for the extremal-t dependence model. The mev package also proposes likelihood functions and unconditional simulation routines for generalized r -Pareto processes (de Fondeville and Davison, 2021). Some other implementations allowing estimation of asymptotic dependence structures use original event data X i (s j ) and can be viewed as working on the interface of max-stable and peaks over threshold models. For example, moment-based estimation of parametric models, based on contrasting empirical and parametric versions of a variant of the so-called tail dependence function, is implemented in the package tailDepfun [START_REF] Einmahl | A continuous updating weighted least squares estimator of tail dependence in high dimensions[END_REF].

Modeling spatially varying marginal distributions

In practice, marginal distributions F s in functional data are usually not stationary, such that variation of marginal extreme value parameters with respect to space and time, or with respect to other available auxiliary variables, has to be captured. In the locationwise maxima setting, we can use use the generalized extreme value distribution and consider its parameters as functions of space, i.e., ξ(s), µ(s), σ(s). Different options exist in the peaks over threshold setting. A common approach is to fix a high, potentially nonstationary threshold u(s), and then estimate the threshold exceedance probability p(s) = Pr{X (s) > u(s)} = 1 -F s {u(s)} and the generalized Pareto parameters ξ(s), σ(s) based on observations of the exceedances X (s)u(s) > 0.

The regression frameworks discussed in Section 4.3 are relevant for modeling marginal extreme value parameters that vary with location in a first modeling step. Generalized additive modeling allows capturing complex nonlinear patterns of spatial nonstationarity using relatively large numbers of parameters. Some care may be required in tuning smoothing hyperparameters since in this step one usually assumes independence of observations X i (s j ), so functional dependence across space or time is disregarded. Specifically, MCMC-based Bayesian estimation of marginal parameters (using Gaussian process priors) for generalized extreme value distributions for maxima is possible through SpatialExtremes (Ribatet, 2022), and hkevp [START_REF] Sebille | hkevp: Spatial extreme value analysis with the hierarchical model of Reich and Shaby[END_REF] includes a similar function. The SpatialGEV package [START_REF] Chen | Fast approximate inference for spatial extreme value models[END_REF] provides a template for fitting latent spatial models with marginal generalized extreme value distributions and Gaussian process priors on the parameters using quadratic approximations to the marginal posterior. The unpublished package SpatGEVBMA, formerly on the CRAN, fits a latent model with generalized extreme value margins whose parameters follow Gaussian process priors with explanatory variables. Its defining functionalities are the use of Laplace approximation for automating proposals and efficient exploration of the posterior, and Bayesian model averaging to account for variable selection uncertainty [START_REF] Dyrrdal | Bayesian hierarchical modeling of extreme hourly precipitation in norway[END_REF].

An important alternative to Monte Carlo methods is to estimate complex integrals arising from Equation (6) through the integrated nested Laplace approximation (INLA). The INLA package proposes computationally convenient representations of the spatial Matérn covariance function through the stochastic partial differential equation approach of [START_REF] Lindgren | An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach[END_REF] for spatial and spatio-temporal latent Gaussian modeling. As mentioned in Section 2.4, INLA provides implementation for generalized extreme value distributions (with covariates and random effects in the location parameter) and the generalized Pareto distribution (with covariates and random effects in a quantile at a probability level α ∈ (0, 1) specified by the user; see [START_REF] Opitz | INLA goes extreme: Bayesian tail regression for the estimation of high spatio-temporal quantiles[END_REF] and Krainski et al. (2018, Chapter 6).

The package further allows joint estimation of several regression designs where some of the random effects can be in common (i.e., shared through a scaling coefficient) among these, which is beneficial to obtain cross-correlation in the posteriors of the predictors of several response types. For example, we could combine a logistic regression for the exceedance probability with a generalized Pareto regression for the excess above the threshold, and a shared random effect with a positive sharing coefficient would entail positive posterior correlation between the exceedance probability and the size of the excess.

Outlook for functional extremes

The coverage of max-stable processes, which remains an area of very active research, is much more comprehensive than others, with the notable exception of composite or full log likelihood inference for maxstable processes. Formulae exist for many partial derivatives of the exponent function V arising in the multivariate max-stable cumulative distribution functions and, in principle, the Stephenson-Tawn likelihood (or a bias-corrected version thereof) could be programmed for full likelihood inference beyond the bivariate case. Most of the models are also implemented with spatial applications in mind, even if temporal or spatio-temporal applications are possible. Max-infinitely divisible models are not covered in software yet, and Bayesian models with latent processes are often not provided with numerical implementation because of the complexity of implementation and also sometimes very long execution times of codes.

There are much fewer implementations for threshold models. Whereas their construction can be viewed as more flexible and intuitive than the one of the corresponding max-stable processes, they are conditional models with respect to threshold exceedance of the summary functional r . In the finitesample setting of statistical practice, this means that observations at some locations may not correspond to marginal exceedances and may therefore not be coherent with the asymptotic model. A common remedy is censoring, but this makes estimation more costly because the likelihood functions now include high-dimensional distribution functions which typically must be calculated via Monte-Carlo methods for each vector of observation. Generic full likelihood estimation procedures have been proposed, and are implemented (though computationally costly) for some models. However, available implementations do not yet come with a comprehensive set of available models and methods for parameter inference, model validation and comparison. An obvious solution to facilitate such implementation, provided that parameters are identifiable from lower-dimensional summaries, would be to use composite likelihood, but this is not implemented. Likewise, Bayesian generalized Pareto models with latent Gaussian process priors could be easily implemented in many probabilistic programming languages outside of R, such as Stan, but no general-purpose routines exist so far.

Simulation algorithms for unconditional simulation from generalized r -Pareto processes with arbitrary risk functionals r are still elusive, as designing efficient accept-reject methods requires case-by-case analysis. Available conditional simulation code typically amounts to simulation of elliptical distributions (log-Gaussian or Student-t ) with linear constraints.

Implementations with documented code are often available as supplementary material to methodological papers but have not been encapsulated in officially validated packages; see Huser and Wadsworth 9: Overview of simulation algorithms for bivariate (b) and multivariate (m) max-stable distributions and for max-stable processes (f), and for angular (a) and Pareto processes (p) with associated references. Some of the listed functions ( * ) are generic and include specific classes for max-stable models, but other models as well.

(2019); [START_REF] Bacro | Hierarchical space-time modeling of asymptotically independent exceedances with an application to precipitation data[END_REF]; [START_REF] Simpson | High-dimensional modeling of spatial and spatio-temporal conditional extremes using INLA and the SPDE approach[END_REF] for recent examples. [START_REF] Huser | Modeling spatial processes with unknown extremal dependence class[END_REF] has companion code for frequentist estimation of a flexible subasymptotic spatial model in the unpublished package spatialADAI. An INLA-based implementation for Bayesian conditional extremes models for spatial and spatio-temporal data is provided for [START_REF] Simpson | High-dimensional modeling of spatial and spatio-temporal conditional extremes using INLA and the SPDE approach[END_REF] as supplementary material. The implementation of many Bayesian extreme value models in the literature is achieved with standard MCMC algorithms that are tailored to the particular data application, but often generic and easily reusable or reproducible code is not provided, which hinders reproducibility.

Specialized topics

While our review has ranged mostly over software providing implementation of relatively generic methods that can be useful in various application contexts, there also has been active development of software libraries targeting specific application fields, and we here cite some of them. Hydrology and climate: Regional frequency analysis using L-moments is possible with the lmomrfa package. The climextRemes package leverages extRemes for climate extremes and implements methods highly relevant for this field, such as local likelihood fitting; the package is also available in Python.

IDF provides intensity-duration-frequency (IDF) curves [START_REF] Ulrich | Estimating IDF curves consistently over durations with spatial covariates[END_REF]. jointPm implements the method of [START_REF] Zheng | Efficient joint probability analysis of flood risk[END_REF] for evaluating bivariate probabilities of exceedance. An example of a highly specialized package is futureheatwaves and facilitates finding, characterizing and exploring heatwaves in climate projections, while the Python package teca is dedicated to tracking extremes of large scale climate models. Renext includes methods for peaks over threshold with a variety of distributions and the possibility to include historical maximum records, along with tests of exponentiality and goodness-of-fit. Financial and actuarial science: Some packages provide implementation of various generic models and methods for extreme values, but make strong use the semantics of those fields in their documentation. The packages QRM (and its successor qrmtools) and ReIns implement various functions to accompany the books [START_REF] Mcneil | Quantitative Risk Management: Concepts, Techniques and Tools -Revised Edition[END_REF] and [START_REF] Albrecher | Reinsurance: Actuarial and Statistical Aspects[END_REF], respectively. The package fExtremes [START_REF] Wuertz | fExtremes: Rmetrics -Modelling extreme events in finance[END_REF] provides functions for financial analysis used by the Rmetrics project.

The package VaRES [START_REF] Nadarajah | VaRES: Computes value at risk and expected shortfall for over 100 parametric distributions[END_REF] provides two popular risk measures (value at risk and expected shortfall) for a large collection of probability distributions, including many heavy-tailed distributions. The extremis package proposes functionalities to cluster multivariate financial time series based on their frequency and magnitude of extreme events. Machine learning: The interface between statistical machine learning and extreme values has been growing in recent years, with proposals encompassing the use of gradient boosting for extremes [START_REF] Velthoen | Gradient boosting for extreme quantile regression[END_REF], unpublished package gbex) and of extremal random forests [START_REF] Gnecco | Extremal random forests[END_REF], package erf) for modeling high quantiles of a univariate distribution. Another area of active research is open-set classification, dealing with classification of observations in categories not observed in training data: the Python package EVM implements the extreme value machine of [START_REF] Rudd | The extreme value machine[END_REF], whereas R package evtclass includes the algorithms described in [START_REF] Vignotto | Extreme value theory for anomaly detection -the GPD classifier[END_REF]. Survival analysis: Presence of censoring or truncation mechanisms, common in survival analysis, require dedicated software implementations because they affect the likelihood contribution of observations. The Matlab LATools [START_REF] Rootzén | Human life is unlimited -but short (with discussion)[END_REF] proposes an interface for interval-truncated generalized Pareto observations, while the longevity package handles more general partial observation schemes.

Discussion and conclusion

We have covered in this review a wide range of available software implementations for extremes, and we sincerely hope without omissions that are considered as important by authors or users. The development of extreme value software is key to extreme value analysis in practice and has become an active area of research, but the availability of implementations tends to lag behind methodological innovations since these are often not accompanied by generic, easily reusable and validated codes. To encourage modelers in applied sciences and in operational services to make use of the most advanced methods and models, off-the-shelf implementations are desirable. However, generic code may be difficult to provide due to the high sophistication of approaches as, for instance, with functional extremes. Designing generic estimation procedures that are flexible enough to be useful while at the same time being robust requires particular care. Writing this review made us aware of how challenging it is for the extreme value community to develop tested and easily reusable software that keeps pace with methodological progress: most software was written a decade ago, there are only handful of active maintainers, and most models proposed in the literature are not put together with software.

Many methods proposed in the last two decades are still not available; lack of implementation of bespoke code is a major impediment for their adoption. The most obvious gap is in software for fitting multivariate max-stable models (with composite likelihood) and multivariate generalized Pareto distributions with censoring in moderate dimensions for the parametric models with suitable tools. The conditional spatial extremes model, which extends the Heffernan-Tawn approach to the spatial setting, has been used in many recent papers but no software has been released.

More methods and more refined tools are also required for the nonstationary exploration and inference of extreme values. In many application fields, physical change processes (e.g., climate change, land-use change) require tools to explore, model and infer nonstationary behavior in extremes, for instance for climate-change detection and attribution. Currently, nonstationary modeling is implemented for marginal distributions through regression designs, but implementations providing dedicated methods for extreme value detection and attribution under climate change are scarce, and easily reusable codes for nonstationary extreme value dependence structures are yet missing.

Figure 1 :

 1 Figure 1: Left: profile log likelihood of η = -ξ/σ for a generalized Pareto distribution with scale σ and shape ξ.The black line shows one data set for which the conditional maximum ξ η = -1 lies on the boundary of the parameter space (black) and one with ξ η > -1 (grey). Right: conditional log likelihood surface for the inhomogeneous Poisson process at σ for simulated data (larger values have darker grey-scale shade). The white dot indicates the maximum likelihood estimate, while the hyperbola defines the feasible region of the parameter space given by the support constraints.

Fawcett

  abba

)p(θ)dθ ,

	package	functions	type		link	model
	ismev	gpd.fit, gev.fit linear	custom	GEV, GP
	evd	fgev	linear	linear	GEV (loc)
	texmex	emv	linear	log (scale)	-
	eva	gevrFit, gpdFit	linear	custom	GEVR, GP
	extRemes fevd	linear	linear/log (scale) -
	ismev	gamGPDfit	GAM		fixed	GP
	GEVcdn	gevcdn.fit	neural network fixed	GEV
	VGAM	gev, gp	GAM		custom	GEV, GP, *
	evgam	fevgam	GAM		fixed	GEVR, GP, *
	Table 1: Functionalities for generalized linear (linear) or generalized additive modelling (GAM) of the pa-
	rameters of extreme value distributions. Model families supported include generalized extreme value distri-
	bution (GEV), generalized Pareto (GP), r -largest extremes (GEVR) and more general families or special cases
	of extreme value distributions, miscellaneous (*).
	2.4. Bayesian modeling
	2.4.1 Generalities
	In the Bayesian paradigm, the likelihood of the data vector Y is combined with prior distributions for
	the model parameters θ = (θ 1 , . . . , θ m ) ∈ Θ, with prior density p(θ); we use the generic notation p(. . .)
	for various conditional and unconditional densities. The distribution of the data given the parameter is
	obtained from the likelihood function p(y | θ) ∝ exp{ (θ; y)}, where ∝ denotes proportionality.
	The posterior distribution,
		p(θ | y) =	p(y | θ)p(θ) p(y | θ

Table 2 :

 2 Comparison of R packages for Bayesian univariate extreme value modeling. Families: generalized extreme value distribution (1), generalized Pareto distribution (2), inhomogeneous Poisson process (3), order statistics/r -largest (4) or custom/other (*). Sampling: random walk Metropolis-Hastings (RWMH), exact sampling ratio-of-uniforms (RU), independent Metropolis-Hastings (IMH); the INLA package uses deterministic Laplace approximations. "PC" priors refer to penalized complexity priors. All packages, except

	package	function	models covariates	sampling prior choice
	evdbayes	posterior 1-4	loc./thresh RWMH	multiple
	extRemes	fevd	1-4,*	all	RWMH	custom
	INLA	inla	1-2,*	loc./thresh -	PC
	MCMC4Extremes ggev, . . .	1-2,*	no	RWMH	fixed
	revdbayes	rpost	1-4	no	RU	custom
	texmex	evm	1-2,*	all	IMH	Gaussian

evdbayes, also provide S3 methods (notably plot and summary). All packages return a matrix of posterior draws.

Table 3 :

 3 but no software implementation exists at the time of writing. Main functionalities of R packages for nonparametric Hill-type estimators of the shape parameter,

	package	estimation	function	features
	evir	-	hill	e,p
	evmix	smoothing	hillplot	p
	evt0	location invariant	gh, PORT.Hill	p, q
	extremefit	weighted, time series	hill, hill.adapt, hill.ts	e, p, q, o
	ExtremeRisks time series, CI	HTailIndex, EBTailIndex	e, o
	fExtremes	-	hillPlot, shaparmHill	e, p
	ptsuite	-	alpha_hills	e
	QRM	-	hill, hillPlot	e, p
	rbm	random block	rbm, rbm.plot	e, p
	ReIns	conditional, censoring (c)Hill, (c)genHill, crHill, . . . e, p
	tea	smoothing	althill, avhill	p
	including functionalities for estimation of the shape or tail index (e), Hill threshold diagnostic plots (p),
	quantile estimates (q) and other methods (o).

2.6. Threshold selection

Many methods are driven by analyses of the most extreme observations. In the univariate case, these are the k largest order statistics or, equivalently, observations that exceed a threshold u as presented MLE -Table 4: Comparison of R packages for classical visual methods. Stability: function name for a threshold stability plot; models: either generalized Pareto (1), inhomogeneous Poisson process (2) or extended generalized Pareto model of Papastathopoulos and Tawn (2013) (3); inference: method of inference, either maximum likelihood estimation (MLE) or Bayesian (B); MRL: mean residual life plot, if applicable.

Table 5 :

 5 Overview of formal threshold selection methods and numerical implementations

	type	methods	package	function
	penultimate	Northrop and Coleman (2014)	mev	NC.diag
		Papastathopoulos and Tawn (2013) mev	tstab.egp
	goodness-of-fit Gerstengarbe and Werner (1989)	tea	ggplot
		Hosking and Wallis (1997)	POT	lmomplot
		Bader et al. (2018)	eva	gpdSeqTests
	sequential	Wadsworth (2016)	mev	W.diag
		Thompson et al. (2009)	tea	TH
		del Castillo and Padilla (2016)	ercv	cvplot, thrselect
	predictive	Northrop et al. (2017)	threshr	ithresh
	mixture	Scarrott and MacDonald (2012)	evmix	-
		Durrieu et al. (2015)	extremefit •paretomix
		Naveau et al. (2016)	mev	fit.extgp

Table 6 :

 6 Overview of some direct estimators of the extremal index with associated references and tuning parameters.

	package	estimator(s) estimation	UQ diagnostics
	evd	runs, FS	exi	no exiplot
	evir	blocks 2	exindex	no exindex
	extRemes runs, FS	extremalindex	yes -
	exdex	ILS	iwls	no -
		K -gaps	kgaps	yes choose_uk
		SPM	spm	yes choose_b
	fExtremes runs	runTheta	no exindexPlot
		blocks 1	clusterTheta	no exindexesPlot
		blocks 2	blocktheta	no
		intervals	ferrosegersTheta no
	mev	ILS, FS	ext.index	no ext.index
		K -gaps	ext.index	no infomat.test, ext.index
	POT	runs	fitexi	no exiplot
	revdbayes K -gaps	kgaps_post	yes -
	texmex	FS	extremalIndex	yes extremalIndexRangeFit
	tsxtreme	runs	thetaruns	yes -

Table 7 :

 7 Comparison of R packages for the direct estimation of the extremal index.

	Estimator(s): name(s) of

preprint, version of May 17, 2022Author's Original Version

Author's Original Version preprint, version of May 17, 2022

Acknowledgements

Funding in partial support of this work was provided by the Natural Sciences and Engineering Research Council (RGPIN-2022-05001, DGECR-2022-00461).