Quentin Sebille

Anne-Laure Fougères

Cécile Mercadier

A comparison of spatial extreme value models. Application to precipitation data

Keywords: Spatial modelling of extreme events, extreme value theory, max-stable processes, hierarchical models, spatial prediction, precipitation data. 2010 MSC: 60G70, 62G32, 62H11, 62P12

In this paper, focus is done on spatial models for extreme events and on their respective efficiency regarding the estimation of two risk measures: one extrapolating marginal distributions and one summarizing the spatial bivariate dependence of extremes. A wide comparison is performed on a simulation plan that has been specifically designed from a daily precipitation data set. The objective of this paper is twofold: firstly, pointing out the inherent properties of each model, and secondly, advising users on how to choose the model depending on the specific type of risk.

Introduction

Analyses of extreme values of environmental variables such as precipitation are of great importance since it involves human lives as well as considerable losses of money when a catastrophic event occurs. Recently for instance, on the 3rd October 2015, the region of Cannes (France) faced such an extreme event with a huge and sudden amount of rainfall that caused twenty deaths and millions of Euros in damages. Accurate risk measures for such extreme phenomena are therefore needed to prevent from this type of scenario.

The risk estimation of these events is challenging because they involve values that are beyond the range of the observations. For this purpose, adapted tools come from extreme value theory. See for instance de [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF], [START_REF] Beirlant | [END_REF], [START_REF] Finkenstädt | Extreme values in finance, telecommunications and the environment[END_REF] and [START_REF] Coles | An introduction to statistical modeling of extreme values[END_REF]. Since precipitation phenomena have a spatial feature and data are generally observed at several stations, the dedicated setting to handle this question is that of max-stable spatial models. Detailed and helpful reviews on these models are [START_REF] Cooley | A survey of spatial extremes: Measuring spatial dependence and modeling spatial effects[END_REF], Davison et al. (2012) or [START_REF] Ribatet | Spatial extremes: Max-stable processes at work[END_REF].

Six max-stable spatial models have been selected among the most popular or the most recent in the literature. This choice includes Bayesian and frequentist concepts, and goes from simple to smooth spatial dependence structure. Within this paper, the main goal is to answer: Which of the six competing models yields the best prediction for extreme behaviors of simulated processes mimicking precipitation data? Addressing this question supposes in particular a careful device of the data sets involved, as well as a relevant choice of performance criteria.

Two comparative criteria adapted to extreme events prediction are evaluated. The estimation of rare events at a location where no data is available is handled first; this induces the capacity of spatial extrapolation of the extreme behavior when looking at marginal information only. A clear and well known way to summarize this marginal information into a concrete risk measure is the return level. Then a second and complementary criterion is the measure of extreme sets for the bivariate distribution at a pair of locations. This aims at capturing the spatial dependence structure of extremes. One could also study higher dimensional indices but most of the models have tractable formulae only in two dimensions.

Several options can be chosen to define the terms of comparison. One possible option could be to start from an expert point of view and compare each model with an a priori value of the previous criteria. To depart from a subjective choice, an intensive simulation study has been preferred. A parametric bootstrap procedure is used to produce simulations as close as possible to a true phenomenon. More precisely, a real precipitation data set is considered over France, on which each of the six maxstable spatial models is fitted. These fitted models are then fixed to play the role of extreme rainfall generators. The two comparative criteria are finally evaluated on each generated sample and faced to the corresponding (known) true value.

The remainder of this article is organized as follows. The theoretical background of extreme value theory is addressed in the following section, with an emphasis on the so-called block maxima approach and on the six max-stable spatial models of interest. Section 3 describes the designed simulation plan, the two criteria of competition and the results. Conclusions are drawn in Section 4, where some recommendations are provided with respect to different purposes.

Notation and models

Definition of max-stable processes

Let S be a compact subset of R d that represents the spatial region of interest, d being a positive integer. Consider a random process

Y (•) = {Y (s)} s∈S defined over S, with continuous sample paths. Write Y 1 (•), . . . , Y T (•) for independent copies of Y (•). The process Y is called max-stable if for each T > 1, there exist continuous functions a T (•) > 0 and b T (•) ∈ R such that: T t=1 Y t (•) -b T (•) a T (•) d = Y (•) ,
where d = denotes equal in distribution. Such max-stable processes arise as non degenerate limits for pointwise maxima of stochastic processes on S, and have been introduced by de Haan (1984). In particular, for each s ∈ S, the random variable Y (s) follows a generalized extreme value distribution GEV(µ(s), σ(s), ξ(s)). Location µ(s) ∈ R, scale σ(s) > 0 and shape ξ(s) ∈ R parameters are indexed by s. Recall that the cumulative distribution function (cdf) of a GEV(µ, σ, ξ) random variable Y is:

P(Y y) = exp -1 + ξ y-µ σ -1/ξ + if ξ = 0 exp -exp -y-µ σ if ξ = 0 ,
where z + denotes the positive part of z. Note that the usual lack of clear spatial pattern for the shape ξ(•) when dealing with precipitation data, jointly with the difficulty of estimating this parameter lead in this paper to consider ξ(•) ≡ ξ 0 . Thanks to the one-to-one mapping:

Z(s) = 1 + ξ 0 Y (s) -µ(s) σ(s) 1/ξ0 + , (1)
one obtains a simple max-stable process Z(•), that is with unit Fréchet margins, corresponding to GEV(1,1,1). The joint cdf of Z(•) at a set of sites {s 1 , . . . , s n } ⊂ S is given by:

P (Z(s 1 ) z 1 , . . . , Z(s n ) z n ) = exp [-V (z 1 , ..., z n )] ,
in terms of an exponent measure V that contains the information about the spatial dependence of extremes (see e.g. (de Haan and Ferreira, 2006, Chapter 9)). The two limit cases are the independence case, with V (z 1 , . . . , z n ) = n i=1 z -1 i and the total positive dependence case, with V (z 1 , . . . , z n ) = n i=1 z -1 i . A common measure of spatial dependence between the set of sites {s 1 , ..., s n } is the extremal coefficient θ of [START_REF] Schlather | A dependence measure for multivariate and spatial extreme values: Properties and inference[END_REF] defined by:

P (Z(s 1 ) z, ..., Z(s n ) z) = P (Z(s 1 ) z) θ .
One can interpret the coefficient θ ∈ [1, n] as the number of components of {Z(s 1 ), ..., Z(s n )} that are independent and it is linked to the exponent measure via the relation θ = V (1, ..., 1). Note that θ = θ(s 1 , ..., s n ) depends on the set of sites {s 1 , . . . , s n }, as does the exponent measure V . Finally, the description of a max-stable process is done in two parts. First, the marginal effect is captured by the processes that represent the GEV parameters. Classically, the inference is done under a linear model involving covariates. Such details are postponed until later (as in Equation ( 6) for instance). Second, the spatial dependence of extremes is measured by V (or summarized through θ). Parametric models for V can be helpful and some examples are presented in the next section.

Spectral representation and parametric models

Max-stable processes can be described thanks to the following spectral representation, due to de [START_REF] De Haan | A spectral representation for max-stable processes[END_REF]. Let {ζ j } j∈N be a Poisson point process on (0, ∞) with intensity measure dζ/ζ 2 and consider independent copies {W j (s), s ∈ S} j∈N of a stationary process W (•) verifying sup s∈S W (s) < ∞ and E [W (s) + ] = 1. Then the process Z(•), defined for each s ∈ S by:

Z(s) = max j 1 ζ j W j (s) , (2)
is max-stable with unit Fréchet margins. Its joint cdf is expressed as:

P (Z(s) z(s), s ∈ S) = exp -E sup s∈S W (s) z(s) .
Different choices for the so-called spectral processes W j (•) in (2) lead to different max-stable models. Four of them are now presented.

The Smith model: GEVP

Consider W j (s) = f (s-v j )
, where f is the d-variate normal density with zero mean and covariance matrix Σ, and {v j } j∈N is an homogeneous Poisson point process on R d . In this construction, Z(s) can be seen as the maximum over an infinite number of storms located at {v j } j∈N with severity {ζ j } j∈N . The effect of these storms at a given point s is then described by ζ j f (s -v j ). This particular form for f has been considered by [START_REF] Smith | Max-stable processes and spatial extremes[END_REF] and leads to one of the most classical max-stable model. It is called Gaussian Extreme Value Process (GEVP) in this paper.

The exponent measure can be written for a pair of sites (s 1 , s 2 ) as:

V GEVP (z 1 , z 2 ) = 1 z 1 Φ a(h) 2 + log(z 2 /z 1 ) a(h) + 1 z 2 Φ a(h) 2 + log(z 1 /z 2 ) a(h) , (3)
in terms of h = s 1 -s 2 , where a 2 (h) = h T Σ -1 h is the Mahalanobis distance between the sites s 1 and s 2 and where Φ(•) is the cdf of the standard normal distribution. The extremal coefficient of this model is given by θ(h) = 2Φ {a(h)/2}. The dependence parameters that need to be estimated for this model are the elements of the variance-covariance matrix Σ: precisely Σ 11 , Σ 12 and Σ 22 when d = 2.

The Schlather model: EGP

A second choice is to take W j (•) as stationary Gaussian processes, with correlation function ρ(•), scaled so that E [W (s) + ] = 1. The resulting max-stable process, defined by [START_REF] Schlather | Models for stationary max-stable random fields[END_REF], is called the Extremal Gaussian Process (EGP). The corresponding exponent measure is:

V EGP (z 1 , z 2 ) = 1 2 1 z 1 + 1 z 2 1 + 1 -2 [ρ(h) + 1]z 1 z 2 (z 1 + z 2 ) 2 , ( 4
)
where h is now the Euclidean distance between the sites s 1 and s 2 . The extremal coefficient of the EGP is θ

(h) = 1 + {1 -ρ(h)}/2.
A shortcoming of the EGP is that the bivariate extremal coefficient does not span over the interval [1,2] as it should do, but over [1, 1.838] instead since ρ(•) is positive definite. In other words, this model never allows independence of extremes, even when the distance between two sites increases indefinitely.

Different correlation functions ρ(•) can be chosen. We work with the powered exponential form

ρ(h) = (1 -η) exp [-(h/λ) ν ]
, where η ∈ [0, 1), ν ∈ (0, 2] and λ > 0 are respectively the nugget, the smooth and the range parameters. They need to be estimated when fitting the EGP. We thus denote them η EGP , ν EGP and λ EGP .

The Brown-Resnick model: BRP

A third possibility is to take,

W (•) = exp [ε(•) -γ(•)],
in terms of a Gaussian process ε with stationary increments and semivariogram γ. Then, the representation (2) leads to the so-called geometric Gaussian process. [START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF] showed that choosing ε as a fractional Brownian motion yields the process introduced by [START_REF] Brown | Extreme values of independent stochastic processes[END_REF], called the Brown-Resnick Process (BRP). Its exponent measure has a similar form to the GEVP one:

V BRP (z 1 , z 2 ) = 1 z 1 Φ a(h) 2 + log(z 2 /z 1 ) a(h) + 1 z 2 Φ a(h) 2 + log(z 1 /z 2 ) a(h) , (5)
where a 2 (h) = 2γ(h) and h is the Euclidean distance between the sites s 1 and s 2 . Its extremal coefficient is therefore given by θ(h) = 2Φ{ γ(h)/2}. When ε(•) is a fractional Brownian motion, the semivariogram γ(•) is given by γ(h) = h λ ν , where λ and ν are respectively the range and nugget parameters. These parameters must be estimated to describe the spatial structure, and they are respectively denoted by λ BRP and ν BRP .

2.2.4. The t-extremal model: TEP Opitz (2013) obtains the only possible max-stable limit for asymptotically dependent elliptical processes, namely the t-extremal process, denoted TEP here. Let δ > 0 and {ζ i } be a Poisson point process on (0, ∞) with intensity dΛ(t) = δt -(δ+1) dt. Let the W i (•)'s be independent copies of a stationary standard Gaussian process with correlation function ρ

(•). Set m δ = 2 (δ-2)/δ √ π -1 Γ δ+1 2
. The TEP has then the following spectral representation:

Z(•) = m -1/δ δ max i 1 ζ i W i (•) ,
and the associated exponent measure

V TEP (z 1 , z 2 ) = 1 z 1 T δ+1 z 2 z 1 1/δ ; ρ(h) , 1 -ρ(h) 2 δ + 1 + 1 z 2 T δ+1 z 1 z 2 1/δ ; ρ(h) , 1 -ρ(h) 2 δ + 1 ,
where T ν (•; µ, σ) is the cdf of the non-standard Student distribution with ν degrees of freedom, ρ is the powered exponential correlation function and h is the Euclidean distance between s 1 and s 2 . Note that both the EGP of [START_REF] Schlather | Models for stationary max-stable random fields[END_REF] and the BRP of [START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF] appear to be special cases of the TEP model. Indeed, the EGP is obtained straightforwardly when δ = 1, while the TEP model converges towards the BRP when δ → ∞. The dependence parameters that need to be estimated are then: δ TEP , the degree of freedom, and the parameters of the correlation function, namely η TEP , λ TEP and ν TEP .

Hierarchical modelling

Hierarchical models suppose that the response variable, say Y (•) in our case, is spatially independent conditionally on unobserved latent processes or variables. The interested reader may find examples in [START_REF] Banerjee | Hierarchical modeling and analysis for spatial data[END_REF] for instance. Models based on hierarchical approaches are usually defined within the Bayesian paradigm. The estimation is generally performed through a Markov Chain Monte Carlo (MCMC) algorithm which handles easily the conditional densities integration. Recent studies exploiting hierarchical modelling for extreme precipitation are for instance [START_REF] Apputhurai | Spatiotemporal hierarchical modelling of extreme precipitation in western australia using anisotropic gaussian random fields[END_REF] and [START_REF] Dyrrdal | Bayesian hierarchical modeling of extreme hourly precipitation in norway[END_REF].

The latent variable model: LVM

Davison et al. ( 2012) introduces a simple hierarchical structure for spatial extremes, called the latent variable model (LVM). The response variables {Y (s)} s∈S are assumed to be independent conditionally on latent processes that describe the GEV parameters; more precisely

Y (s)|{µ, σ, ξ} indep ∼ GEV(µ(s), σ(s), ξ(s)) , µ(s) = β T µ c(s) + ε µ (s) , σ(s) = β T σ c(s) + ε σ (s) , ξ(s) ≡ ξ 0 .
The

ε• (h) = δ • exp(-h/λ • ).

The Reich and Shaby model: HKEVP

From Davison et al. (2012), one knows that the LVM is particularly appealing when the estimation of the marginal distributions is of interest, as it focuses on modelling the GEV parameters. The main drawback is that the dependence structure of extremes is not considered since spatial independence is assumed.

The aim of this subsection is to present a model that describes both the marginal effect and the dependence structure within a Bayesian framework. The Hierarchical Kernel Extreme Value Process (HKEVP) has been introduced by Reich and Shaby (2012) and further developed in Shaby and Reich (2012) and [START_REF] Reich | A hierarchical model for serially-dependent A study of heat waves in the western us[END_REF]. It is defined as follows. Suppose that Y (s) ∼ GEV(µ(s), σ(s), ξ(s)) and model the margins by:

   µ(s) = β T µ c(s) + ε µ (s) , log[σ(s)] = β T σ c(s) + ε σ (s) , ξ(s) ≡ ξ 0 , (6)
The Whittle-Matérn form is taken rather than exponential one for the correlation function of ε µ (•) and ε σ (•), that is:

ρ(h) = 2 1-ν Γ(ν) h λ ν K ν h λ
, where ν > 0 and λ > 0 are respectively the smooth and range parameters, and Γ(•) and K ν are the Gamma and modified Bessel functions.

Consider now Z(•) the associated simple max-stable process. Assume that Z(•) = U (•)ϑ(•) where:

• U (•) is a spatially-independent process with common marginal distribution GEV(1, α, α). The parameter α ∈ (0, 1] is unknown.

• ϑ(s) = L =1 A ω 1/α (s) α
, with the same α as in U , describes the spatial dependence structure driven by some kernels {ω (•)} =1,...,L and some random effects {A } =1,...,L . These last variables are assumed to be independent and identically distributed from PS(α), the positive stable distribution with characteristic exponent α. In the definition of ϑ, {ω (•)} =1,...,L stands for a family of deterministic kernel functions normalized so that L =1 ω (s) = 1, ∀s ∈ S. These kernels have been chosen according to Reich and Shaby (2012)'s choice, namely as rescaled Gaussian densities (see [START_REF] Smith | Max-stable processes and spatial extremes[END_REF]). More precisely, let V := {v 1 , ..., v L } ⊂ S be a set of knots, and

ω (s) = K(s|v , τ ) L j=1 K(s|v j , τ )
, where

K(s|v , τ ) = 1 2πτ 2 exp - 1 2τ 2 (s -v ) T (s -v ) ,
for τ a bandwidth parameter that needs to be estimated, along with the dependence parameter α. The hierarchical formulation for the process Y (•) follows:

Y (s)|µ, σ, ξ, α, A, V indep ∼ GEV(µ * (s), σ * (s), ξ * 0 ) , µ * (s) = µ(s) + σ(s) ξ(s) ϑ(s) ξ(s) -1 , σ * (s) = ασ(s)ϑ(s) ξ(s) , ξ * 0 = αξ 0 .
The exponent measure of the HKEVP is given by:

V HKEVP (z 1 , . . . , z n ) = L =1 n i=1 ω (s i ) z i 1/α α .
To illustrate the central role of α, consider its boundary values. If α tends to 0, the so-called nugget process U (•) will be approximately equal to 1, so that Z(•) ≈ ϑ(•). The max-stable process is then spatially smooth and the spatial dependence strong. Now, if the parameter α is close to 1, the positive stable distribution PS(α) tends to be the Dirac at the singleton {1}. Combined with the normalized condition on the kernel functions, one gets ϑ(•) ≡ 1. The process Z(•) = U (•) will be spatially independent, which is then similar to the LVM. Note also that if, in addition, the set of knots V is equally-spaced and its cardinal L tends to infinity, then one obtains the definition of the GEVP.

Comparison of the spatial models

This section is threefold. Firstly, the design of the simulation study is settled. Secondly, the maxstable models are fitted on the resulting simulated data, and results are analyzed. Finally, additionnal characteristics of these models are discussed.

Designing the simulation as a precipitation data set

In this section, the six models described in Section 2, namely the LVM, HKEVP, GEVP, EGP, BRP, and TEP are applied to a set of precipitation data recorded in France. Along the paper, the models are sorted this way to respect an increasing "smoothness within dependence modelling", going from conditional independence to a spatial dependence structure with finite conditioning and then to four continuous max-stable dependence structures.

The precipitation data

This data set is extracted from the European Climate Assessment & Dataset (ECA&D) website (http://eca.knmi.nl/). The thirty-nine stations we consider1 are viewed with corresponding elevation on Figure 1. At each site, focus is done on the period 1960-1999. The yearly maximum of daily cumulative precipitation is then computed for each of these forty years. If there is more than 30 days of missing values in a year a missing value is allocated. The resulting data set contains 1.47% of missing values.

-4 -2 0 2 4 6 8 42 44 46 48 50 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0 250 500

Fitting procedures

The collected yearly maxima are considered as the realizations of some spatial max-stable process Y (•). The six spatial models can thus be fitted through the freely available software R2 . Most of the fitting functions come from the SpatialExtremes package [START_REF] Ribatet | SpatialExtremes: Modelling Spatial Extremes[END_REF]). More precisely, its function fitmaxstab (resp. rmaxstab) allows to fit (resp. simulate) from several max-stable models, as the GEVP of [START_REF] Smith | Max-stable processes and spatial extremes[END_REF], the EGP of [START_REF] Schlather | Models for stationary max-stable random fields[END_REF], the BRP of [START_REF] Brown | Extreme values of independent stochastic processes[END_REF] and the TEP of [START_REF] Opitz | Extremal t processes: Elliptical domain of attraction and a spectral representation[END_REF], and its function latent can be used to fit spatial extremes with the LVM of Davison et al. (2012).

As far as we know, there is no package including the fitting procedure for the HKEVP of [START_REF] Reich | A hierarchical max-stable spatial model for extreme precipitation[END_REF]. The only connected function is the routine abba included in the package extRemes of [START_REF] Gilleland | extRemes: New software to analyze how extremes change over time[END_REF] and related to the recent paper of [START_REF] Stephenson | Estimating spatially varying severity thresholds of a forest fire danger rating system using max-stable extreme-event modeling[END_REF]. Therein, the model uses a CAR prior over a network of thousands of gridded locations and is therefore not designed to make prediction outside the observed set of sites. The authors of the HKEVP propose nonetheless an open code available on Reich's website3 . Sebille (2016) recently published on CRAN the R package hkevp that contains in particular a routine (hkevp.fit) fitting the HKEVP. This function is widely inspired by Reich & Shaby's code, and the main changes are listed in the reference manual of hkevp.

Fitting the models on the ECA&D precipitation data set

The fitting functions introduced previously are now applied to the data set of yearly maxima of daily precipitation described in Section 3.1.1. The marginal GEV parameters are estimated at each of the thirty-nine meteorological stations. Available covariates are used in their estimation. More precisely, the following linear forms are used through the routine fitmaxstab for the GEVP, EGP, BRP and TEP fits:

   µ(s) = β µ,0 + lon(s)β µ,1 + lat(s)β µ,2 + alt(s)β µ,3 + msp(s)β µ,4 , σ(s) = β σ,0 + lon(s)β σ,1 + lat(s)β σ,2 + alt(s)β σ,3 + msp(s)β σ,4 , ξ(s) = β ξ,0 , (7)
for all s ∈ S where lon(s), lat(s), alt(s) and msp(s) stands respectively for the longitude, the latitude, the altitude and the mean seasonal precipitation (msp) associated to the position s. The choice of the msp as spatial covariate is motivated by [START_REF] Cooley | Bayesian spatial modeling of extreme precipitation return levels[END_REF]. The spatial form ( 7) is consistent with the one used in Davison et al. (2012), though it adds the altitude and msp covariates in the location and scale parameters. For the LVM and the HKEVP, the estimation of marginal parameters is included in the routines. The same spatial covariates in ( 7) are used to model the mean function of the latent processes.

Estimations of µ(•), σ(•) and ξ(•) are very similar from one model to another. Their values are not presented here for seek of concision. For a visual information, the location and scale parameters obtained using the BRP are given on Figure 2. To complete these patterns, note that the shape parameter varies from -0.04 to 0.27.

-4 -2 0 2 4 6 8 42 44 46 48 50 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 24 28 32 36 40 44 48 µ -4 -2 0 2 4 6 8 42 44 46 48 50 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 5 7 9 11 13 15 17 19 σ to independence. For the GEVP, the covariance matrix is almost degenerate. In the EGP case, the nugget η EGP is near one, which results in a process close to pure noise. For the BRP model, the range parameter λ BRP is lower than 0.1622, the minimum distance between observed positions. For the TEP, the high value of δ TEP is a sign for independence since the bivariate extremal coefficient θ tends to 2 as δ tends to infinity. Finally, the HKEVP reveals a structure close to the conditional independence with α near 0.94 . Note that standard deviations are not displayed for the GEVP and the EGP since maximization procedures of the log-likelihood functions did not converged. It can also be remarked that standard deviations for δ TEP and λ TEP are abnormally large. This is actually observed on more general fits: estimation of these two dependence parameters with the TEP on simulations from simple max-stable processes are highly uncertain, with very large standard deviations. The near-independence pattern observed in Table 1 is consistent with the fact that the meteorological network used in this article is rather sparse over the region of interest (see Figure 1). Extreme precipitation data are known to show a low range of spatial dependence. For instance, [START_REF] Fawcett | Estimating the probability of simultaneous rainfall extremes within a region: a spatial approach[END_REF] found out that the dependence between extreme precipitation in Great Britain is genuine for a distance h < 100km. A brief analysis on the data set used in the present paper shows that there is roughly 6% of all the inter-distances between the meteorological stations that are under this limit case. Our approach remains nevertheless justifiable since the case of independence is allowed by most of the spatial models considered in this paper. The only exception is for the EGP, though this disadvantage could be solved by using the ideas of [START_REF] Davison | Geostatistics of extremes[END_REF]. In their max-stable model, the Gaussian spectral processes {W j (•)} j∈N are defined over compact random sets {B j } j∈N and equal 0 elsewhere. In this case, the independence is obtained when the distance between two sites becomes greater than the volume of the B j 's.

Comparative study of the spatial models

The six spatial models will compete through the simulation design presented below.

Simulated data sets

The aim of this section is to present the way to produce simulations of annual maxima of daily precipitation via a parametric bootstrap procedure. For this purpose, samples are drawn from the six models with parameters fixed from Section 3.1.3. Routines rmaxstab and hkevp.rand are used from packages SpatialExtremes and hkevp respectively. We denote by:

Y (k) M S :=     Y (k) 1 (s 1 |M S ) . . . Y (k) 1 (s n |M S ) . . . . . . . . . Y (k) ny (s 1 |M S ) . . . Y (k) ny (s n |M S )     ny×n ,
the k-th replicate, for k ∈ {1, ..., K = 50}, coming from the simulation model M S . Data are generated on the n = 39 positions {s 1 , ..., s n } coinciding with the sites of Figure 1. The number of simulations n y = 40 is the maximal length observed of yearly maxima series per station on this data set.

Since the generated data come from the fitted models, the marginal and dependence parameters used are the ones estimated on real data. This procedure aims first at producing simulations as close as possible to real annual maxima of precipitation and second to depart from one particular simulation model M S that may favor one of the six spatial models put in competition.

However, numerical issues are encountered when simulating from the GEVP and the BRP:

-the non-definite form of the covariance matrix Σ prevents from simulating data with the GEVP, -the low value of the range λ BRP generates a bug inside the simulation function rmaxstab: the margins of the obtained process Z(•) are no longer unit Fréchet.

For these reasons, simulations from these two models are not taken into account in the next sections. Therefore, we have M S ∈ {LVM, HKEVP, EGP, TEP}.

Two comparative criteria

The performance of the competing spatial models is measured on the two criteria described below.

The first criterion focuses on the quality of spatial extrapolation of the GEV marginal distribution. It consists in the estimation of the T -year return level at a "target" site s * where no data is available. For any period of time T , this value is a quantile (denoted by y T (s * )) of the marginal distribution that we expect to be exceeded once over T years. It is therefore defined by P(Y (s * ) y T ) = 1 -1/T . Knowing the marginal parameters µ, σ and ξ evaluated at s * , it is possible to explicitly compute its value:

y T (s * ) = µ(s * ) + σ(s * ) ξ(s * ) log T T -1 -ξ(s * ) -1 . ( 8
)
The second criterion provides a measure of the spatial dependence structure of a specific model. It is simply the estimation of a joint probability at a pair of sites (s 1 , s 2 ):

p := P {Y (s 1 ) y T (s 1 ), Y (s 2 ) y T (s 2 )} ,
where y T (s 1 ) and y T (s 2 ) are respectively the T -year return levels of precipitation at sites s 1 and s 2 . It is easy to see that this is equivalent to computing the exponent measure V (z 1 , z 2 ) evaluated at (s 1 , s 2 ) with

z 1 = z 2 = log T T -1 -1 .
For this criterion, the pair (s 1 , s 2 ) is chosen as the two closest sites of {s 1 , ..., s n } in order to capture some spatial correlation, because of the low range of dependence observed on precipitation data. This pair is referred by triangles in Figure 1. With T = 100, we have p ∈ [0.9801, 0.99], spanning the different cases between independence (for p = 0.9801) and complete dependence (for p = 0.99). Since we decided to generate n y = 40 years of data per simulation set, the same as in the real data set, it is reasonable to apply these criteria for the fixed value T = 100. and L 1 -distances

1 K K k=1 y 100 (k) (s i , M S |M F ) -y 100 (s i , M S ) i=1,...,n
have also been drawn, but they are not provided since they showed similar results to Figures 3 and4. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Bias LVM HKEVP GEVP EGP BRP TEP -100 -50 0 50 100 Simulations from LVM q q q q q q q q q q q Bias LVM HKEVP GEVP EGP BRP TEP -100 -50 0 50 100 Simulations from HKEVP q q q q q q q q q q q q q q q q q q Bias LVM HKEVP GEVP EGP BRP TEP -100 -50 0 50 100 Simulations from EGP q q q q q q q q q q q Bias LVM HKEVP In Figure 3, 4 and 5, the model M S used to produce the simulated data is indicated at the top of each of the four plots. The results are then displayed by fitting models M F , whose name are indicated at the bottom of each boxplot. The dashed red line corresponds to 0 in Figure 3 and 4 and to the exact value p(M S ) in Figure 5. From these three figures, it is possible to order roughly the six spatial models of extreme values by their performances over the two criteria defined in Section 3.2.2.

The "best" one may be the TEP of [START_REF] Opitz | Extremal t processes: Elliptical domain of attraction and a spectral representation[END_REF]. It shows indeed great precision toward the extrapolation of the 100-years return level and has a good robustness on the simulation model M S . q q q q q q q q q q q q q q q q q q q q q RMSE LVM HKEVP Simulations from EGP q q q q q q q q RMSE LVM HKEVP Moreover, in terms of spatial dependence structure, this model is the most flexible and the most precise to estimate the joint probability p.

The LVM of Davison et al. (2012) outperform the five others if the goal is the extrapolation of the marginal effect, represented in this comparative study by the 100-years return level at a target site. This is not surprising since it is the only one whose construction focuses only on the spatial modelling of the marginal parameters. However, this model fails down if the goal of the user is to estimate p or any measure linked to the dependence structure, because of the assumption that the whole spatial dependence is explained by the margins. It can thus be observed in Figure 5 that the joint probability is always estimated by p = 0.9801 = 0.99 2 which means total independence. The BRP of [START_REF] Brown | Extreme values of independent stochastic processes[END_REF] and the HKEVP of [START_REF] Reich | A hierarchical max-stable spatial model for extreme precipitation[END_REF] may be classified as quite good compromises for the two criteria. The former tends to underestimate the 100-years return level, in particular if M S is the LVM, while the latter seems to overestimate this criterion, in particular if M S is the EGP. Since this criterion is a well-used risk measure, the safest choice for this q q q q q q q q q q q q q q q q q q q q q q q q p LVM HKEVP GEVP EGP BRP TEP 0.980 0.982 0.984 0.986 0.988 0.990

Simulations from LVM q q q q q q q q q q q q p LVM HKEVP GEVP EGP BRP TEP 0.980 0.982 0.984 0.986 0.988 0.990

Simulations from HKEVP q q q q q q q q p LVM HKEVP GEVP EGP BRP TEP 0.980 0.982 0.984 0.986 0.988 0.990

Simulations from EGP q q q q q q q q q q q q q q q q p LVM HKEVP criterion between these two models may then be the HKEVP over the BRP. As for the estimation of the joint probability, the performances of these two models are quite comparable and depend strongly on the model M S used to produce the simulations. One drawback of the HKEVP comes from its Bayesian nature: the independence may never be totally reached numerically (see Section 3.3), as it can be seen in upper-left plot of Figure 5, the best choice to describe a bivariate structure between these two models is then the BRP.

The EGP of [START_REF] Schlather | Models for stationary max-stable random fields[END_REF] and the GEVP of [START_REF] Smith | Max-stable processes and spatial extremes[END_REF] are the most unsatisfactory models. Indeed, in Figure 4 the GEVP has always a strong RMSE and his bias in Figure 3 is reasonably low only when M S is the LVM. The EGP tends to overestimate the extrapolation of the 100-years return level, except when M S is the EGP itself, showing a lack of robustness toward the simulated data set. Concerning the estimation of the second criterion of comparison which represents the spatial dependence structure, it can be seen from Figure 5 that the GEVP always find spatial independence of the data, with a joint probability p estimated at 0.9801. This model seems to have difficulties to detect any trace of spatial dependence, even when M S is the TEP where the exact value of p is relatively high. The estimations of this joint probability with the EGP are nearly always degenerated at p ≈ 0.98299. This enlightens the main drawback of this model discussed in Section 2.2.2, that is, the fact that θ(h) 1.838 and therefore that the EGP cannot account for spatial independence even if the distance h tends to infinity. The only exception where estimations are not degenerated at the lower bound case p is when M S is the EGP itself, which is of limited use.

3.3.

Beyond the scope of this study: some additional comparisons In the previous section, a comparative study has been driven between six spatial models based on two risk measures. The first one assesses the marginal behavior through the extrapolation of a 100-years return level, while the second evaluates the dependence structure with a joint probability. While these two criteria summarize the efficiency of the models for spatial extreme value analyses, it may seem insufficient to just look at these values for a valid competition. In this section, we discuss other properties of the six spatial models that have not been fully taken into consideration yet. The aim is to guide the practitioner for a better choice of model, depending on the pursued goal.

Two Bayesian models have been used: the LVM and the HKEVP. The results obtained when fitting these two models are Markov chains that, if convergence is assumed, represent a posterior distribution for each parameter or for a functional of parameters like the 100-years return level. To be feasible, the competition of the previous section needed to focus on point estimations of the two criteria and it has been chosen to take the median of the posterior distributions. This choice may therefore seem rather restrictive for such models that can provide much more information on a given parameter. An example is the case of the HKEVP, when the model M S used to produce the simulation is the LVM and therefore the spatial independence is assumed. To accept independence, the parameter α in the HKEVP must be equal to 1, its upper bound. Numerically, the Markov chain cannot afford for α = 1 exactly since all values of the random effect, which is updated at each iteration, should then be degenerated to 1. Taking the median of the chain as point estimation for α implies in this case a slight underestimation of its true value, which can explain the behaviour of this model on upper-left plot of Figure 5.

However, the Bayesian models present two main drawbacks at the inference step. First, they have generally more entry parameters than for the non-Bayesian ones. For instance, the user of the function latent must provide values that control initial steps, prior distributions, random walks of proposal distributions and, above all, the number of iterations to use in the algorithm for assessing convergence (burn-in period) and to provide a satisfactory sample of the posterior distribution. In the function hkevp.fit, default values are available for these arguments, though they should be considered with care. Second, these models are more time consuming, especially when the length desired of the posterior samples is high. The HKEVP is the heaviest in this case and parallel computing is advised when using this method. The choice of the number L of knots used in the HKEVP has to be seen as a trade-off between efficient estimation and computational burden: all the values of the random effect are updated at each MCMC step, which represents L × T parameters. This has been studied in [START_REF] Reich | A hierarchical max-stable spatial model for extreme precipitation[END_REF] and the conclusions were roughly that too few knots may lead to a larger bias in the estimation of the GEV parameters, while too many knots than necessary does not improve significantly these estimations.

The HKEVP has another drawback: when the exact value of α is near 0, that is, the case of very strong spatial dependence. In this case, convergence of the Markov chains are very slow due to the fact that the values of the random effect A are nearly uniformly distributed over R + \{0}. However, it has to be noted that this feature is purely theoretical, since annual maxima of a natural phenomenon such as precipitation has apparently no reason to show such strong dependence.

Finally, it can be remarked that the joint probability, i.e. the second criterion of comparison, has been evaluated on a set of only two sites. The reason for this is that closed-form for the exponent measure are not always available. For instance, [START_REF] Genton | On the likelihood function of gaussian max-stable processes[END_REF] shows that the joint cdf of the GEVP cannot be explicitly given for n > d + 1 sites, where d is the dimension of the space. In our case we have d = 2, which allows then to compute a trivariate joint probability at the most for this model. For the EGP, the BRP and the TEP, the joint cdf of an arbitrary number of sites n involves a multivariate Student or normal cdf that has to be numerically computed. Conversely, the joint cdf of the HKEVP has an explicit and simple form in any dimension.

As a summary, Table 2

Discussion and conclusion

In this article, six models for spatial extreme values have been put in competition over two risk measures that represent usual interest in application: the extrapolation of a 100-years return level at an ungauged site and the estimation of an extreme joint probability.

Results from Section 3.2.4 show dissimilarities between models and tend to discard some of them, depending on what is the main objective. On the one hand, if the interest lies in the estimation of the marginal effect, one should prefer the LVM of Davison et al. (2012) and may also consider the TEP of [START_REF] Opitz | Extremal t processes: Elliptical domain of attraction and a spectral representation[END_REF]. On the other hand, if the goal involves the modelling of the joint dependence structure, the TEP is the best choice but the BRP and the HKEVP may also provide reliable estimations. In any case, the GEVP of [START_REF] Smith | Max-stable processes and spatial extremes[END_REF] and the EGP of [START_REF] Schlather | Models for stationary max-stable random fields[END_REF] show results that are generally unsatisfactory on the two criteria defined for the competition.

However, this comparison of models has been made under given circumstances that may influence the general conclusions. Namely, the set of precipitation data we used is rather sparse over a large domain, which results in weak spatial dependence as already pointed out in Section 3.1. But we took care of avoiding subjective choices for the simulation plan, relying only on the real precipitation data set by using the parametric bootstrap procedure.

As discussed in Section 3.3, the comparison between models may generally be more complex because of several features that characterize each of them. For instance, the HKEVP is the only one which can give an explicit formulation of the dependence structure for an arbitrary set of sites, therefore allowing a conditional sampling of the yearly maxima process. However, his inference is less tractable: it involves a lot of arbitrary choices like the positions of the knots, and it demands more computational resources than the other five models to be properly fitted.

This paper can be regarded as a practical guide when fitting annual maxima of precipitation data. Depending on the question of interest, the user has to choose between max-stable (TEP and BRP) or hierarchical max-stable (LVM and HKEVP) models. Moreover, the simulation plan is based on real precipitation data so that the comparison made here does not suffer from subjectivity.

  mean function of the latent processes µ(•) and σ(•) is written as a linear combination of covariates c(•) with coefficients β . (•) ∈ R. Example of covariates are the longitude, latitude and elevation of position s. The random part ε • (•) is assumed to be a stationary zero-mean Gaussian process with correlation function ρ(•). As an example, Davison et al. (2012) considers the exponential form ρ

Figure 1 :

 1 Figure 1: Map of France with locations of the meteorological stations and a color code indicating the elevation. The two stations represented by triangles (see south west) are considered later on to compute a joint probability.

Figure 2 :

 2 Figure 2: Map of France with positions of the meteorological stations and a color code indicating the values of the GEV location (left) and scale (right) parameters, estimated with the BRP.

Figure 3 :

 3 Figure 3: Boxplots per site of the biases {B(s i , M S |M F )} i=1,...,n for the estimation of the 100-years return level.

Figure 4 :

 4 Figure 4: Boxplots per site of the RMSEs {R(s i , M S |M F )} i=1,...,n for the estimation of the 100-years return level.

Figure 5 :

 5 Figure 5: Boxplots of the estimated joint probabilities { p (k) (M S |M F )} k=1,...,K .

Table 1 :

 1 Table1summarizes the estimated values of the unknown dependence parameters for the different models described in Sections 2.2 to 2.3.2. Recall that LVM does not appear in the list since it assumes independence. We observe that, under each model, the estimated spatial dependence structure is close Estimated dependence parameters (and standard deviation when available) for the HKEVP, GEVP, EGP, BRP and TEP when fitted on the ECA&D precipitation data set described in Section 3.1.1.

	Model Parameter Estimation (sd)
	HKEVP	α	0.846 (0.021)
		τ	1.21 (0.096)
	GEVP	Σ 11	0.064
		Σ 12	0.035
		Σ 22	0.02
	EGP	η EGP	0.999
		λ EGP	2.965
		ν EGP	0.839
	BRP	λ BRP	0.103 (0.029)
		ν BRP	0.723 (0.092)
	TEP	δ TEP	15.899 (17.09)
		η TEP	0.047 (0.085)
		λ TEP	3.819 (5.699)
		ν TEP	0.941 (0.35)

  provides a visual assessment of the six models over the two criteria and the points discussed above.

		LVM HKEVP GEVP EGP BRP TEP
	Marginal extrapolation	≈	≈
	Joint probability	≈
	Bayesian approach
	Fast program	≈	≈
	Explicit multivariate cdf

Table 2 :

 2 Sketch of the characteristics of the six spatial max-stable models.

A preprocessing of the data set aimed at suppressing doubloons, merging data from the same locations, and keeping stations with elevation less than 800 meters to avoid very specific extreme behaviour.

Software R http://www.r-project.org/

Homepage http://www4.stat.ncsu.edu/∼reich/

Details for HKEVP and LVM fit: 30.000 iterations run, the first 10.000 burned to assess convergence. Median of the last 20.000 is taken as point estimate.

Acknowledgements

We would like to thank many contributors without whom this paper would never exist. First of all, we thank the authors of the SpatialExtremes R package, in particular Mathieu Ribatet for his very explicit functions around the GEVP, EGP, TEP and BRP models and the great help he provided us.

16

We also thank a lot Brian Reich and Benjamin Shaby for their model (HKEVP) and the indications they have given, allowing moreover the implementation of the R package hkevp.

This paper has been written during the PhD thesis of the first author. His thesis has been financed by EDF. We would like to thank several researchers from EDF R&D, namely Anne Dutfoy, Marie Gallois, Thi Thu Huong Hoang and Sylvie Parey, for numerous fruitful discussions that improved substantially this work.

Finally, we are grateful to the ECA&D website and to R project which both provide free material.

Concerning the models GEVP, EGP, BRP and TEP, the linear expression (7) allows to extrapolate directly the GEV parameters at s * . The associated 100-years return level can then be computed thanks to (8). For the two hierarchical models, namely the LVM and the HKEVP, there are several ways to estimate this criterion. We choose to extrapolate the value of the GEV parameters at s * for each MCMC step r, using a simple kriging estimator. Formula (8) leads to the associated 100-years return level for each MCMC step r. The point estimate is then taken as a functional (the median in our case) of the sample {y (r) 100 (s * )} r=1,...,R , where R = 20.000 is the number of iterations kept in the fit.

Estimating the comparative criteria under the six models

Let M F be a fitting model chosen among LVM, HKEVP, GEVP, EGP, BRP, TEP. The two criteria defined in Section 3.2.2 are estimated via M F over each set of simulated data Y (k) M S . We denote by y 100 (s i , M S ) the exact value of the 100-years return level evaluated at s i for i ∈ {1, ..., n} with the parameters according to M S and by p(M S ) the joint probability corresponding to the second criterion evaluated with the model M S over real data. Following this notation, we respectively denote by y 100 (k) M S . Note that the first criterion y 100 (s i , M S ) is estimated at one of the meteorological stations s i of the ECA data set, though it has been decided in Section 3.2.2 to study the spatial extrapolation of the GEV margins. Therefore, to obtain the k-th estimation of this criterion at s i , the six spatial models are fitted over the simulated set

M S minus the observations at s i , which is then treated as a target site. This procedure is repeated for each i ∈ {1, ..., n}. In constrast, the estimation of the joint probability p is obtained by fitting M F over the complete set of 39 stations. This choice has been motivated by the fact that different conclusions can be drawn depending on target position s * where the first criterion is defined, while it is not the case for the joint probability. Indeed, this second criterion only depends on the distance between the two sites (s 1 , s2 ) for four of the six models, that are said isotropic. The two exceptions are the GEVP and the HKEVP, but no genuine differences have been observed by choosing a different pair (s 1 , s2 ) to compute the joint probability for these two models.

Results

After a quite intensive simulation from M S ∈ {LVM, HKEVP, EGP, TEP} and fitting procedure via M F ∈ {LVM, HKEVP, GEVP, EGP, BRP, TEP}, a pseudo-sample of estimated criteria is produced. To display the results of the estimated 100-years return levels, we focus on biases {B(s i , M S |M F )} i=1,...,n and root-mean-square errors (RMSE) {R(s i , M S |M F )} i=1,...,n per site, respectively defined by:

for i ∈ {1, ..., n}. These are respectively given in Figure 3 and