Clément Dombry

Mathieu Ribatet

Stilian Stoev

♦ March

Probabilities of concurrent extremes

Keywords: Max-stable process, Poisson point process, Slyvniak formula, Concurrence, Kendall's τ, Temperature

The statistical modelling of spatial extremes has recently made major advances. Much of its focus so far has been on the modelling of the magnitudes of extreme events but little attention has been paid on the timing of extremes. To address this gap, this paper introduces the notion of extremal concurrence. Suppose that one measures precipitation at several synoptic stations over multiple days. We say that extremes are concurrent if the maximum precipitation over time at each station is achieved simultaneously, e.g., on a single day. Under general conditions, we show that the finite sample concurrence probability converges to an asymptotic quantity, deemed extremal concurrence probability. Using Palm calculus, we establish general expressions for the extremal concurrence probability through the max-stable process emerging in the limit of the componentwise maxima of the sample. Explicit forms of the extremal concurrence probabilities are obtained for various max-stable models and several estimators are introduced. In particular, we prove that the pairwise extremal concurrence probability for max-stable vectors is precisely equal to the Kendall's τ. The estimators are evaluated by using simulations and applied to study the concurrence patterns of temperature extremes in the United States. The results demonstrate that concurrence probability can provide a powerful new perspective and tools for the analysis of the spatial structure and impact of extremes.

Introduction

While most of the time extreme value analysis focuses on the magnitude of extreme events, i.e., how large extremes events are, little interest has been paid to their genesis. This paper tries to fill in this gap by looking at what we shall call concurrency of extremes, e.g., have two locations been impacted by the same extreme event or was it a consequence of two different ones? For example, one could observe daily rainfall at various weather stations and would like to quantify the risk that the rainfall extremes over a spatial domain are due to a single extreme event, i.e., a large storm, affecting the entire area. Although potentially rare, such events have great socio-economic consequences and their probabilities should be assessed precisely.

More formally, given a sequence X 1 , . . . , X n of independent copies of a stochastic process X defined on a compact set X ⊂ R d , d ≥ 1, we say that extremes are sample concurrent at locations s 1 , . . . , s k ∈ X, k ≥ 2, if max i=1,...,n X i (s j ) = X (s j ), j = 1, . . . , k,

for some ∈ {1, . . . , n}. Clearly this means that only the observation X contributes to the pointwise maxima at locations s 1 , . . . , s k . It occurs with probability p n (s 1 , . . . , s k ) = P for some ∈ {1, . . . , n} : max i=1,...,n X i (s j ) = X (s j ), j = 1, . . . , k ,

1 henceforth referred to as sample concurrence probability.

Provided that X has continuous margins it is not difficult to see that p n (s 1 , . . . , s k ) = nE F {X(s 1 ), . . . , X(s k )} n-1 , where F is the multivariate cumulative distribution of {X(s 1 ), . . . , X(s k )}. Interestingly, concurrence of extremes event is invariant under increasing transformations of the marginals so that the concurrence probability does not depend on the marginal distributions of X but only on its dependence structure, i.e., the copula C associated to F.

One drawback of the sample concurrence probability p n (s 1 , . . . , s k ) is that it varies with the number of observations n. Surprisingly, however, we prove in Theorem 1 below that, under mild regularity conditions, this quantity stabilizes to a universal large sample limit p n (s 1 , . . . , s k ) -→ p(s 1 , . . . , s k ), n → ∞.

(3)

Throughout this paper, we will call the limiting probability p(s 1 , . . . , s k ) the extremal concurrence probability. This asymptotic quantity is naturally expressed in terms of a max-stable process η emerging in the limit of the normalized maxima in (1), as n → ∞. A direct, intuitive, and equivalent definition of the extremal concurrence probability can be given in terms of the spectral representation of this max-stable process η. Following de [START_REF] De Haan | A spectral representation for max-stable processes[END_REF], [START_REF] Penrose | Semi-min-stable processes[END_REF] and [START_REF] Schlather | Models for stationary max-stable random fields[END_REF], let

η(s) = max i≥1 ζ i Y i (s) s ∈ X, (4)
where {ζ i : i ≥ 1} are the points of a Poisson process on (0, ∞) with intensity measure ζ -2 dζ, Y i are independent copies of a non negative stochastic process with continuous sample paths such that E{Y(s)} = 1 for all s ∈ X and E{sup s∈X Y(s)} < ∞. It is often more convenient to rewrite (4) into

η(s) = max ϕ∈Φ ϕ(s), s ∈ X, (5)
where Φ = {ϕ i : i ≥ 1} with ϕ i = ζ i Y i is a Poisson point process on C 0 , the space of non negative continuous functions on X. Within this framework, we now say that extremes are concurrent at s 1 , . . . , s k ∈ X if η(s j ) = ϕ (s j ), j = 1, . . . , k,

for some ≥ 1, and similarly to the definition of the sample concurrence probability (2), the extremal concurrence probability is defined by p(s 1 , . . . , s k ) = P for some ≥ 1 : η(s j ) = ϕ (s j ), j = 1, . . . , k .

When k = 2 the extremal concurrence probability p(s 1 , s 2 ) coincides with the dependence measure considered in [START_REF] Weintraub | Sample and ergodic properties of some min-stable processes[END_REF] to study mixing properties of max-stable processes. Another well known measure of dependence is the pairwise extremal coefficient [START_REF] Schlather | A dependence measure for multivariate and spatial extremes: Properties and inference[END_REF][START_REF] Cooley | Variograms for spatial max-stable random fields[END_REF]

θ(s 1 , s 2 ) = -log P{η(s 1 ) ≤ 1, η(s 2 ) ≤ 1}, s 1 , s 2 ∈ X. (8)
Interestingly, the extremal concurrence probability and the pairwise extremal coefficient share connections. For instance Proposition 5.1 in [START_REF] Stoev | On the ergodicity and mixing of max-stable processes[END_REF] implies

1 2 {2 -θ(s 1 , s 2 )} ≤ p(s 1 , s 2 ) ≤ 2{2 -θ(s 1 , s 2 )}, (9)
and we shall see later that the properties of the extremal concurrence probability are similar to that of the pairwise extremal coefficient. The structure of the paper is as follows. In Section 2 we make connections between sample concurrence probabilities and their extremal counterparts and derive their properties. Section 3 gives closed forms for various parametric max-stable models, and Section 4 introduces various estimators for the sample/extremal concurrence probabilities. The proposed estimators are then analyzed in a simulation study in Section 5 and applied to US continental temperature extremes in Section 6.

Concurrence of extremes

In this section we show that sample concurrence probabilities converge to extremal concurrence ones under rather mild domain of attraction conditions. We then provide formulas for the extremal concurrence probability based on the spectral representation of the associated max-stable process and establish their basic properties.

Sample and extremal concurrence

Concurrence of extremes can be defined through the more general notion of a hitting scenario, which reflects precisely how many different events contribute to the componentwise maximum. Let X 1 , . . . , X n be a sequence of independent copies of a stochastic process X defined on X and s 1 , . . . , s k ∈ X be different locations. We suppose that X has continuous marginals to ensure that {X 1 (s j ), . . . , X n (s j )} has no ties almost surely and that the maximum is uniquely reached. Let M n (s) = max i=1,...,n X i (s) be the componentwise maximum and consider the sets C i = { j : M n (s j ) = X i (s j )}, i = 1, . . . , n, that account for the location where the i-th component X i dominates the rest. Some of these sets may be empty, but from the above discussion, with probability one the non-empty ones are disjoint and form a random partition of {1, . . . , k}. This partition π n = {C i : C i ∅} will be referred to as the sample hitting scenario.

By analogy with extremal concurrence, one can define an extremal hitting scenario associated to a max-stable process by using the underlying Poisson point process [START_REF] Wang | Conditional sampling for spectrally discrete max-stable random fields[END_REF]Dombry et al., 2013;Dombry and Éyi-Minko, 2013]. More precisely, for η as in ( 5), the extremal hitting scenario π is defined as the random partition of {1, . . . , k} such that two indices j 1 , j 2 ∈ {1, . . . , k} are in the same component of π if and only if arg max

i≥1 ϕ i (s j 1 ) = arg max i≥1 ϕ i (s j 2 ).
Whatever type of concurrence is considered, i.e., sample concurrence (1) or extremal concurrence (6), extremes are said concurrent if and only if π n = {1, . . . , k} or π = {1, . . . , k}. The next theorem shows the convergence of the sample hitting scenario to the extremal one.

Theorem 1. Assume that [h 1 {X(s 1 )}, . . . , h k {X(s k )}] belongs to the maximum domain of attraction of {η(s 1 ), . . . , η(s k )}, for some strictly increasing deterministic functions h i , i = 1, . . . , k. Then, the sample hitting scenario π n converges weakly as n → ∞ to the extremal hitting scenario π associated to the max-stable process η.

Proof. Let X 1 , X 2 , . . . be a sequence of independent copies of X observed at some locations s = (s 1 , . . . , s k ) ∈ X k , k ≥ 2. For brevity we will use componentwise algebra (sum, maximum, etc.) and write X(s) = {X(s 1 ), . . . , X(s k )}. Since the hitting scenario is invariant to strictly increasing deterministic transformations of the marginals, we can assume without loss of generality that h i (x) = x, i = 1, . . . , k. By assumptions the componentwise maxima M n (s) = {max i=1,...,n X i (s 1 ), . . . , max i=1,...,n X i (s k )} converge in distribution

M n (s) -b n (s) a n (s) -→ η(s), n → ∞,
where a n (s) > 0 and b n (s) ∈ R. It is well known that the above convergence is equivalent to the convergence of the point process

Φ n = X i (s) -b n (s) a n (s) : i = 1, . . . , n ⊂ R k . to a Poisson point process Φ on [0, ∞) k \{0}
where the convergence is meant in the space of point measures M p ([0, ∞] k \ {0}) equipped with the metric of vague convergence [START_REF] Resnick | Extreme Values, Regular Variation and Point Processes[END_REF]. Consider the mapping

Π : M p ([0, ∞] k \ {0}) -→ P k Ψ -→ Π(Ψ),
where P k is the set of all possible partitions of {1, . . . , k} and Π(Ψ) is the hitting scenario associated to the collection of functions Ψ = {ψ i : i ∈ I}, I ⊆ N. Clearly, π n = Π(Φ n ) and π = Π(Φ). Since the map Π is well defined and continuous at each point Ψ ∈ M p ([0, ∞] k \ {0}) for which the maxima are uniquely defined, see Appendix A.1, we can apply the continuous mapping theorem to show that the weak convergence Φ n → Φ entails the weak convergence π n → π.

Remark. A direct consequence of Theorem 1 is that, provided the stochastic process X is in the domain of attraction of some max-stable process η, the sample concurrence probability converges to its extremal counterpart, i.e.,

p n (s 1 , . . . , s k ) = P [π n = {1, . . . , k}] -→ P [π = {1, . . . , k}] = p(s 1 , . . . , s k ), n → ∞,
and proves (3).

General formulas and properties of extremal concurrence probabilities

The following theorem gives an expression for the extremal concurrence probability p(s 1 , . . . , s k ).

Theorem 2. We have

p(s 1 , . . . , s k ) = E Y        E Ỹ max j=1,...,k Ỹ(s j ) Y(s j ) -1        , (10)
where Y and Ỹ are independent copies of the stochastic process appearing in (4).

Proof. Let Λ denotes the intensity mesaure of the C 0 -valued Poisson point process Φ in (5) given by

Λ(A) = ∞ 0 P(ζY ∈ A)ζ -2 dζ,
for all Borel set A ⊂ C 0 . We have

p(s 1 , . . . , s k ) = P {∃ϕ ∈ Φ : ϕ(s 1 ) = η(s 1 ), . . . , ϕ(s k ) = η(s k )} = C 0 P f (s j ) > η(s j ), j = 1, . . . , k Λ(d f ) (11) = E Y ∞ 0 P η(s j ) < ζY(s j ), j = 1, . . . , k ζ -2 dζ = E Y ∞ 0 exp -E Ỹ max j=1,...,k Ỹ(s j ) ζY(s j ) ζ -2 dζ = E Y        E Ỹ max j=1,...,k Ỹ(s j ) Y(s j ) -1       
where η and Ỹ are independent copies of η and Y respectively. Note that the second equality uses Slyvniak's formula, the fourth one the cumulative distribution of the max-stable process η and the last one the expectation of an inverse exponential random variable.

Remark. By the seminal paper of de Haan [1984] (see also [START_REF] Stoev | Extremal stochastic integrals: a parallel between max-stable processes and α-stable processes[END_REF]; [START_REF] Kabluchko | Spectral representations of sum-and max-stable processes[END_REF]), any continuous in probability max-stable process can be represented as

{η(s) : s ∈ X} d = max i≥1 ζ i f s (u i ) : s ∈ X , (12)
where { f s : s ∈ X} is a collection of non-negative integrable functions on the space (U, U, ν).

Here {(ζ i , u i ) : i ≥ 1} is a Poisson point process on (0, ∞) × U with intensity ζ -2 dζν(du).
The functions { f s : s ∈ X} are known as spectral functions of η and (12) as de Haan's spectral representation. When ν is a probability measure, one can view Y(s) = f s as random variables on the probability space (U, U, ν) and then (12) becomes (4). Conversely, any representation (12) can be cast in the form (4) with a change of variables. Depending on the context one representation may be more convenient than the other. In terms of (12), the concurrence probability formula in (10) becomes

p(s 1 , . . . , s k ) = U U max j=1,...,k f s j (ũ) f s j (u) ν(dũ) -1 ν(du), (13)
and the proof is essentially the same.

One could expect from Definition 6 and Theorem 2 that the extremal concurrence probability depends on the distribution of the spectral process Y in (4) or the choice of spectral functions in (12). These representations are not unique, but we will see in the theorem below that the extremal concurrence probability p(s 1 , . . . , s k ) depends only on the distribution of the max-stable process η and not on the choice of the specific spectral representation.

Theorem 3. For s 1 , . . . , s k ∈ X, k ≥ 2, we have

p(s 1 , . . . , s k ) = k r=1 (-1) r J⊆{1,...,k} |J|=r E η log P η η(s j ) ≤ η(s j ), j ∈ J , ( 14
)
where η is an independent copy of η. In particular when k = 2,

p(s 1 , s 2 ) = 2 + E η log P η η(s j ) ≤ η(s j ), j = 1, 2 . ( 15
)
Proof. Starting from (11) and applying the inclusion-exclusion formula, we have

p(s 1 , . . . , s k ) = E η Λ f (s j ) > η(s j ) for some j ∈ {1, . . . , k} = E η               k r=1 J⊆{1,...,k} |J|=r (-1) r+1 Λ f (s j ) > η(s j ), j ∈ J               . Since the cumulative distribution function of η is P η η(s j ) ≤ η(s j ), j ∈ J = exp -Λ f (s j ) > η(s j ), j ∈ J , we get p(s 1 , . . . , s k ) = k r=1 (-1) r J⊆{1,...,k} |J|=r E η log P η η(s j ) ≤ η(s j ), j ∈ J .
The simplification when k = 2 is straightforward because when |J| = 1 we have

E η log P η η(s j ) ≤ η(s j ), j ∈ J = E η -η(s j ) -1 = -1,
as η(s j ) -1 is a standard exponential random variable.

In the remaining part of this section, we investigate some properties of the extremal concurrence probabilities. Surprisingly, although the two notions are different, we encounter strong similarities with the extremal coefficient (8). We recall that the extremal coefficient θ(s 1 , s 2 ) takes values in [1,2], the lower and upper bounds correspond to perfect dependence and independence respectively. The next proposition states a similar result for the extremal concurrence probability.

Proposition 1. For all s 1 , s 2 ∈ X, we have i) p(s 1 , s 2 ) = 0 if and only if η(s 1 ) and η(s 2 ) are independent;

ii) p(s 1 , s 2 ) = 1 if and only if η(s 1 ) and η(s 2 ) are almost surely equal.

The proof uses the following generalization and improvement of the upper bound in (9).

Lemma 1. For all s 1 , . . . , s k ∈ X, k ≥ 2, we have p(s 1 , . . . , s k ) ≤ E min j=1,...,k Y(s j ) .

Proof. In the context of Theorem 2, we have (by conditioning on Y)

E Ỹ max j=1,...,k Ỹ(s j ) Y(s j ) ≥ max j=1,...,k Y(s j ) -1 E Ỹ Ỹ(s j ) = min j=1,...,k Y(s i ) -1
, since E Ỹ { Ỹ(s j )} = 1. This, in view of (10) implies the desired result.

Proof of Proposition 1. Equation (9) implies that p(s 1 , s 2 ) = 0 if and only if θ(s 1 , s 2 ) = 2 which is equivalent to the independence of η(s 1 ) and η(s 2 ). When p(s 1 , s 2 ) = 1, Lemma 1 entails Y(s 1 ) = Y(s 2 ) almost surely so that η(s 1 ) = η(s 2 ) almost surely. It is easy to prove the converse implication: if η(s 1 ) and η(s 2 ) are almost surely equal, the same holds for Y(s 1 ) and Y(s 2 ) so that p(s 1 , s 2 ) = 1.

Interestingly p(s 1 , . . . , s k ) can be expressed via the extremal coefficients of another max-stable process.

Proposition 2. Let η be a simple max-stable process as defined in (4) and η, η1 , η2 , . . . independent copies of it.

Consider the simple max-stable process

ξ(s) = max i≥1 ζ i Y i (s) ηi (s) , s ∈ X, then p(s 1 , . . . , s k ) = k r=1 (-1) r+1 J⊆{1,...,k} |J|=r θ ξ (s j , j ∈ J),
where

θ ξ (s j , j ∈ J) = -log P{ξ(s j ) ≤ 1, j ∈ J}. In particular p(s 1 , s 2 ) = 2 -θ ξ (s 1 , s 2 ).
Proof. Clearly ξ is a simple max-stable process since both Y and η are non negative and E{Y(s)/ η(s)} = 1 for all s ∈ X.

We have

E η log P η η(s j ) ≤ η(s j ), j ∈ J = -E η E Y max j∈J Y(s j ) η(s j ) = log P ξ ξ(s j ) ≤ 1, j ∈ J = -θ ξ (s j , j ∈ J).
The next corollary lists some properties of the extremal concurrence probability function that closely parallel those of the extremal coefficient function. In view of Proposition 2, the proof follows as in [START_REF] Schlather | A dependence measure for multivariate and spatial extremes: Properties and inference[END_REF] or [START_REF] Cooley | Variograms for spatial max-stable random fields[END_REF].

Corollary. Let p : h → p(o, h) be an extremal concurrence probability function associated to a stationary max-stable process in X for some arbitrary origin o ∈ X and h ∈ X. Then the following assertions hold.

i) The function h → p(h) is positive semidefinite;

ii) The function h → p(h) is not differentiable at the origin unless p(h) = 1 for all h ∈ X;

iii) If d ≥ 1 and if η is isotropic, then h → p(h) has at most a jump at the origin and is continuous elsewhere;

iv) {2 -p(h 1 + h 2 )} ≤ {2 -p(h 1 )}{2 -p(h 2 )} for all h 1 , h 2 ∈ X; v) {2 -p(h 1 + h 2 )} α ≤ {2 -p(h 1 )} α + {2 -p(h 2 )} α -1 for all h 1 , h 2 ∈ X and 0 ≤ α ≤ 1; vi) {2 -p(h 1 + h 2 )} α ≥ {2 -p(h 1 )} α + {2 -p(h 2 )} α -1 for all h 1 , h 2 ∈ X and α < 0.
We conclude this section with an unexpected result that relates the bivariate extremal concurrence probability with the well known Kendall's τ.

Theorem 4. For any max-stable process η, we have p(s 1 , s 2 ) = τ where τ = E sign{η(s 1 ) -η * (s 1 )}sign{η(s 2 ) -η * (s 2 )} is the Kendall's τ of {η(s 1 ), η(s 2 )} and η * is an independent copy of η.

Proof. Let W = F{η(s 1 ), η(s 2 )} where F is the bivariate cumulative distribution function of {η(s 1 ), η(s 2 )}. From (15) we have p(s 1 , s 2 ) = 2 + E log W . But since {η(s 1 ), η(s 2 )} is a bivariate max-stable random vector, we know that [START_REF] Ghoudi | Propriétés statistiques des copules de valeurs extrêmes bidimensionnelles[END_REF]] and hence, after some simple calculations, p(s 1 , s 2 ) = τ.

P(W ≤ w) = w -(1 -τ)w log w, 0 ≤ w ≤ 1 [

Formulas for extremal concurrence probabilities

In this section we gather formulas for the extremal concurrence probabilities for some popular models of max-stable random vectors and processes. As we will see, it is not always possible to get explicit formulas, and in such situations, we propose to use Monte-Carlo methods.

Closed forms

Example 1 (Logistic model). The concurrence probability for the k-variate logistic model, i.e., with cumulative distribution

F(z 1 , . . . , z k ) = exp          -         k j=1 z -1/α j         α          , 0 < α ≤ 1, z 1 , . . . , z k > 0, is p(s 1 , . . . , s k ) = k-1 j=1 (1 -α/ j).
Recall that for this model independence is reached when α = 1 while perfect dependence occurs as α ↓ 0 and, as expected, for such situations we have p(s 1 , . . . , s k ) = 0 and p(s 1 , . . . , s k ) = 1 respectively.

Proof. It is not difficult to see that the multivariate logistic model corresponds to the case where Y in ( 4) is a pure noise process with margins such that P{Y(s)

< y} = exp[-{Γ(1 -α)y} -1/α ]
where Γ is the Gamma function. Using Theorem 2, we have

p(s 1 , . . . , s k ) = E Y        E Ỹ max j=1,...,k Ỹ(s j ) Y(s j ) -1        = E Y -log F{ Ỹ(s 1 ), . . . , Ỹ(s k )} -1 = E Y                   k j=1 Y(s j ) -1/α          -α          = Γ(k -α) Γ(k)Γ(1 -α) = k-1 j=1 (1 -α/ j),
where the fourth equality used the fact that k j=1 Y(s j ) -1/α is a Gamma random variable with scale Γ(1 -α) 1/α and shape k for which negative moments are known.

Example 2 (Max-linear model). Consider the max-linear model η(s) = max m=1,...,n ϕ m (s)Z m , where Z 1 , . . . , Z n are independent unit Fréchet random variables and some functions {s → ϕ m (s) ∈ C 0 , m = 1, . . . , n} such that n m=1 ϕ m (s) = 1 for all s ∈ X. We have i) The concurrence probability equals

p(s 1 , . . . , s k ) = n =1 p (s 1 , . . . , s k ), p (s 1 , . . . , s k ) =        n m=1 max j=1,...,k ϕ m (s j ) ϕ (s j )        -1 , (16)
with the convention that 0/0 = 0, a/0 = ∞ if a > 0, and 1/∞ = 0.

ii) The probability that component dominates at sites s 1 , . . . , s k is given by the term p in ( 16), i.e.,

p (s 1 , . . . , s k ) = P η(s j ) = ϕ (s j )Z , j = 1, . . . , k (17)
Proof. Part i) is an immediate consequence of (13). Indeed, let U = {1, . . . , n} be equipped with the counting measure

ν{1} = • • • = ν{n} = 1. By taking f s (m) = ϕ m (s), m ∈ U,
we obtain that the max-linear model has the representation (12). The integral expression of the concurrence probability in (13) then becomes a sum of the terms p in (16).

Part ii) shows an intriguing fact that the concurrence probability p(s 1 , . . . , s k ) for the max-linear model is the sum of the probabilities that one of the n components dominates the rest. Indeed, by the max-stability property and the independence of the unit Fréchet random variables Z m 's, we have that the right-hand side of (17) equals

P max m max j=1,...,k ϕ m (s j ) ϕ (s j ) Z m ≤ Z = P(aZ 1 ≤ Z 2 ),
where a = m max j=1,...,k ϕ m (s j )/ϕ (s j ). Equation ( 17) follows from the fact that P(aZ

1 ≤ Z 2 ) = (1 + a) -1 , a ≥ 0.
Example 3 (Chentsov random fields). Suppose that the process Y(s) = 1 A (s), where A ⊂ R d is a random set. Then, by analogy with the theory of symmetric α-stable process [Samorodnitsky and Taqqu, 1994, Chap. 8], the max-stable process η(s) = max

i≥1 ζ i 1 A i (s), s ∈ X,
where Y i ≡ 1 A i are independent copies of Y ≡ 1 A , will be referred to as a Chentsov max-stable random field on X.

For a Chentsov-type max-stable process we have

p(s 1 , . . . , s k ) = P ({s 1 , . . . , s k } ⊂ A | {s 1 , . . . , s k } ∩ A ∅) , (18)
or less formally that the extremal concurrence probability is the conditional probability that the sites s 1 , . . . , s k are covered by the random set A given that at least one of the sites is covered.

Proof. The result is an immediate consequence of Theorem 2. Since the outer expectation in (10) can be restricted to the event {min j=1,...,k Y(s j ) > 0}, we have

p(s 1 , . . . , s k ) = E Y        E Ỹ max j=1,...,k 1 Ã(s j ) 1 A (s j ) -1 1 {min j=1,...,k 1 A (s j )>0}        = E Y        E Ỹ max j=1,...,k 1 Ã(s j ) -1 min j=1,...,k 1 A (s j )        (19) = E min j=1,...,k 1 A (s j ) E max j=1,...,k 1 A (s j )
,

where in the second relation we used the fact that 1 {min j=1,...,k 1 A (s j )} = min j=1,...,k 1 A (s j ) and that the event min i=1,..

.,k 1 A (s i ) > 0 is equivalent to {1 A (s 1 ) = 1, . . . , 1 A (s k ) = 1}. This proves (18).
Example 4 (Extremal processes). The max-stable process {η(s) : s ∈ [0, 1]} is an extremal process if it has stationary and independent max-increments, i.e., {η(s 1 ), . . . , η(s

k )} d = [s 1 Z 1 , max{s 1 Z 1 , (s 2 -s 1 )Z 2 }, . . . , max{s 1 Z 1 , . . . , (s k -s k-1 )Z k }] ,
where 0 < s 1 < • • • < s k and Z 1 , . . . , Z k are independent unit Fréchet random variables. It can be shown that

η(s) d = max i≥1 ζ i 1 [U i ,1] (s), s ∈ [0, 1],
where U i 's are independent U(0, 1) random variables. Using our previous result on Chentsov random fields, we have for all 0

< s 1 < • • • < s k ≤ 1 p(s 1 , . . . , s k ) = P ({s 1 , . . . , s k } ⊂ [U, 1]) P ({s 1 , . . . , s k } ∩ [U, 1] ∅) = P(U ≤ s 1 ) P(U ≤ s k ) = s 1 s k ,
where U ∼ U(0, 1). This result is not surprising since for this simple case, extremes are concurrent at locations 0

< s 1 < • • • < s k < 1 if η(s) has no jumps in the interval [s 1 , s k ].
Hence using the independence and stationarity of the max-increments, the probability of the latter event is

P{s 1 Z 1 ≥ (s k -s 1 )Z 2 } = s 1 /s k ,
where Z 1 and Z 2 are two independent unit Fréchet variables.

Example 5 (Indicator moving maxima). In the context of ( 12), if f s (u) = 1 A s (u), for some sequence of measurable deterministic sets A s , by using (13), we obtain as in (19) that

p(s 1 , . . . , s k ) = ν(∩ j=1,...,k A s j ) ν(∪ j=1,...,k A s j ) . ( 20
)
In the simple case f s (u) = 1 A (us), i.e. A s = s + A with some deterministic set A, where ν is the Lebesgue measure on R d , (20) implies

p(s, s + h) = p(h) = |A ∩ (h + A)| |A ∪ (h + A)| = c A (h) 2|A| -c A (h) ,
where |A| denotes the d-dimensional volume of A and c A (h) = |A ∩ (h + A)|. The latter function and hence the extremal concurrence probability function p(h) can then be obtained in closed form for many different sets. For example, in the case η is isotropic, i.e., A = {s ∈ R d : x ≤ r} is the centered ball of radius r > 0 in Euclidean space, using the formula for the volume of the cap, we obtain

c A ( h ) = C d r d B (d+1)/2,1/2 h (2r -h ) 2r 2 , C d = π d/2 Γ(1 + d/2) , where B a,b (x) = B(a, b) -1 x 0 u a-1 (1 -u) b-1 du is the cumulative distribution function of a Beta(a, b) random variable.

Monte-Carlo methods

It may happen that for some parametric max-stable models, explicit forms for extremal concurrence probabilities are not available but hopefully it is often possible to use Monte-Carlo methods to approximate the theoretical extremal concurrence probabilities with arbitrary precision. A naive strategy would consist in using ( 10) to devise a Monte-Carlo estimator, but it is wiser to take advantage of the closed forms of max-stable processes cumulative distributions, i.e.,

P η(s j ) ≤ z j , j = 1, . . . , k = exp{-V s 1 ,...,s k (z 1 , . . . , z k )}, z 1 , . . . , z k > 0,
where V s 1 ,...,s k is an homogeneous function of order -1. Rewriting (10), we found

p(s 1 , . . . , s k ) = E Y -log P η η(s j ) ≤ Y(s j ), j = 1, . . . , k -1 = E Y V s 1 ,...,s k [Y(s 1 ), . . . , Y(s k )] -1 (21)
which can easily be estimated by sampling independent copies of Y and computing the sample mean. We can often make use of antithetic variables to get more precise estimates. Note that specific choice of the spectral process Y can lead to better strategies as we will illustrate in the following examples.

Example 6 (Brown-Resnick model). Let η be a Brown-Resnick stationary random field on X driven by a Gaussian process [START_REF] Kabluchko | Stationary max-stable fields associated to negative definite functions[END_REF]. That is, the processes Y i in (4) are equal in distribution to

Y(s) = exp{W(s) -γ(s)}, s ∈ X, ( 22
)
where W is a zero mean Gaussian random field with stationary increments and semi-variogram γ, i.e. 2γ(h

) = E{W(h) 2 } = E[{W(s + h) -W(s)} 2 ], s, h ∈ X.
For this model, the bivariate extremal concurrence probability function is given by

p(o, h) = E Φ(Z) + exp γ(h) -2γ(h)Z Φ 2γ(h) -Z -1 , (23)
where Z ∼ N(0, 1) has the standard normal distribution with cumulative distribution function Φ. As expected p(o) = 1 and p(h) → 0 as h → ∞ provided that the semi-variogram is unbounded, i.e., γ(h) → ∞ as h → ∞.

Proof. Without loss of generality, for Brown-Resnick processes we can assume that in (4) we have Y(o) = 1 almost surely. The bivariate cumulative function is given by P{η(0

) ≤ z 1 , η(h) ≤ z 2 } = exp{-V h (z 1 , z 2 )}, z 1 , z 2 > 0 with V h (z 1 , z 2 ) = 1 z 1 Φ        γ(h)/2 + 1 2γ(h) log z 2 z 1        + 1 z 2 Φ        γ(h)/2 + 1 2γ(h) log z 1 z 2        .
Equation ( 21) together with the fact that {Y(o), Y(h)} has the same distribution as {1, e

√ 2γ(h) Z-γ(h) } entails p(o, h) = E {V h [Y(o), Y(h)]} -1 = E V h 1, exp 2γ(h) Z -γ(h) -1 .
Equation ( 23) follows after straightforward simplifications.

A popular special case of the the Brown-Resnick family of models is the moving maximum storm model introduced by [START_REF] Smith | Max-stable processes and spatial extreme[END_REF] known also as the Gaussian extremal process. Consider the spectral representation (12), where U = R d and ν is the Lebesgue measure. Taking f s (u) = ϕ Σ (su), where ϕ Σ is the multivariate Normal density with zero mean and covariance matrix Σ, we obtain the max-stable process

η(s) = max i≥1 ζ i ϕ Σ (s -u i ), s ∈ X. (24)
Then, the following are true:

i) The process η belongs to the family of degenerate Brown-Resnick models in ( 22) with W(s) = s Z, s ∈ R d and Z ∼ N(0, Σ -1 ).

ii) Consequently, the concurrence probability function of η is given by ( 23) with

γ(h) = h Σ -1 h/2.
Proof. Let η be a Brown-Resnick process with variogram 2γ(h) = h Σ -1 h, that is, we can assume without loss of generality that in the spectral characterization we have

Y(s) = exp(s Σ -1 Z -s Σ -1 s) with Z ∼ N(0, Σ)
. Then for all s 1 , . . . , s k ∈ X and z 1 , . . . , z k > 0 we have

-log P{η(s j ) ≤ z j , j = 1, . . . , k} = E Z max j=1,...,k z -1 j exp s j Σ -1 Z - 1 2 s j Σ -1 s j = (2π) -k/2 |Σ| -1/2 R d max j=1,...,k z -1 j exp - 1 2 y Σ -1 y + s j Σ -1 y - 1 2 s j Σ -1 s j dy = R d max j=1,...,k z -1 j ϕ Σ (s j -y)dy.
The last relation equals the negative log cumulative distribution function of the moving maxima in (24).

Example 7 (Schlather and extremal-t processes). Let η be an extremal-t process on X, i.e., the processes Y i in (4) are equal in distribution to

Y(x) = c ν max{0, W(s)} ν , c ν = √ π2 -(ν-2)/2 Γ ν + 1 2 -1 s ∈ X,
where ν ≥ 1 and W is a stationary standard Gaussian process with correlation function ρ. The Schlather process is obtained when ν = 1. The corresponding extremal concurrence probability function is

p(o, h) = E        T ν+1 (T ) + {ρ(h) + σ(h)T } -ν T ν+1 - ρ(h) σ(h) + 1 σ(h)(ρ(h) + σ(h)T ) -1 1 {ρ(h)+σ(h)T >0}        where σ(h) = {1 -ρ(h) 2 }/(1 + ν)
and T is a Student random variable with ν + 1 degrees of freedom and cumulative distribution function T ν+1 .

Proof. For the notational convenience, we write shortly ρ = ρ(h), σ = σ(h). To obtain the desired result, we use a different spectral representation. Comparing the two cumulative distribution functions, one can show that

{η(o), η(h)} d = max i≥1 ζ i { Ỹi (o), Ỹi (h)}
with Ỹi , i ≥ 1, i.i.d. copies of the bivariate random vector

{ Ỹ(o), Ỹ(h)} =        2(1, max{0, ρ + σT } ν ), with probability 1/2 2(0, c), with probability 1/2 with c = 1 -E[max(0, ρ + σT ) ν ] such that E( Ỹ) = (1, 1). Equation (21) yields p(o, h) = E V h { Ỹ(o), Ỹ(h)} -1 = E V h {1, (ρ + σT ) ν } -1 1 {ρ+σT >0} ,
with [START_REF] Davison | Statistical modelling of spatial extremes[END_REF]]

V h (z 1 , z 2 ) = 1 z 1 T ν+1        - ρ σ + 1 σ z 2 z 1 1/ν        + 1 z 2 T ν+1        - ρ σ + 1 σ z 1 z 2 1/ν       
, and straightforward simplifications give the announced result.

4 Statistical inference and asymptotic properties

Sample concurrence probability estimators

In this section we define a sample concurrence probability estimator by blocking the data and study its basic properties as well as the optimal choice of the block-size. We conclude with a methodological improvement of the estimator based on permutation bootstrap. Let X i = {X i (s j ) : j = 1, . . . , k}, i = 1, . . . , n, be random vectors in R k , k ≥ 2. Partition the data into non-overlapping blocks of size m < n, and define the sample concurrence probability estimator

pm ≡ pm (X 1 , . . . , X n ) = 1 [n/m] [n/m] r=1 max =1,...,m 1 {maxi=1,...,m X i+(r-1)m ≤X +(r-1)m} . (25)
The max statistic above is simply an indicator of whether or not we have concurrence in the r-th block, i.e., whether one of the vectors dominates the componentwise maximum of the rest in the r-th block of size m.

Assuming that X 1 , . . . , X n are independent and identically distributed, the above estimator is the sample mean of [n/m] independent Bernoulli(p m ) random variables, where p m is as in (2) with n replaced by m. Therefore,

E( pm ) = p m , Var( pm ) = p m (1 -p m ) [n/m] ,
i.e., pm is a unbiased estimator for p m . As argued in the introduction, a major drawback of the sample concurrence probability p m is that it depends on the sample size m and it is thus more sensible to focus on the limiting extremal concurrence probability p = p(s 1 , . . . , s k ). The sample concurrence probability estimator pm is biased for p with mean squared error

MSE( pm ) = (p m -p) 2 + p m (1 -p m ) [n/m] . (26)
Although the bias term p mp is difficult to estimate in general, for the max-stable case a precise expression is available.

Proposition 3. Assume that X is a max-stable process, then p m is non-increasing in m and satisfies

0 ≤ p m -p ≤ (1 -p) m , m ≥ 1.
Furthermore, if p < 1, then there exists an integer r ∈ {1, . . . , k -1} and a positive constant c r ∈ (0, 1p], such that (p mp) ∼ c r /m r as m → ∞.

The proof expresses (p mp) via the distribution of the extremal hitting scenario and is postponed to Appendix A.2. The above result suggests that an asymptotically optimal choice of the block size can be made to minimize the rate of the mean squared error in (26). In view of Proposition 3, for the general case 0 < p < 1,

MSE( pm ) = c r m r 2 + p(1 -p)m n + o(m -2r ) + O(m/n)
Taking the derivative with respect to m we see that the optimal rate corresponds to 2rc 2 r m -2r-1 ∼ p(1p)n -1 , as n → ∞. Hence the block size that asymptotically minimizes the mean squared error is

m MSE (n) ∼ 2rc 2 r n p(1 -p) 1/(2r+1) , n → ∞. (27)
For this rate-optimal mean squared error we obtain MSE( pm ) ∝ n -2r/(2r+1) .

Remark. The constants r, c r and p in ( 27) are unknown. The precise expressions for r and c r involve multiple concurrence probabilities, i.e., when two or more events contribute to the maximum. In principle, pilot estimates of these parameters could be obtained and used as plug-ins in ( 27). Since the cases r ≥ 2 correspond to very specific dependence structures, in practice, we recommend using the conservative choice r = 1 and c r = 1.

The following result establishes the asymptotic behavior of the sample concurrence probability estimator. The proof is given in Appendix A.3.

Theorem 5. Suppose that 0 < p < 1 and let m = m(n) be such that n/m(n) → ∞ and m(n)/n 1/(2r+1) → λ ∈ (0, ∞] as n → ∞. Then n/m( pm -p) -→ N c r λ r+1/2 , p(1 -p) , n → ∞
where 1/∞ is interpreted as zero.

Remark. In Theorem 5, we encounter a typical tradeoff between rate-optimality and bias. In particular, the mean squared error optimal choice of m as in ( 27) corresponding to λ = (2rc r ) 2/(2r+1) {p(1p)} -1/(2r+1) yields the limit distribution

N        p(1 -p) 2r , p(1 -p)       
, which has non-zero mean. On the other hand, the rate sub-optimal choices where λ = ∞ yield Normal unbiased limits.

In the case p = 0, we have no optimal block size or optimal rate estimation. Still, the following asymptotic result can be useful. Theorem 6. Suppose that p = 0 and let m = m(n) be such that n/m(n) → ∞ and m(n

)/n 1/(r+1) → λ ∈ (0, ∞] as n → ∞. If λ < ∞, then (n/m) pm -→ Poisson(c r /λ r+1 ), n → ∞.
Otherwise, if λ = ∞, then P( pm = 0) → 1 as n → ∞ and hence a n pm → 0 in probability as n → ∞ for any sequence a n > 0.

Remark. When k = 2 it is possible to get an unbiased estimator for p based on a slight modification of pm . Indeed for this specific case, (36) implies that p = (mp m -1)/(m -1) and hence the estimator

pm = m pm -1 m -1 (28)
is unbiased for p as pm is unbiased for p m .

We conclude this subsection with a brief methodological improvement of the sample concurrence probability estimator pm based on permutation bootstrap. The idea is to compute the estimator pm for several independent random permutations of the sample X 1 , . . . , X n . Then the average of the resulting estimator would have a lower variance and the same mean p m .

Formally, this procedure is justified by the following simple observation based on Rao-Blackwellization. Consider the lexicographic linear order in R k , denoted ≺, and let X (1) ≺ X (2) ≺ • • • ≺ X (n) be the sorted sample obtained from X 1 , . . . , X n . The independence of the X i 's and the continuity of their marginals entails that the above ordering is strict with probability one.

Let T {X 1 , . . . , X n } = (X (1) , . . . , X (n) ). It can be shown that T is a sufficient statistic for the parameter p m = p m (s 1 , . . . , s k ) and the Rao-Blackwell theorem implies the following propostion.

E p * m -p m 2 ≤ E ( pm -p m ) 2 . ( 29
)
Moreover we have

p * m = 1 n! σ∈S n pm {X σ(1) , . . . , X σ(n) }, (30)
and where S n denotes the set of all permutations of {1, . . . , n}. An alternative expression for p * m is

p * m = 1 n m n i=1 d i m -1 ( 31
)
where

d i = n k=1 1 {X k <X i } and d i m-1 = 0 if d i < m -1.
Proof. Let f (x; p m ), x ∈ R k , be the density of X = {X(s 1 ), . . . , X(s k )} with respect to some dominating measure λ. The fact that for the likelihood, we have

L(p m ; X i , i = 1, . . . , n) = n i=1 f (X i ; p m ) = n i=1 f (X (i) ; p m ),
shows that T = {X (1) , . . . , X (n) } is a sufficient statistic for p m . The inequality in (29) follows by appealing to the Rao-Blackwell Theorem or simply applying the conditional form of Jensen's inequality. The independence of the X i 's and the lack of ties (with probability one) in the lexicographic order imply

P{X σ(i) = X (i) , i = 1, . . . , n} = 1 n! ,
for all σ ∈ S n . This shows that p * n = E( pm | T ) is expressed as in ( 30). By the definition of the sample concurrence probability estimator, we get

p * m = 1 n! σ∈S n 1 [n/m] [n/m] r=1 1 {sample concurrence occurs within (X σ(i+(r-1)m) : i = 1, . . . , m)} = 1 n m S ∈P m (n)
1 {sample concurrence occurs within (X i : i ∈ S )}

where P m (n) is the collection of all subsets S ⊂ {1, . . . , n} of m elements. Given a subset S ∈ P m (n) and i 0 ∈ S , it is easy to see that X i 0 dominates (X i ) i∈S if and only if S \ {i 0 } is included in the set {k = 1, . . . , n; X k < X i 0 }. We deduce that for a given index i 0 , the number of subset S ∈ P m (n) such that X i 0 dominates (X i ) i∈S is equal to

d i 0
m-1 . Formula (31) follows easily.

The above result shows that the estimator p * m is superior to pm in terms of mean squared error. From a numerical point of view, formula ( 31) is much more computationally efficient than (30). We shall refer to p * m as to the sample concurrence probability bootstrap estimator. It is significantly better, in practice, than the simple sample concurrence probability estimator pm and therefore in applications we recommend using only p * m . As indicated above, in the case of bivariate concurrence, the bias of pm can be removed. As in Relation (28), we obtain the following unbiased modification of p * m for the case of pairwise concurrence

p * m = m p * m -1 m -1 . ( 32
)
The respective performances of pm , p * m and p * m are analyzed in Section 5.

Extremal concurrence probability estimators

Although the sample concurrence probability can be easily estimated, deriving an estimator for the extremal concurrence probability seems at first sight difficult since in (10) the stochastic processes Y and Ỹ are not observable. Recall that for statistical purposes we often assume that the pointwise block maxima, e.g. pointwise annual maxima, are distributed according to some max-stable process and thus we observe η but not Y. Fortunately, Theorem 3 enables us to estimate p(s 1 , . . . , s k ) without the need of observing Y.

Based on η 1 , . . . , η n independent copies of η, one possible estimator for p(s 1 , . . . , s k ) is to consider the sample counterpart of ( 14), i.e., p(s 1 , . . . ,

s k ) = k r=1 (-1) r J⊆{1,...,k} |J|=r 1 n n i=1 log        1 n n k=1 1 {η k (s j )≤η i (s j ), j∈J}        . ( 33
)
In the innermost summation, we include the index k = i to ensure that the logarithm is always well defined. Although estimator (33) seems natural it has undesirable properties since it is not linear in the data and is thus likely to show some bias for small sample sizes. Also the study of its asymptotic properties seems delicate. However, to reduce the bias of this estimator, it is always possible to use a Jackknife procedure.

Fortunately, in the bivariate case, it is possible to get an unbiased estimator based on Theorem 4. This theorem suggests the simple estimator p ≡ p(s 1 , s 2 ) = τ, where τ is the Kendall's τ statistic, i.e.,

p ≡ τ = 2 n(n -1) 1≤i< j≤n sign{η i (s 1 ) -η j (s 1 )}sign{η i (s 2 ) -η j (s 2 )}, ( 34
)
that is well known to be unbiased and such that [Dengler, 2010, Theorem 4

.3] √ n(τ -τ) -→ N(0, σ 2 τ ), σ 2 τ = 15Var[F{η(s 1 ), η(s 2 )} -F{η(s 1 )} -F{η(s 2 )}].
Although the asymptotic variance σ 2 τ is hard to evaluate as it requires knowledge of the dependence structure, in practice it can be accurately and consistently estimated using Jacknife [START_REF] Schemper | Nonparametric estimation of variance, skewness and kurtosis of the distribution of a statistic by jackknife and bootstrap techniques[END_REF].

Integrated concurrence probabilities and area of concurrence cell

As we will see in Section 6, the above methodology can be used to provide bivariate concurrence probability maps s → p(s 0 , s) centered at a given location s 0 . Such maps show how fast the dependence in extremes decreases when moving away from s 0 . A drawback of this approach is that one may produce one such map for every choice of an origin s 0 and the choice of an origin is hence quite arbitrary. To bypass this issue, we propose to consider the integrated concurrence probability

I(s 0 ) = s∈X p(s 0 , s)ds, s 0 ∈ X.
Intuitively, this quantity measures how fast the dependence in extremes decreases when moving away from s 0 . Interestingly, it can be related to the notion of concurrence cell and its area. Consider the Poisson process representation of the max-stable random field η = {η(s) : s ∈ X} in (5). Recall that we have a concurrence of extremes at sites s 0 and s if η(s 0 ) = φ(s 0 ) and η(s) = φ(s), for the same φ ∈ Φ. Let C(s 0 ) denotes the random set of all sites s that are in a concurrence relation with s 0 . This set will be referred to as the concurrence cell containing the site s 0 .

Proposition 5. For any site s 0 ∈ X, we have I(s 0 ) = E{|C(s 0 )|} where |C(s 0 )| is the d-dimensional volume of C(s 0 ).

Proof. Observe that the concurrence probability satisfies p(s 0 , s) = E 1 C(s 0 ) (s) and that the volume of the concurrence cell is given by

|C(s 0 )| = X 1 C(s 0 ) (s)ds.
The result follows by applying the Tonelli-Fubini's theorem. recommend using slightly larger values of m so as to ensure that the block-maxima are closer to a max-stable model but also to take into account that data usually exhibit serial dependence, e.g., daily observations. We now investigate the performance of the extremal concurrence estimator p. Figure 2 shows the evolution of the root mean squared error as the number of spectral functions n 0 in (35) and the theoretical extremal concurrence probability increase. As expected, as the sample size n grows, the estimator p becomes much more efficient. Interestingly, for small sample sizes, p appears to be fairly robust to the lack of max-stability in the data, i.e., n 0 < ∞. This is not true anymore for larger sample sizes since, as expected, p becomes increasingly more efficient as the number of spectral functions increases.

Finally, we compare the performance of the sample concurrence probability estimators p * m and p * m with their extremal concurrence counterpart p. To compare the two types of estimators on a fair basis, we analyze their behaviour when the simulated data are either perfectly max-stable or non max-stable, but in the domain of attraction of a maxstable distribution.

Figure 3 shows boxplots of the sample p * m , unbiased sample p * m and extremal concurrence probability estimators p, based on 2000 Monte-Carlo realizations of both a Brown-Resnick and extremal-t models. Recall that we focus here on pairwise concurrence probabilities. In this case, the extremal concurrence probability coincides with Kendall's τ and therefore, the estimator p in (34) is in fact unbiased for the case of max-stable data. This is confirmed by the results in Figure 3. As expected, the variability of all estimators decreases as the sample size grows; the extremal concurrence probability estimator being the most precise one. Since the simulated data are max-stable, we can see that the sample concurrence probability estimator is biased even when the sample size is large while the remaining two estimators are, as expected, unbiased. Overall the extremal concurrence probability appears to be the best estimator provided that the data are max-stable.

To corroborate this finding, Table 1 reports Monte-Carlo sample means and standard deviations of these estimators as the assumption of max-stability becomes more accurate, i.e., as the number n 0 of spectral functions in (35) grows. As expected, when the max-stability assumption is most unreasonable, i.e., n 0 = 1, all estimators show a substantial bias with the extremal concurrence probability estimator p having the largest bias while the unbiased sample one p * m the lowest. As the assumption of max-stability becomes increasingly more accurate, the bias of the unbiased sample concurrence and extremal concurrence estimators improve. When this assumption holds exactly (indicated by n 0 = ∞), both of these estimators exhibit essentially no bias as stipulated by the theory and seen in Figure 3. The sample concurrence probability estimator appears to be biased in all situations-the bias being less significant as the number of spectral functions is larger. Interestingly, whatever the estimator considered, the bias and variance appear to increase as the theoretical extremal concurrence probability value p becomes smaller. Overall the extremal concurrence probability estimator p in (34) has the lowest variability. 6 Concurrence of temperature extremes in continental USA

In this section, we apply the developed methodology to estimate the probabilities of concurrence associated with extreme temperatures-both extreme cold and hot events. The data consists of daily temperature minima and maxima recorded at 424 weather stations over the period 1911-2010. The spatial distribution of these stations is given in Figure 4. This data set, as a subset of the United States Historical Climatological Network USHCN [2014], was chosen as it meets very high data quality standards and involves fewer than 2.4% missing values while spanning the entire territory of continental US. It can be freely downloaded from http://cdiac.ornl.gov/.

To avoid any seasonal influence on our results we decided to analyze minima and maxima for each season sepa- rately. We focus on the concurrence of extreme cold (minima) during the Fall and Winter seasons-generally colorcoded in blue; and extreme hot (maxima) during the Spring and Summer seasons-generally color-coded in red. The right panel of Figure 4 shows the times series of these seasonal extrema for four selected weather stations. These stations were selected as they recorded the top seasonal records over the whole spatial and temporal domains. We can see that all four time series of seasonal extremes (cold in blue and hot in red) at these stations appear to be stationary without any clear temporal trend. This is in contrast with the generally accepted trend of about 0.2 • C per decade for average temperatures [START_REF] Stocker | Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change[END_REF].

Figure 5 plots the estimated spatial distribution of the extremal concurrence probabilities function for each season, relative to the chosen station. More precisely, for a given origin location s 0 , the maps display estimates of the pairwise concurrence probability p(s 0 , s) as a function of s. These maps were obtained by first computing the estimator (34) over all 423 pairs of stations (s 0 , s) and then interpolated using thin plate splines (on logit scale) provided by the R package fields [START_REF] Nychka | fields: Tools for spatial data[END_REF]. As expected, the highest concurrence probability occurs in the neighbourhood of the selected stations independently of the season. The areal extent of high concurrence probabilities, however, seem to be larger for minimum temperatures (cold extremes) than for maximum temperatures (hot extremes). This finding is consistent with the physical notion of entropy, i.e., when the ambient temperature is higher (Spring and Summer seasons), the entropy is greater and hence involves less spatial dependence than for cooler temperatures leading to smaller probability of simultaneous extremes. This difference can be also attributed to the fact that extreme cold temperatures are often due to high-pressure systems, which tend to linger longer and cover a larger spatial area than warm fronts giving rise to concurrence of extreme hot events.

Although Figure 5 displays interesting patterns, it has the drawback of being dependent on the choice of the origin, i.e., the selected station. As stated in Section 5, it is possible to bypass this hurdle by considering the area of concurrence cell. Figure 6 plots the estimated spatial distribution of the concurrence cell area for the preindustrial period, i.e., 1910-1975, and the postindustrial one, i.e., 1976-2010. To emphasize the possible impact of anthropogenic influences, the ratio of these two cell areas is also reported. We can see that during the last sixty years the expected cell area for winter minima have increased of about 30% over the whole USA while there is a decrease of about the same amount for summer maxima. These findings indicates that today's climate shows cold spells that have a larger impact than in the beginning of the 20th century while hot spells are more localized. Our results agree with the conclusions drawn by [START_REF] Field | IPCC, 2012: Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation. A Special Report of Working Groups I and II of the Intergovernmental Panel on Climate Change[END_REF] who states that "there is evidence from observations gathered since 1950 of change in some extremes". These changes in the concurrence patterns of summer extremes can be attributed to global warming since an increase in entropy generally leads to more "mixing" in the system and hence less dependence leading to smaller areas of concurrence. The changes in concurrence patterns of extreme cold events, however, are harder to explain. They may be triggered by structural changes in important climatological mechanisms such as the Arctic : 1910-1950, 1951-2010, and their ratio (1951-2010 Oscillation. Finally, we consider another cut of the data by stratifying according to an important climate phenomenon known as the El Niño Southern Oscillation (ENSO). Positive ENSO (El Niño) refers to the event of a warm-up of the surface water in the central and east-central equatorial Pacific ocean. It is well known that years with high ENSO have a general warming effect in North America during the winter season. The opposite effect of negative ENSO (La Niña) is characterized by a cool-down in the same area of the Pacific and it generally leads to unusually cold winters in the northwestern part of the US, northern California and the north-central states [START_REF] Graham | NASA: Earth Observatory[END_REF]. Figure 7 shows estimates of concurrence cell areas anomalies for winter minima and summer maxima. These anomalies were defined as pointwise deviations from the expected cell area obtained from La Nada seasons, i.e., neither El Niño nor La Niña seasons. We can see that La Niña does not seem to have an impact on the spatial coverage of winter minima but that El Niño seems to induce more massive cold extremes over the whole USA. For the summer season, La Niña seems to reduce the spatial extent of heat waves over the whole USA while El Niño has a less pronounced impact-although it generally yields larger spatial coverages especially along the East coast.

Discussion

In this paper we introduced a new framework for the analysis of dependence of extremes: the extremal/sample concurrence probability. This tool plays a similar role to that of the extremal coefficient but has the benefit, as a probability, of being more interpretable and intuitive. Theoretical properties and closed forms of these concurrence probabilities have been established and several estimators have been proposed. A simulation study has shown that the proposed estimators work well in practice and that they give a new insight about the dependence of extremes, such as the spatial distribution of the expected concurrence cell area of extreme temperature in the continental US.

with Φ 1 , . . . , Φ m independent copies of Φ. By max-stability, η = m -1 max 1≤i≤m Z i has the same distribution as η and Φ = ∪ 1≤i≤m {m -1 φ, φ ∈ Φ i } has the same distribution as Φ. We consider the following events A = {sample concurrence occurs for Z 1 , . . . , Z m } B = {extremal concurrence occurs for η} Clearly, P(A) = p m and P(A) = P(A ∩ B) + P(A ∩ B c ) with B c the complementary set of B. We analyze the two terms separately.

Observe first that if extremal concurrence occurs then we also have sample concurrence. Indeed, if B occurs, then one function of Φ dominates all the others at (s 1 , . . . , s k ). This function is of the form m -1 φ with φ ∈ Φ i , for some i = 1, . . . , m, showing that Z i dominates Z 1 , . . . , Z m , i.e., we have sample concurrence. Hence B ⊂ A and P(A ∩ B) = P(B) = p. We now consider the second term and the event A ∩ B c , i.e., sample concurrence occurs in Z 1 , . . . , Z m but not extremal concurrence for η. Let π be the hitting scenario of η. We know that B c is equivalent to |π| ≥ 2, i.e. the maximum at locations (s 1 , . . . , s k ) is attained by at least two functions in Φ. These functions are of the form m -1 φ j , 1 ≤ j ≤ , with = |π| and φ j ∈ Φ i j for some 1 ≤ i j ≤ m. If A is also realized, i.e., some Z i dominates Z 1 , . . . , Z m , then we must have i 1 = • • • = i = i. Note, however, that since the point processes Φ i , i = 1, . . . , m are independent and identically distributed, any given function m -1 φ ∈ Φ ≡ ∪ m i=1 m -1 Φ i , independently from the others, has equal chance of coming from any one of the m point processes m -1 Φ i , i = 1, . . . , m. Therefore, the probability that all functions contributing to the maximum at sites (s 1 , . . . , s k ) are assigned to component i is m -. Since there are m possible choices for the index i, we deduce

P(A ∩ B c ) = k =2 P(|π| = )m 1-.
Equation (36) follows.

A.3 Proofs of Theorems 5 and 6

In the context of these two theorems, we have

F n (x) = P B -1/2 n {S n -E(S n )} ≤ x ≡ P          [n/m] p m (1 -p m ) ( pm -p m ) ≤ x          , x ∈ R.
The Berry-Essen theorem (see e.g. Theorem V.2.3 in [START_REF] Petrov | Sums of independent random variables[END_REF]) implies that

sup x∈R |F n (x) -Φ(x)| ≤ AL n , ( 38
)
where A is an absolute constant, Φ(x) denotes the standard Normal cumulative distribution function and (39)

Proposition 4 .

 4 For p * m = E( pm | T ) we have E p * m = p m and

Figure 2 :

 2 Figure 2: Evolution of the root mean squared error for p as the theoretical extremal concurrence probability p and the number of spectral function n 0 in (35) increase. These estimates were obtained from 2000 Monte-Carlo samples of size n with, from left to right, n = 25, 50, 100, 500.

Figure 3 :

 3 Figure 3: Boxplots of the sample (red / middle), unbiased sample (green / right) and extremal (grey / left) concurrence probability estimators at distance lags h = 1, 2, 3, 4. The boxplots were obtained from 2000 independent estimates. From left to right: the sample size is respectively 25, 50, 100 and 500. The top panel corresponds to an extremal-t model with ν = 5, and correlation function ρ(h) = exp(-h/10). The bottom panel corresponds to a Brown-Resnick model with semi variogram γ(h) = h/3. For each panel, the solid line represents the corresponding theoretical extremal concurrence probability function.

Figure 4 :

 4 Figure 4: Left: Spatial distribution of the 424 weather stations. The triangles indicate the selected stations for the analysis-upward: daily maxima, downward: daily minima. Right: The seasonal extrema time series of the selected stations.

Figure 5 :Figure 6 :

 56 Figure 5: Maps of the extremal concurrence probability for the four selected stations. Top left: Fall (September, October November), top right: Winter (December, January, February), bottom left: Spring (March, April, May) and bottom right: Summer (June, July, August).

Figure 7 :

 7 Figure7: Spatial distribution of the estimated concurrence cell areas anomalies in squared degree, i.e., around 1000 km 2 for winter minima (top) and summer maxima (bottom). The data were stratified into three classes: El Niño, La Niña and the base class "La Nada". The left panels show the anomalies for El Niño, the right ones La Niña.

  p mp ∼ c r m -r , m → ∞,(37)for some r ∈ {1, . . . , k -1} and c r > 0-cf. Proposition 3. Let S n = [n/m] pm = [n/m] i=1 ξ i,n , where ξ i,n are iid Bernoulli(p m ). Proof of Theorem 5. Let B n = Var(S n ) = [n/m]p m (1p m ) and introduce the cumulative distribution function

E

  |ξ i,np m | 3 .Using that E(|ξ i,np m | 3 ) ≤ p m (1p m ) and straightforward algebra, we obtainL n ≤ (n/m)p m (1p m ) -1/2 .

Table 1 :

 1 Performance of the sample ( p * m ), unbiased sample ( p * m ) and extremal ( p) concurrence probability estimators. The table report the sample mean and the standard deviation in paren based on 2000 Monte-Carlo replicates. The data are either simulated from an extremal-t model with correlation function ρ(h) = exp(-h/10) and ν = 5 degrees of freedom or from its truncated representation with n 0 extremal functions. Throughout this simulation study the block size is held fixed to m = 10, independently of the sample size n.

			p = 0.25			p = 0.50			p = 0.75
		p * m	p * m	p	p * m	p * m	p	p * m	p * m	p
					Sample size n = 20
	n 0 = 1	0.41 (0.24)	0.35 (0.26)	0.47 (0.13)	0.64 (0.22)	0.60 (0.24)	0.71 (0.09)	0.83 (0.15)	0.81 (0.17)	0.87 (0.05)
	n 0 = 10	0.34 (0.24)	0.27 (0.25)	0.31 (0.14)	0.57 (0.23)	0.52 (0.26)	0.58 (0.12)	0.79 (0.17)	0.77 (0.19)	0.80 (0.07)
	n 0 = 15	0.33 (0.24)	0.27 (0.25)	0.30 (0.15)	0.56 (0.23)	0.52 (0.26)	0.56 (0.12)	0.78 (0.17)	0.76 (0.19)	0.78 (0.07)
	n 0 = ∞	0.33 (0.24)	0.27 (0.25)	0.25 (0.15)	0.55 (0.24)	0.50 (0.26)	0.50 (0.13)	0.77 (0.18)	0.75 (0.20)	0.75 (0.08)
					Sample size n = 50
	n 0 = 1	0.41 (0.13)	0.35 (0.14)	0.47 (0.08)	0.65 (0.10)	0.61 (0.12)	0.71 (0.05)	0.84 (0.07)	0.82 (0.07)	0.87 (0.03)
	n 0 = 10	0.34 (0.13)	0.26 (0.14)	0.31 (0.09)	0.57 (0.12)	0.52 (0.13)	0.57 (0.07)	0.79 (0.08)	0.76 (0.09)	0.80 (0.04)
	n 0 = 15	0.33 (0.13)	0.25 (0.14)	0.29 (0.09)	0.56 (0.12)	0.51 (0.13)	0.56 (0.07)	0.78 (0.08)	0.76 (0.09)	0.79 (0.04)
	n 0 = ∞	0.32 (0.13)	0.25 (0.14)	0.24 (0.09)	0.54 (0.12)	0.49 (0.14)	0.50 (0.08)	0.77 (0.09)	0.74 (0.09)	0.75 (0.05)
					Sample size n = 100
	n 0 = 1	0.41 (0.08)	0.35 (0.09)	0.46 (0.06)	0.65 (0.07)	0.61 (0.07)	0.71 (0.03)	0.83 (0.04)	0.82 (0.04)	0.87 (0.02)
	n 0 = 10	0.34 (0.09)	0.26 (0.10)	0.31 (0.06)	0.57 (0.08)	0.52 (0.09)	0.57 (0.05)	0.78 (0.05)	0.76 (0.05)	0.80 (0.03)
	n 0 = 15	0.33 (0.09)	0.26 (0.10)	0.29 (0.06)	0.56 (0.08)	0.51 (0.09)	0.55 (0.05)	0.78 (0.05)	0.76 (0.06)	0.78 (0.03)
	n 0 = ∞	0.33 (0.09)	0.25 (0.10)	0.25 (0.07)	0.55 (0.08)	0.50 (0.09)	0.50 (0.05)	0.78 (0.05)	0.75 (0.06)	0.75 (0.03)

Acknowledgements

M. Ribatet was partly funded by the MIRACCLE-GICC and McSim ANR projects. The authors gratefully acknowledge the help of Prof. Paul H. Whitfield with the interpretation of the results on concurrence for temperature extremes and also for the suggestion to stratify by El Niño/La Niña effect.

m (right) as the block size m and the sample size n increase. These estimates were obtained from 2000 Monte-Carlo samples sampled from a Brown-Resnick model with semivariogram γ(h) = h/1.627. This semivariogram was chosen such that the theoretical extremal concurrence probability is p = 0.5. The red circles indicate the optimal block sizes as defined by (27) and their corresponding optimal root mean squared error (26).

We will provide and discuss in Section 6 some maps of the integrated concurrence probability s 0 → I(s 0 ) that allow to evaluate at each location s 0 the dependence in extremes around s 0 . For a detailed study of the properties of the concurrence cells associated to a max-stable random field and of the tessellation of the entire domain generated by the concurrence cells, please refer to the recent work of [START_REF] Dombry | Random tessellations associated with max-stable random fields[END_REF].

Simulation study

In this section, we analyze the performance of the pairwise sample concurrence probability estimators pm , p * m and p * m defined in (25), ( 31) and (32) respectively, and that of their extremal counterpart p in (34). Since the latter estimator relies on the max-stability assumption while the former three assume that observations belong to the max-domain of attraction, we need to handle both situations. The first one is well known and consists in sampling from max-stable processes using the methodology of [START_REF] Schlather | Models for stationary max-stable random fields[END_REF]. In the second situation, to be able to control the degree to which the model differs from a max-stable one, we consider the following partial maxima

where Y i are as in (4), U 1 , . . . , U n 0 independent U(0, 1) random variables and for some suitable n 0 ∈ N. By construction, η belongs to the max-domain of attraction of η in (4) and in some sense can be viewed as a truncation of the spectral representation in (4) (see, e.g., the proof of Proposition 3.1 in [START_REF] Stoev | Extremal stochastic integrals: a parallel between max-stable processes and α-stable processes[END_REF].)

We first focus only on the sample concurrence probability estimators, i.e., pm and p * m , and analyze their performance with respect to the block size m and the sample size n. Based on a Monte-Carlo simulation, Figure 1 shows the evolution of the root mean squared error as the block size grows. As expected, both estimators become increasingly more efficient as the sample size grows and, as seen from ( 29), the permutation estimator p * m is more efficient than pm -independently of the block size m and the sample size n. The circles on the plot indicate the asymptotically optimal block size in ( 27), which are valid olnly for max-stable data. As expected the observed optimal block sizes are in good agreement with the theoretical ones. In practice, however, since the data are not exactly max-stable, we

A Proofs

A.1 Proof of the continuity of Π in Theorem 1

has finitely many points ψ 1 , . . . , ψ N ε and no point of Ψ lies on the boundary ∂K ε . The convergence Ψ n → Ψ entails that for n large enough, Ψ n ∩ K ε has the same number of points ψ n,1 , . . . , ψ n,N ε that can be reordered in such a way that ψ n,i → ψ i as n → ∞, i = 1, . . . , N ε . Assume that the maxima max 1≤i≤N ε ψ i (s j ), j = 1, . . . , k, are uniquely attained. Then the hitting scenario Π(Ψ) is well defined and depends only on {ψ i : i = 1, . . . , N ε }. By the convergence ψ n,i → ψ, for large n the maxima max 1≤i≤N ε ψ n,i (s j ) are uniquely attained so that the hitting scenario Π(Ψ n ) is well defined and depends only on {ψ n,i : i = 1, . . . N ε }. It is not difficult to see (although tedious to write formally) that the convergence {ψ n,i : i = 1, . . . , N ε } → {ψ i : i = 1, . . . , N ε } implies Π(Ψ n ) = Π(Ψ) for large n. This proves the announced continuity for the mapping Π.

A.2 Proof of Proposition 3

We shall prove below the following formula, which may be of independent interest. Lemma 2. In the context of Proposition 3, we have

where π is the extremal hitting scenario and |π| its number of components.

Proposition 3 follows directly from (36). Indeed, the monotonicity of p m is immediate and since p = P(|π| = 1), we have

If p < 1, then at least one of the probabilities P(|π| = ), = 2, . . . , k is non-zero and by (36) the asymptotic equivalence (p mp) ∼ c r /m r , m → ∞ holds where r ∈ {1, . . . , k -1} is the smallest integer, such that c r = P(|π| = r + 1) > 0.

Proof of Lemma 2. Since Z 1 , . . . , Z m are independent with the same distribution as η, we can suppose from (5) that

By using (37) and Slutsky's theorem, we obtain the final result.

Proof of Theorem 6. In case λ < ∞, since p = 0, from (37), we have p m ∼ c r /m r , n → ∞, and hence

which completes the proof.