Mikael Escobar-Bach

Yuri Goegebeur
email: yuri.goegebeur@imada.sdu.dk

Armelle Guillou

Alexandre You

Bias-corrected and robust estimation of the bivariate stable tail dependence function

Keywords: Multivariate extreme value statistics, stable tail dependence function, robustness, bias-correction AMS Subject classification: 62G05, 62G20, 62G32

We consider the estimation of the bivariate stable tail dependence function and propose a bias-corrected and robust estimator. We establish its asymptotic behavior under suitable assumptions. The finite sample performance of the proposed estimator is examined on a simulation study involving both uncontaminated and contaminated samples.

Introduction

Modelling dependence among variables is a challenging topic in multivariate extreme value theory. For instance, in the case of environmental data, it can be of interest to link still water levels and wave heights in order to avoid flooding, whereas in finance measuring the dependence between risky asset returns can be crucial. Extremal dependence can be measured by coefficients that give a representative picture of the dependence, like e.g. the coefficient of tail dependence.

Alternatively one could use functions that give a complete characterisation of the extreme dependence, e.g. the spectral distribution function or the Pickands dependence function. We refer to [START_REF] Beirlant | Statistics of Extremes -Theory and Applications[END_REF] and de [START_REF] De Haan | Extreme Value Theory -An Introduction[END_REF] for recent accounts on the available approaches for describing tail dependence. In this paper, we focus on the bivariate stable tail dependence function, firstly introduced by Huang (1992) as follows. Let (X, Y ) be a bivariate vector with continuous marginal distribution functions F X and F Y . The stable tail dependence function is defined for (x, y) ∈ R 2 + as L(x, y) := lim

t↓0 1 t P (1 -F X (X) ≤ tx or 1 -F Y (Y ) ≤ ty)
which can be rewritten as

L(x, y) = lim t↓0 1 t [1 -F (F ← X (1 -tx), F ← Y (1 -ty))] (1)
with F the bivariate distribution function of the vector (X, Y

L m (x, y) = 1 m n i=1 1l {X i ≥X n-[mx]+1,n or Y i ≥Y n-[my]+1,n } . (2)
The weak convergence of this estimator, after proper normalisation, was established by [START_REF] Huang | Statistics of bivariate extremes[END_REF]. Other references include [START_REF] Qi | Almost sure convergence of the stable tail empirical dependence function in multivariate extreme statistics[END_REF], [START_REF] Drees | Best attainable rates of convergence for estimators of the stable tail dependence function[END_REF], [START_REF] Schmidt | Non-parametric estimation of tail dependence[END_REF]. This estimator can suffer from bias as illustrated in Peng (2010), [START_REF] Fougères | Bias correction in multivariate extremes[END_REF] or [START_REF] Beirlant | Bias-corrected estimation of stable tail dependence function[END_REF] where bias-corrected estimators have been proposed. However, in practical data analysis, it can also happen that observations are contaminated in the sense that some outliers may have a disturbing effect on the estimators. Such a situation is illustrated in Figure 1 where we plot the mean and the mean squared error (MSE) of L m (x, y) and some recently introduced bias-corrected estimators for L(x, y) as a function of m. The estimates are based on 500 samples of size 1000 from a Logistic distribution with 5% of contamination (as described in Section 4). As is clear from this figure, the classical estimators behave very poorly. On the left: mean and on the right: MSE based on 500 samples of size 1000. The horizontal line in the left display is the true value of L(x, y).

To solve this issue, it is important to propose robust estimators which also keep the nice property of being asymptotically unbiased as those proposed by [START_REF] Fougères | Bias correction in multivariate extremes[END_REF] and [START_REF] Beirlant | Bias-corrected estimation of stable tail dependence function[END_REF]. This is the aim of the present paper, where we propose to use the minimum density power divergence (MDPD) criterion introduced by [START_REF] Basu | Robust and efficient estimation by minimizing a density power divergence[END_REF]. This idea consists of defining the density power divergence between density functions f and h as follows

∆ α (f, h) :=              R h 1+α (y) -1 + 1 α h α (y)f (y) + 1 α f 1+α (y) dy, α > 0, R log f (y) h(y) f (y)dy, α = 0.
Remark that for α = 0 one recovers the Kullback-Leibler divergence, whereas setting α = 1 leads to the L 2 divergence. Assume that the density function h depends on a parameter vector θ, and let Y 1 , . . . , Y n be a sample of independent and identically distributed random variables according to the density function f . The minimum density power divergence estimator (MDPDE) is the point θ minimizing the empirical density power divergence

∆ α (θ) :=                  R h 1+α (y)dy -1 + 1 α 1 n n i=1 h α (Y i ), α > 0, - 1 n n i=1 log h(Y i ), α = 0.
Note that in case α = 0, the empirical density power divergence corresponds with minus the log-likelihood function. The parameter α controls the trade-off between efficiency and robustness of the MDPDE: it becomes more efficient but less robust against outliers as α gets closer to zero, whereas by increasing α the robustness increases and the efficiency decreases.

In the present paper we introduce a robust and bias-corrected estimator of the stable tail dependence function L. To the best of our knowledge, robust and bias-corrected estimation of L has not been considered in the extreme value literature. The remainder of the paper is organized as follows. In Section 2, we introduce our estimator of L(x, y), while in Section 3 we establish its asymptotic properties. The finite sample performance of our procedure is examined in Section 4 on a simulation study. All the proofs of our results are postponed to the appendix.

Construction of the estimators

Let (X, Y ) be a bivariate random vector with continuous marginal distribution functions F X and F Y satisfying

P (1 -F X (X) < x, 1 -F Y (Y ) < y) = x d 1 y d 2 g(x, y) 1 + 1 η δ(x, y) , (x, y) ∈ [0, 1] 2 \ {(0, 0)}, (3)
where

d 1 , d 2 are positive constants such that d 1 +d 2 = η -1 ∈ (0,
for all (x, y) ∈ [0, ∞) 2 \ {(0, 0)}. Also, ξ is assumed to be continuous, homogeneous of order τ > 0, and the convergence is uniform on {(x, y) ∈ [0, ∞) 2 |x 2 + y 2 = 1}. Note that model ( 3) is essentially a condition on the copula function C of F . Indeed, one easily verifies that

P (1 -F X (X) < x, 1 -F Y (Y ) < y) = x + y -1 + C(1 -x, 1 -y).
If η ∈ (0, 1) then we are in the asymptotic independence case, and in that context L(x, y) = x+y for x, y > 0. As such, estimating L is not relevant, and main interest is in estimating η, as has been done in [START_REF] Beirlant | Some comments on the estimation of a dependence index in bivariate extreme value statistics[END_REF] For convenience we suppose that the marginal distributions are unit Pareto. In that case, model

(3) with η = 1 becomes

P (X > x, Y > y) = x -d 1 y -d 2 g 1 x , 1 y 1 + δ 1 x , 1 y , x, y > 1. ( 5
)
Note that one can write

P (X > x, Y > y) = P X > x, ω 1 -ω Y > x ,
where ω := x/(x + y) ∈ (0, 1) can be interpreted as being a radial parameter. Thus the transformed variable Z ω := min(X, ω 1-ω Y ) admits the following survival function

F Zω (x) = 1 -ω ω -d 2 g 1 x , ω 1 -ω 1 x x -1 1 + δ 1 x , ω 1 -ω 1 x = 1 -ω ω -d 2 g 1, ω 1 -ω x -1 1 + δ 1 x , ω 1 -ω 1 x =: K ω x -1 (1 + δ ω (x)) (6)
using the homogeneity of order 0 of g. Now, if we come back to our initial problem which is estimating the bivariate stable tail dependence function, recall that

L(x, y) = lim t↓0 1 t P (1 -F X (X) ≤ tx) + 1 t P (1 -F Y (Y ) ≤ ty) - 1 t P (1 -F X (X) ≤ tx, 1 -F Y (Y ) ≤ ty) = lim t↓0 x + y - 1 t P Z 1-ω ≥ 1 tx = lim t↓0 x + y -xK 1-ω 1 + δ 1-ω 1 tx
by [START_REF] Beirlant | Second-order refined peaks-over-threshold modelling for heavy-tailed distributions[END_REF]. This leads to the simple estimator

L m (x, y) = x + y -x K 1-ω,m , (7)
which requires an estimate for K 1-ω . This estimator K 1-ω,m can now be obtained as follows.

Consider the tail distribution of Z 1-ω given in [START_REF] Beirlant | Second-order refined peaks-over-threshold modelling for heavy-tailed distributions[END_REF].

H(y; δ 1-ω (u), τ ) :=    1 -y -1 1 + δ 1-ω (u)(1 -y -τ ) -1 , y > 1, 0, y ≤ 1,
where δ 1-ω (u) > max{- In practice, the sample of pairs at our disposal, say (X 1 , Y 1 ), ..., (X n , Y n ), have unknown margins.

Thus the first step is to transform the margins into unit Pareto distributions. To this aim, we can use the empirical distribution functions of the X and Y observations. This gives

Z 1-ω,i := min n + 1 n + 1 -R X i , 1 -ω ω n + 1 n + 1 -R Y i with R X
i and R Y i denoting the rank of X i and Y i , i = 1, ..., n, in the respective samples. Then, in a second step, the parameter δ 1-ω (u) of the EPD is estimated by fitting the density function h to the relative excesses E j := Z 1-ω,n-j+1,n / Z 1-ω,n-m,n , j = 1, ..., m, where 1 ≤ m ≤ n -1 using the minimum density power divergence criterion. Note that, as usual in extreme value statistics, we use for the threshold u the n-m largest observation Z 1-ω,n-m,n within the sample.

Concerning the parameter τ of the EPD, we fix it at some value, either the true one or a misspecified one, or we estimate it externally in a consistent way. External estimation of second order parameters is quite common in the univariate framework, see e.g. [START_REF] Gomes | Asymptotically unbiased" estimators of the tail index based on external estimation of the second order parameter[END_REF], [START_REF] Gomes | Tail index estimation for heavy-tailed models: accommodation of bias in weighted log-excesses[END_REF] and [START_REF] Dierckx | An asymptotically unbiased minimum density power divergence estimator for the Pareto-tail index[END_REF]. Joint estimation of δ 1-ω (u) and τ has been considered, but leads to theoretical difficulties. In particular, we would need to introduce √ mq 1 (k/n) → ∞, which is not in accordance with the assumptions of our Theorem 1 where the quantile process of Z 1-ω is studied. To get around this one would need to introduce a third order model with strict unit Pareto margins for deriving the asymptotic properties. The fact that we cannot use the empirical transformation of the margins (as we do in our methodology) due to the inversion procedure in the proof of Theorem 1, implies that such a method with strict Pareto margins would be very restrictive and without practical interest. Moreover, from a practical perspective, one can expect that joint estimation of δ 1-ω and τ will lead to an increase in the asymptotic variance, as was already observed in the univariate context. Therefore, this approach is not pursued in this paper. Finally, in a last step, an estimator of K 1-ω follows

K 1-ω,m = m n Z 1-ω,n-m,n 1 + δ 1-ω ( Z 1-ω,n-m,n ) (8)
from which the estimator of L(x, y) can be deduced by applying [START_REF] Beirlant | Some comments on the estimation of a dependence index in bivariate extreme value statistics[END_REF]. The estimator (8) can be intuitively motivated as follows. Consider the tail distribution of Z 1-ω /u given that Z 1-ω > u:

F Z 1-ω (ux) F Z 1-ω (u) = K 1-ω (ux) -1 (1 + δ 1-ω (ux)) F Z 1-ω (u) .
Take u = Z 1-ω,n-m,n , replace the denominator of the above display by its empirical value, being m/n, and set x = 1, to obtain

1 ≈ n m K 1-ω Z -1 1-ω,n-m,n (1 + δ 1-ω (Z 1-ω,n-m,n )),
and thus

K 1-ω ≈ m n Z 1-ω,n-m,n 1 + δ 1-ω (Z 1-ω,n-m,n ) .
The estimator ( 8) is then obtained by replacing Z 1-ω,n-m,n by its empirical analogue Z 1-ω,n-m,n , and

δ 1-ω ( Z 1-ω,n-m,n ) by its MDPDE, denoted δ 1-ω ( Z 1-ω,n-m,n ).

Asymptotic properties

The crucial point in the methodology consists of computing the estimate for δ 1-ω ( Z 1-ω,n-m,n ) by using the MDPD criterion. In the sequel, although this parameter depends on the threshold Z 1-ω,n-m,n , we do not make this dependence explicit but prefer to use the notation δ 1-ω for simplicity. As we already mentioned, the parameter τ of the EPD is not estimated jointly with δ 1-ω by the MDPD method but it is fixed at some value or estimated externally. The MDPDE of δ 1-ω satisfies thus the estimating equation

∞ 1 h α (y; δ 1-ω , τ ) ∂h(y; δ 1-ω , τ ) ∂δ 1-ω dy - 1 m m j=1 h α-1 (E j ; δ 1-ω , τ ) ∂h(E j ; δ 1-ω , τ ) ∂δ 1-ω = 0, (9)
where the density function h of the EPD is given by

h(y; δ 1-ω , τ ) = y -2 1 + δ 1-ω (1 -y -τ ) -2 1 + δ 1-ω (1 -(1 -τ )y -τ ) , y > 1,
where τ > 0 and δ 1-ω > max{-1, -1/τ }. Note that the estimating equation ( 9) depends only on the data through statistics of the form (1/m) m j=1 E s j for s < 0. Thus, considering the tail quantile process Q n (t) := Z 1-ω,n-[mt],n , 0 < t < n/m, these statistics can be expressed as the

lim t↓0 q 1 (t) -1 P (1 -F X (X) < tx, 1 -F Y (Y ) < ty) P (1 -F X (X) < t, 1 -F Y (Y ) < t) -c(x, y) =: c 1 (x, y) (
exists for all x ≥ 0, y ≥ 0 with x + y > 0, a function q 1 tending to zero as t ↓ 0, and c 1 a function neither constant nor a multiple of c. Moreover, we assume that the convergence

is uniform on {(x, y) ∈ [0, ∞) 2 |x 2 + y 2 = 1}.
In case η = 1, the specific functions of Condition (SO) are for model [START_REF] Beirlant | Bias-reduced estimators for bivariate tail modelling[END_REF] given by As a preliminary result, we establish in our first theorem the behavior of the tail quantile process Q n (t), when correctly normalized. To this aim, set

           c(x, y) = x d 1 y d 2 g(x, y) g(1, 1) , q 1 (t) = τ δ(t, t), c 1 (x, y) = c(x, y) ξ(x, y) -1 τ , (11)
c 1-ω := c(1, (1 -ω)/ω), ξ 1-ω := ξ(1, (1 -ω)/ω), q(t) := P(1 -F X (X) < t, 1 -F Y (Y ) < t), let k/n = q ← (m/n),

and denote

= lim n→∞ m/k. Remark that this limit exists, see (26) infra.

Theorem 1. Let (X 1 , Y 1 ), ..., (X n , Y n ) be independent copies of the random vector (X, Y ) which has a joint distribution satisfying (3) with η = 1 such that the function c given in [START_REF] Dutang | Robust and bias-corrected estimation of the coefficient of tail dependence[END_REF] has

continuous first order partial derivatives. For m, k → ∞, as n → ∞ such that √ mq 1 (k/n) →
λ ∈ R we have that there exist suitable versions of Q n such that for all t 0 , ε > 0 sup

0<t≤t 0 t 3 2 +ε √ m k n Q n (t) - c 1-ω t - c 1-ω t 2 W t c 1-ω -λ c 1-ω t ( t c 1-ω ) τ ξ 1-ω -1 τ = o P (1),
where

W (x) := W x, x 1 -ω ω with W (x, y) := 1 √ [W (x, 0) + W (0, y) -W (x, y)] - √ [c x (x, y)W (x, 0) + c y (x, y)W (0, y)] ,
where c x and c y denote the derivatives of c with respect to x and y, respectively, and W a Gaussian process with mean zero and covariance structure given by

E (W (x 1 , y 1 )W (x 2 , y 2 )) = x 1 ∧ x 2 + y 1 ∧ y 2 -[c(x 1 , y 1 ) + c(x 2 , y 2 ) -c(x 1 ∨ x 2 , y 1 ∨ y 2 )] .
Note that this theorem is in spirit close to the one derived in Dutang et al. (2014, Theorem 1) but here we consider the case η = 1 and we have a process W much more complicated. The covariance structure of W can be obtained after tedious computations involving mainly the use of the covariance structure of W and the homogeneity properties of c x and c y . It is given by

E W t 1 c 1-ω W t 2 c 1-ω = t 1 + t 2 -t 1 ∨ t 2 + t 1 ∧ t 2 c 1-ω c 2 x 1, 1 -ω ω + c 2 y 1, 1 -ω ω 1 -ω ω -c x 1, 1 -ω ω t 1 ∧ t 2 + c(t 1 ∧ t 2 , t 1 ∨ t 2 1-ω ω ) c 1-ω -c y 1, 1 -ω ω t 1 ∧ t 2 + c(t 1 ∨ t 2 , t 1 ∧ t 2 1-ω ω ) c 1-ω + 2 c 1-ω c x 1, 1 -ω ω c y 1, 1 -ω ω × c t 1 ∧ t 2 , t 1 ∨ t 2 1 -ω ω + c t 1 ∨ t 2 , t 1 ∧ t 2 1 -ω ω .
We are now able to derive the asymptotic behavior of our main statistic appearing in the estimating equation ( 9).

Theorem 2. Under the assumptions of Theorem 1, for any s < 0, we have that

1 m m j=1 Z 1-ω,n-j+1,n Z 1-ω,n-m,n s = 1 1 -s + s √ m 1 0 t -s 1 t W t c 1-ω -W 1 c 1-ω dt - λ s (1 -s)(1 -s + τ ) ξ 1-ω c τ 1-ω 1 √ m + o P 1 √ m .
From now on we will denote the true value of δ 1-ω by δ 0 . Note that δ 0 = δ 0 ( Z 1-ω,n-m,n ). We want to prove the convergence in distribution of L m (x, y). This requires to establish the weak convergence of K 1-ω,m which is a function of the MDPDE δ 1-ω and Z 1-ω,n-m,n = Q n (1). Consequently, in the next theorem, we derive the convergence in distribution of the vector

  Q n (1) δ 1-ω  
correctly normalized, in case where τ is replaced by the true value τ 0 , by a mis-specified value, say τ , possibly different from τ 0 , or by a suitable external estimator τ .

Theorem 3. Under the assumptions of Theorem 1 we have that

√ m   k n Q n (1) -c 1-ω δ 1-ω -δ 0   d -→   A Q A τ  
where

A Q := c 1-ω W 1 c 1-ω + λc 1-ω c -τ 0 1-ω ξ 1-ω -1 τ 0 and • in case τ = τ 0 A τ 0 := (1 + 2α)(1 + 2α + τ 0 )(1 + 2α + 2τ 0 ) τ 2 0 (1 + τ 0 + 2ατ 0 + 4α 2 ) 2α 1 0 t 2α 1 t W t c 1-ω -W 1 c 1-ω dt -(2α + τ 0 )(1 + τ 0 ) 1 0 t 2α+τ 0 1 t W t c 1-ω -W 1 c 1-ω dt
• in case τ = τ , a mis-specified value

A τ := (1 + 2α)(1 + 2α + τ )(1 + 2α + 2 τ ) τ 2 (1 + τ + 2α τ + 4α 2 ) 2α 1 0 t 2α 1 t W t c 1-ω -W 1 c 1-ω dt -(2α + τ )(1 + τ ) 1 0 t 2α+ τ 1 t W t c 1-ω -W 1 c 1-ω dt +λ ξ 1-ω c τ 0 1-ω -2α (1 + 2α)(1 + 2α + τ 0 ) + (2α + τ )(1 + τ ) (1 + 2α + τ )(1 + 2α + τ 0 + τ ) - τ 2 (1 + τ + 2α τ + 4α 2 ) τ 0 (1 + 2α)(1 + 2α + τ )(1 + 2α + 2 τ ) • in case τ = τ , an external estimator satisfying √ k q 1 k n ( τ -τ ) = O P (1) with k an interme- diate sequence such that m k /q 1 (k/n) -→ 0 A τ = A τ 0 .
As is clear from Theorem 3, our estimator δ 1-ω is asymptotically unbiased in the sense that the expectation of A τ 0 is zero, whatever the value of λ. This is also the case for

A τ but not for A τ if τ = τ 0 .
As a corollary of Theorem 3, we can deduce the convergence in distribution of L m (x, y).

Corollary 4. Under the assumptions of Theorem 3 we have that

√ m L m (x, y) -L(x, y) d -→ -x A Q + xK 1-ω A τ + xK 1-ω λ c -τ 0 1-ω ξ 1-ω -1 τ 0 .
Again we can easily check that in case τ = τ 0 or if τ is replaced by an external estimator satisfying

√ k q 1 k n ( τ -τ ) = O P (1)
with k an intermediate sequence such that m k /q 1 (k/n) -→ 0, our estimator of L(x, y) is asymptotically unbiased. When τ is mis-specified, one possibly loses the asymptotic unbiasedness. However, as illustrated in Section 4, in that case the estimators still perform quite well with respect to bias and they outperform estimators which are not corrected for bias. Note that our result in Corollary 4 is a pointwise weak convergence result. This could be extended to finite dimensional convergence, though this does not yield interesting additional insights. Obtaining weak convergence results for the stochastic process

√ m L m (x, y) -L(x, y)
is in our context very complicated as we do not have an explicit estimator for δ 0 . This will though be the subject of future research.

Simulation study

In this section, we illustrate the finite sample properties of our robust and bias-corrected estimator L m (x, y) through a simulation study. Using the homogeneity property, we consider only the estimation of L(t, 1 -t) for t ∈ (0, 1), corresponding to the Pickands dependence function, and the following distributions:

• the symmetric logistic model, for which L(x, y) = (x 1/s + y 1/s ) s . We set s = 1/2;

• the bivariate Pareto of type II model, called BPII(p), for which L(x, y) = x+y-(x -p +y -p ) -1/p .

We set p = 3 and 4;

• the Cauchy distribution, for which L(x, y) = (x 2 + y 2 ) 1/2 .

These distributions have already been considered in

τ k (x, y) := 1 log r log ∆ k,a (rx, ry) ∆ k,a (x, y) -1 ∨ 0 ( 12
)
where

∆ k,a (x, y) := a -1 Ľk (ax, ay) -Ľk (x, y) with Ľk (x, y) := 1 k k j=1 a -1 j K(a j ) L k (a j x, a j y)
and a j := j k+1 and K a kernel such that Concerning the contamination, we use the following algorithm, applied to each of the above mentioned distributions:

1 0 K(u)u -1 2 du < ∞. According to
• We simulate n pairs (X i , Y i ), i = 1, ..., n, independently with unit Fréchet margins;

• We draw two lines both with origin the point (X 990,n , Y 990,n ) and to the points (X n,n , Y j ) and (X k , Y n,n ), respectively, where j and k correspond to the indices associated to the maximum of the other component (note that the point (X 990,n , Y 990,n ) is still in the main cloud, see Figure 2 for an illustration). We put n 0 /2 := (nε)/2 variables on each of these lines, according to a Fréchet distribution. Our contaminated sample has a size n * = n + n 0 ;

• We transform the margins of the contaminated sample empirically (with the ranks) into unit Pareto;

• We apply our MDPD criterion on the relative excesses. 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 q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q qq q q q qq 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 qq q q q q q q q q q q q q qq 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 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 q q q q q q q q q q q q q q q q q qq 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 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 q q q q q q q q q q q q q q q q q q q q q q qq q q qq 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 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 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 qq 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 qq q q q q q q q q q q q qq 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 qq 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 qq 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 q q q q q q q q q q q q q q q qq 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 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 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 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 q q q q q q q q qq 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 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 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 qq q q q q q q q q q q q q q q q q q qq q q q q q q qq 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 qq q q qq qq 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 0 200 400 600 800 1000 Thus this sample is the one on which our MDPD criterion will be applied.

The percentage of contamination is set to ε = 0% and 5%, while n = 1000 and the procedure is repeated 500 times. Two positions have been considered: (x, y) = (0.5, 0.5) and (x, y) = (0.4, 0.6).

Note that other percentages of contamination have been used, and, as expected, increasing the percentage of contamination deteriorates all the estimators but it increases the superiority of our estimator compared to the other ones in particular in terms of MSE. Again, to keep the length of the paper reasonable, they are omitted. All the simulations have been done with Fortran and the NAG library. The code is available upon request from the authors.

First, a comparison between Figures 3-4 and Figures 5-6 allows us to see the impact of estimat-

ing the second order parameter τ rather than fixing it at the value one. Indeed, in Figures 34, corresponding to 0 and 5% of contamination respectively, we show the behavior of our robust and bias-corrected estimator L m (0.5, 0.5) after replacing τ by the external estimator defined in [START_REF] Feuerverger | Estimating a tail exponent by modelling departure from a Pareto distribution[END_REF], whereas in Figures 56we fixed τ at 1. Different values of α have been used: 0.25 (full orange line), 0.5 (dashed red line) and 1 (dotted red line). Note that for the Logistic and BPII(3) distributions τ 0 = 1, whereas for the Cauchy τ 0 = 2 and for the BPII(4) τ 0 = 1/2. On the left of each figure we show the mean over the 500 replications and on the right the MSE as a function of m. The horizontal reference line on the left corresponds to the true value of L(x, y). The difference between our estimators as a function of α is rather small for α less than 1 in both cases, τ estimated or fixed. Increasing α further (i.e. α > 1) leads to results which are again acceptable but worse than those shown in Figures 3 till 6. For this reason, we recommend a value of α less than one allowing a good efficiency of our estimates. This was also suggested in the original Basu et al. (1998) paper and more recently in papers devoted to this topic. A careful study of these figures indicates however that, in case τ estimated or fixed, α = 0.25 seems to be close to the best value, especially with respect to the MSE, although the differences between the values of α are rather small. Concerning the influence of estimating τ , except for the BPII(4) distribution, it is clear from all the other figures that using a fixed value for τ is the best approach, even if this fixed value is not the true one. Indeed, in general in that context, the bias is reduced and also the MSEs are smaller and more stable for a long range of values of m. The fact that this is not the case for the BPII( 4) is probably due to the proximity between the true value τ 0 and the lower bound (0.3) imposed to avoid the instabilities. At first sight, the worse results obtained by estimating τ compared to mis-specifying it, can appear a bit surprising, but this has already been observed in the literature where also τ is often set at the value 1 (see e.g. [START_REF] Feuerverger | Estimating a tail exponent by modelling departure from a Pareto distribution[END_REF] and is due, in our context, to the fact that our estimator τ k (x, y) defined in ( 12) is not robust (see also Figure 8 In order to evaluate the performance of our robust and bias-corrected estimator L m (x, y), we compare it with three other estimators already introduced in the literature:

• the empirical estimator L m (x, y);

• the bias-corrected estimator L m,0. but we show here only the best one. The results are displayed in Figure 7 for 0% of contamination and in Figure 8 for 5% of contamination. Again (x, y) = (0.5, 0.5) in these two figures.

We can observe that in the presence of contamination, our estimator L m (x, y) is always the best one in terms of bias and MSE, whatever the distribution. This is obviously expected since it is the only estimator which is robust, the other ones being very poor in that context. In case of no contamination, as illustrated in Figure 7, our robust and bias-corrected estimator L m (x, y) still performs quite well, always better than the empirical one L m (x, y) and often similarly to L m,990 (x, y). Also note that in that context, L m,0.4,990 (x, y) has nice performance in terms of bias, but is worse than the others in terms of MSE. Moreover, unlike the latter estimator, L m,990 (x, y) and our bias-corrected and robust estimator L m (x, y) have very smooth sample paths.

Finally, in Figure 9, we consider another position, (0.4, 0.6), again for the Logistic and BPII [START_REF] Beirlant | Bias-reduced estimators for bivariate tail modelling[END_REF] distributions with τ = 1 and 0% and 5% of contamination. In that context, the estimation is slightly more difficult than in case x = y with a reduced stable part of the mean and the MSE, but again the influence of α is rather small, the value of α = 0.25 being again the best one.

Appendix: Proofs

Proof of Theorem 1

Let U i := 1 -F X (X i ) and V i := 1 -F Y (Y i ), i = 1, ..., n, with the convention that U 0,n = V 0,n = 0 and U n+1,n = V n+1,n = 1. Now define            S 1 (x, y) := n i=1 1l {U i ≤x and V i ≤y} S 2 (x, y) := n i=1 1l {U i ≤x or V i ≤y} .
We have

1 m S 1 k n x, k n y 1 -ω ω = P(U i ≤ k n x, V i ≤ k n y 1-ω ω ) q(k/n) + k m S 2 ( k n x, 0) k -x + k m S 2 (0, k n y 1-ω ω ) k -y 1 -ω ω - k m S 2 k n x, k n y 1-ω ω k - n k P U i ≤ k n x or V i ≤ k n y 1 -ω ω .
According to convergences (6.3) and (6.4) in Draisma et al. (2004) and using the fact that

m k → > 0 we have 1 m S 1 k n x, k n y 1 -ω ω = P(U i ≤ k n x, V i ≤ k n y 1-ω ω ) q(k/n) + 1 √ m 1 √ W (x, 0) + W 0, y 1 -ω ω -W x, y 1 -ω ω + o 1 √ m
a.s. uniformly on [0, b] 2 , for some positive real b. Thus, using the second order condition (SO)

and the fact that

√ mq 1 (k/n) → λ, we have a.s. uniformly on [0, b] 2 that 1 m S 1 k n x, k n y 1 -ω ω = c x, y 1 -ω ω + q 1 k n c 1 x, y 1 -ω ω + 1 √ m 1 √ W (x, 0) + W 0, y 1 -ω ω -W x, y 1 -ω ω + o 1 √ m . (13)
Define now

F n (r) := 1 n n i=1 1l { k n+1 Z 1-ω,i >r} = 1 n n i=1 1l 1- R X i n+1 < k n+1 1 r and 1- R Y i n+1 < k n+1 1 r 1-ω ω = 1 n n i=1 1l U i <U k r ,n and V i <V 1-ω ω k r ,n = 1 n S 1 U k r ,n-, V 1-ω ω k r ,n- (14)
where S 1 (x -, y -) denotes the left-hand side limit of S 1 at (x, y). We want to study the a.s.

convergence to zero of

sup b * ≤r<∞ √ m F n (r) q(k/n) -r -1 c 1-ω -λc 1 1 r , 1 r
1 -ω ω -W 1 r , 1 r
1 -ω ω , (15)
where b * is a positive real. To this aim, we use the following decomposition

√ m F n (r) q(k/n) -r -1 c 1-ω -λc 1 1 r , 1 r
1 -ω ω -W 1 r , 1 r
1 -ω ω ≤ √ m F n (r) q(k/n) -c n k U k r ,n , n k V 1-ω ω k r ,n -λc 1 n k U k r ,n , n k V 1-ω ω k r ,n - 1 √ W n k U k r ,n , 0 + W 0, n k V 1-ω ω k r ,n -W n k U k r ,n , n k V 1-ω ω k r ,n + √ m c n k U k r ,n , n k V 1-ω ω k r ,n -c 1 r , 1 r
1 -ω ω + √ c x 1 r , 1 r
1 -ω ω W 1 r , 0 + c y 1 r , 1 r
1 -ω ω W 0, 1 r
1 -ω ω +λ c 1 n k U k r ,n , n k V 1-ω ω k r ,n -c 1 1 r , 1 r
1 -ω ω + 1 √ W n k U k r ,n , 0 + W 0, n k V 1-ω ω k r ,n -W n k U k r ,n , n k V 1-ω ω k r ,n -W 1 r , 0 -W 0, 1 r
1 -ω ω + W 1 r , 1 r
1 -ω ω =: T 1 + T 2 + T 3 + T 4 .
A key argument to study these terms consists to use the following stochastic process property concerning uniform order statistics from Drees and Huang (1998): for any fixed value T > 0 one

has sup 0≤x≤T √ k n k U [kx],n -x + W (x, 0) → 0 a.s., (16)
sup 0≤y≤T √ k n k V [ky],n -y + W (0, y) → 0 a.s.. (17)
Combining ( 13) and ( 14) with ( 16) and ( 17), we deduce that T 1 → 0 a.s. uniformly for r ∈ [b * , ∞). Now, concerning T 2 , by applying the mean value theorem with x * , resp. y * , a random value between 1 r and n k U k r ,n , resp.

1 r 1-ω ω and n k V 1-ω ω k r ,n , we have T 2 ≤ √ m c x (x * , y * ) n k U k r ,n - 1 r + √ c x 1 r , 1 r 1 -ω ω W 1 r , 0 + √ m c y (x * , y * ) n k V 1-ω ω k r ,n - 1 r 1 -ω ω + √ c y 1 r , 1 r
1 -ω ω W 0, 1 r
1 -ω ω ≤ m k - √ c x (x * , y * ) √ k n k U k r ,n - 1 r + √ c x (x * , y * ) -c x 1 r , 1 r
1 -ω ω √ k n k U k r ,n - 1 r + √ c x 1 r , 1 r
1 -ω ω √ k n k U k r ,n - 1 r + W 1 r , 0 + m k - √ c y (x * , y * ) √ k n k V 1-ω ω k r ,n - 1 -ω ω 1 r + √ c y (x * , y * ) -c y 1 r , 1 r
1 - ω √ k n k V 1-ω ω k r ,n - 1 r 1 -ω ω + √ c y 1 r , 1 r
1 -ω ω √ k n k V 1-ω ω k r ,n - 1 r 1 -ω ω + W 0, 1 r 1 -ω ω .
Using the continuity of c x and c y , the convergence m k → together with ( 16) and ( 17), we deduce that T 2 → 0 a.s. uniformly for r ∈ [b * , ∞). Finally, for T 3 and T 4 we use the continuity of c 1 and W together with [START_REF] Gomes | Asymptotically unbiased" estimators of the tail index based on external estimation of the second order parameter[END_REF] and [START_REF] De Haan | Extreme Value Theory -An Introduction[END_REF] to deduce that T 3 → 0 a.s. and T 4 → 0 a.s. uniformly on r ∈ [b * , ∞). Altogether this implies the a.s. convergence to 0 of [START_REF] Gomes | Tail index estimation for heavy-tailed models: accommodation of bias in weighted log-excesses[END_REF]. Thus we have

√ m F n (r) q(k/n) -r -1 c 1-ω b * ≤r<∞ -→ W 1 r , 1 r
1 -ω ω + λc 1 1 r , 1 r
1 -ω ω b * ≤r<∞ weakly in D[b * , ∞), from which we deduce that √ m F n (r -1 ) q(k/n) -rc 1-ω 0≤r<b -→ W r, r 1 -ω ω + λc 1 r, r 1 -ω ω 0≤r<b
weakly in D[0, b ), for some positive real b . Finally, using the explicit expression of c 1 given in [START_REF] Dutang | Robust and bias-corrected estimation of the coefficient of tail dependence[END_REF] and applying Vervaat's lemma [START_REF] Vervaat | Functional limit theorems for processes with positive drift and their inverses[END_REF], we derive that

√ m F ← n (tq(k/n)) -1 - t c 1-ω 0<t≤c 1-ω b -→ - 1 c 1-ω   W t c 1-ω , t c 1-ω 1 -ω ω + λt t c 1-ω τ ξ 1-ω -1 τ   0<t≤c 1-ω b =: - 1 c 1-ω W t c 1-ω + λt ( t c 1-ω ) τ ξ 1-ω -1 τ 0<t≤c 1-ω b (18
F ← n (tq(k/n)) = k n Q n (t)+ O(1/m) uniformly, we deduce that √ m k n Q n (t) - c 1-ω t - c 1-ω t 2 W t c 1-ω -λ c 1-ω t t c 1-ω τ ξ 1-ω -1 τ = o P (1).
Now, following the lines of proof of Lemma 6.2 in [START_REF] Draisma | Bivariate tail estimation: dependence in asymptotic independence[END_REF], with some adjustments due to our different model, we can deduce that for all t 0 , ε > 0: sup

0<t≤t 0 t 3 2 +ε √ m k n Q n (t) - c 1-ω t - c 1-ω t 2 W t c 1-ω -λ c 1-ω t t c 1-ω τ ξ 1-ω -1 τ = o P (1),
which achieves the proof of Theorem 1.

Proof of Theorem 2. Recall that 1/m m j=1 ( Z 1-ω,n-j+1,n / Z 1-ω,n-m,n ) s can be rewritten as

1 0 (Q n (t)/Q n ( 1
)) s dt. Thus the main argument to prove Theorem 2 is to use the representation

Q n (t) = n k c 1-ω t   1 + 1 t √ m W t c 1-ω + λ √ m t c 1-ω τ ξ 1-ω -1 τ + o P t -( 1 2 +ε) √ m   (19)
which follows from Theorem 1. To this aim, we use a Taylor expansion of the integrand and thus we need that all the terms appearing in the square parenthesis in [START_REF] Peng | A practical way for estimating tail dependence functions[END_REF] tend to 0 uniformly.

To reach this goal we have to split the integrals into two parts, from 0 to m -κ and from m -κ to 1, for a suitable κ ∈ (0, 1/(1 + 2ε)).

Assuming κ > 1/2, we clearly have

m -κ 0 Qn(t) Qn (1)
s dt = o P (1/ √ m). Also, for ε < 0.5

1 m -κ Q n (t) Q n (1) s dt = 1 m -κ t -s dt + s √ m 1 m -κ t -s 1 t W t c 1-ω -W 1 c 1-ω dt + λs √ m ξ 1-ω c τ 1-ω 1 m -κ t -s t τ -1 τ dt + o P (1) √ m 1 m -κ t -s-1 2 -ε dt + o P 1 √ m = 1 0 t -s dt + s √ m 1 0 t -s 1 t W t c 1-ω -W 1 c 1-ω dt + λs √ m ξ 1-ω c τ 1-ω 1 0 t -s t τ -1 τ dt + o P (1) √ m 1 0 t -s-1 2 -ε dt + o P 1 √ m
from which Theorem 2 follows.

Proof of Theorem 3. The expression of A Q follows from Theorem 1. Below we give the details about the proof which leads to A τ 0 . Concerning A τ and A τ , the proofs are similar. Thus we will limit ourselves to the main differences.

Recall that ∂ ∆α(δ 1-ω ,τ 0 )

∂δ 1-ω δ 1-ω = δ 1-ω
= 0 by definition. Thus a Taylor series expansion around δ 0 combined with the boundedness of the third derivative of ∆ α (δ 1-ω , τ 0 ) with respect to δ 1-ω

leads to 0 = ∂ ∆ α (δ 1-ω , τ 0 ) ∂δ 1-ω δ 1-ω =δ 0 + ∂ 2 ∆ α (δ 1-ω , τ 0 ) ∂δ 2 1-ω δ 1-ω =δ 0 ( δ 1-ω -δ 0 )(1 + o P (1))
from which we deduce that

√ m( δ 1-ω -δ 0 ) = -   ∂ 2 ∆ α (δ 1-ω , τ 0 ) ∂δ 2 1-ω δ 1-ω =δ 0   -1 √ m ∂ ∆ α (δ 1-ω , τ 0 ) ∂δ 1-ω δ 1-ω =δ 0 (1 + o P (1)). ( 20
)
Since we only need the dominant term of the second derivative, it is sufficient to compute it at δ 1-ω = 0. Doing so, we obtain

∂ 2 ∆ α (δ 1-ω , τ 0 ) ∂δ 2 1-ω δ 1-ω =0 = 2(1 + α)(3 + 2α) ∞ 1 y -2(1+α) (1 -y -τ 0 ) 2 dy -4(1 + α) 2 ∞ 1 y -2(1+α) (1 -y -τ 0 )(1 -(1 -τ 0 )y -τ 0 )dy +α(1 + α) ∞ 1 y -2(1+α) (1 -(1 -τ 0 )y -τ 0 ) 2 dy -2(α + 1)(2α + 1) 1 m m j=1 E -2α j (1 -E -τ 0 j ) 2 +4α(1 + α) 1 m m i=1 E -2α j (1 -E -τ 0 j )(1 -(1 -τ 0 )E -τ 0 j ) -(α -1)(1 + α) 1 m m j=1 E -2α j (1 -(1 -τ 0 )E -τ 0 j ) 2 = τ 2 0 (1 + α)(1 + τ 0 + 2ατ 0 + 4α 2 ) (1 + 2α)(1 + 2α + τ 0 )(1 + 2α + 2τ 0 ) (1 + o P (1)) (21)
by applying Theorem 2. Now, concerning the first derivative

∂ ∆ α (δ 1-ω , τ 0 ) ∂δ 1-ω δ 1-ω =δ 0 = -2(1 + α) ∞ 1 y -2(1+α) (1 -y -τ 0 )(1 + δ 0 (1 -y -τ 0 )) -3-2α (1 + δ 0 (1 -(1 -τ 0 )y -τ 0 )) 1+α dy +(1 + α) ∞ 1 y -2(1+α) (1 -(1 -τ 0 )y -τ 0 )(1 + δ 0 (1 -y -τ 0 )) -2(1+α) (1 + δ 0 (1 -(1 -τ 0 )y -τ 0 )) α dy +2(1 + α) 1 m m j=1 E -2α j (1 -E -τ 0 j )(1 + δ 0 (1 -E -τ 0 j )) -1-2α (1 + δ 0 (1 -(1 -τ 0 )E -τ 0 j )) α -(1 + α) 1 m m j=1 E -2α j (1 -(1 -τ 0 )E -τ 0 j )(1 + δ 0 (1 -E -τ 0 j )) -2α (1 + δ 0 (1 -(1 -τ 0 )E -τ 0 j )) α-1 = - 2α(1 + α) √ m 1 0 t 2α 1 t W t c 1-ω -W 1 c 1-ω dt + (1 + α)(2α + τ 0 )(1 + τ 0 ) √ m 1 0 t 2α+τ 0 1 t W t c 1-ω -W 1 c 1-ω dt -q 1 k n ξ 1-ω c τ 0 1-ω τ 0 (1 + α)(1 + τ 0 + 2ατ 0 + 4α 2 ) (1 + 2α)(1 + 2α + τ 0 )(1 + 2α + 2τ 0 ) +δ 0 τ 2 0 (1 + α)(1 + τ 0 + 2ατ 0 + 4α 2 ) (1 + 2α)(1 + 2α + τ 0 )(1 + 2α + 2τ 0 ) + O P (δ 2 0 ) + o P 1 √ m
by again an application of Theorem 2 and using expansions of the integrands. Now recall that

δ 0 = δ 0 ( Z 1-ω,n-m,n ) = δ 1 Z 1-ω,n-m,n , 1 -ω ω 1 Z 1-ω,n-m,n by (6) = ξ 1, 1 -ω ω δ 1 Z 1-ω,n-m,n , 1
Z 1-ω,n-m,n (1 + o P (1)) by (4) = ξ 1-ω τ 0 q 1 1 Z 1-ω,n-m,n (1 + o P (1)) by (11) = ξ 1-ω τ 0 1 c τ 0 1-ω q 1 k n (1 + o P (1)),
using the regularly varying property of q 1 and Theorem 1. Consequently

∂ ∆ α (δ 1-ω , τ 0 ) ∂δ 1-ω δ 1-ω =δ 0 = - 2α(1 + α) √ m 1 0 t 2α 1 t W t c 1-ω -W 1 c 1-ω dt + (1 + α)(2α + τ 0 )(1 + τ 0 ) √ m 1 0 t 2α+τ 0 1 t W t c 1-ω -W 1 c 1-ω dt +o P 1 √ m . (22)
Combining ( 19) and ( 20) with ( 21) and ( 22) yields the expression of A τ 0 . Now, concerning A τ , the limit in probability of the second derivative is the same as in ( 21) with τ 0 replaced by τ . For the first derivative, similar arguments as for ( 22) lead to

∂ ∆ α (δ 1-ω , τ ) ∂δ 1-ω δ 1-ω =δ 0 = - 2α(1 + α) √ m 1 0 t 2α 1 t W t c 1-ω -W 1 c 1-ω dt + (1 + α)(2α + τ )(1 + τ ) √ m 1 0 t 2α+ τ 1 t W t c 1-ω -W 1 c 1-ω dt +q 1 k n ξ 1-ω c τ 0 1-ω 2α(1 + α) (1 + 2α)(1 + 2α + τ 0 ) - (1 + α)(2α + τ )(1 + τ ) (1 + 2α + τ )(1 + 2α + τ 0 + τ ) + τ 2 (1 + α)(1 + τ + 2α τ + 4α 2 ) τ 0 (1 + 2α)(1 + 2α + τ )(1 + 2α + 2 τ ) +o P 1 √ m
which leads to the expression of A τ .

Finally, concerning A τ , the limit in probability of the second derivative is again the same as in [START_REF] Schmidt | Non-parametric estimation of tail dependence[END_REF] since τ is a consistent estimator of τ . For the first derivative, similar arguments as for [START_REF] Vervaat | Functional limit theorems for processes with positive drift and their inverses[END_REF], combined with the mean value theorem, lead to

∂ ∆ α (δ 1-ω , τ ) ∂δ 1-ω δ 1-ω =δ 0 = - 2α(1 + α) √ m 1 0 t 2α 1 t W t c 1-ω -W 1 c 1-ω dt + (1 + α)(2α + τ 0 )(1 + τ 0 ) √ m 1 0 t 2α+τ 0 1 t W t c 1-ω -W 1 c 1-ω dt +o P 1 √ m + O P 1 √ kq 1 ( k n )
.

Using our assumption on the intermediate sequence k, the expression of A τ follows.

Proof of Corollary 4. Note that

√ m L m (x, y) -L(x, y) = -x √ m K 1-ω,m -K 1-ω . (23)
Consequently, we have to look at

√ m K 1-ω,m -K 1-ω = √ m m k k n Q n (1) 1 + δ 1-ω -K 1-ω = √ m 1 + δ 0 1 + δ 1-ω m k k n Q n (1) -c 1-ω 1 + δ 0 - k m 1 + δ 1-ω -δ 0 1 + δ 0 K 1-ω + c 1-ω 1 + δ 0 = √ m k n Q n (1) -c 1-ω (1 + o P (1)) - √ m δ 1-ω -δ 0 K 1-ω (1 + o P (1)) - √ m K 1-ω - m k c 1-ω 1 + δ 0 (1 + o P (1)). ( 24
)
We have the joint distribution of the two first terms, so we only have to look at the last one.

Recall that, by definition (see ( 6))

K 1-ω = 1 -ω ω d 2 g 1, 1 -ω ω .
Now using [START_REF] Dutang | Robust and bias-corrected estimation of the coefficient of tail dependence[END_REF], we have

c 1 x , 1 -ω ω 1 x = x -1 1 -ω ω d 2 g 1 x , 1-ω ω 1 x g(1, 1) = x -1 1 -ω ω d 2 g 1, 1-ω ω g(1, 1)
by homogeneity of order 0 of g. This implies that

Combining ( 25) and (26), the quantity of interest can be rewritten as

√ m K 1-ω - m k c 1-ω 1 + δ 0 = g(1, 1)c 1-ω √ m 1 - 1 + δ k n , k n 1 + δ 0 = g(1, 1)c 1-ω √ m δ 0 -δ k n , k n + o P (1) = K 1-ω √ mq 1 k n c -τ 0 1-ω ξ 1-ω -1 τ 0 + o P (1) (27)
by [START_REF] Dutang | Robust and bias-corrected estimation of the coefficient of tail dependence[END_REF] and using the facts that δ 0 = c -τ 0 1-ω ξ 1-ω 1 τ 0 q 1 k n (1 + o P (1)) and

√ m q 1 (k/n) → λ ∈ R.

  ) and F ← (t) := inf{z : F (z) ≥ t} where denotes either X or Y . This convergence (1) is relevant in multivariate extreme value theory: indeed if F X and F Y are in the max-domains of attraction of some extreme value distributions G X and G Y , respectively, and if (1) is satisfied, then F is in the max-domain of attraction of an extreme value distribution G with marginals G X and G Y and with copula function determined by L. Estimating this function L is thus of interest and is the subject of this paper. With this aim in mind, consider a sample of size n drawn from F and an intermediate sequence m = m n , i.e. m → ∞ as n → ∞ with m/n → 0. Let X m,n , resp. Y m,n , denote the m-th order statistic among n realisations of the margins X, resp. Y . The empirical estimator of L is then given by

Figure 1 :

 1 Figure 1: Logistic distribution, (x, y) = (0.5, 0.5). Comparison between the empirical estimator L m (x, y) (full black line), a bias-corrected estimator L m,0.4,990 (x, y) from Fougères et al. (2015) (dashed-dotted blue line), and the bias-corrected estimator L m,990 (x, y) proposed in Beirlant et al. (2016) (dotted green line) in case of 5% of contamination (see Section 4 for their definitions).

,

  Beirlant et al. (2011), Goegebeur and Guillou (2013), Dutang et al. (2014). On the contrary, if η = 1 then we are in the case of asymptotic dependence, and in that context estimating L is the subject of interest. Our starting point is thus model (3) with d 1 + d 2 = 1.

  This distribution function belongs to the class of distribution functions introduced in Beirlant et al. (2009) and called Condition (R) in Dierckx et al. (2013). As shown in Beirlant et al. (2009), the distribution function of the relative excesses Z 1-ω /u given that Z 1-ω > u can, for u large, be approximated by an extended Pareto distribution (EPD) function given by

  For the sequel it is instructive to keep in mind that (3) satisfies the following second order multivariate regular variation condition introduced in Draisma et al. (2004): Condition (SO): Let (X, Y ) be a random vector with joint distribution function F and continuous marginal distribution functions F X and F Y such that

as shown in Lemma 1 of

 1 [START_REF] Dutang | Robust and bias-corrected estimation of the coefficient of tail dependence[END_REF].

  [START_REF] Fougères | Bias correction in multivariate extremes[END_REF]. They have shown in particular that they satisfy their second order condition. Since the latter implies our model[START_REF] Beirlant | Bias-reduced estimators for bivariate tail modelling[END_REF], these distributions also satisfy our model. Note that other types of distributions have been studied but since the results are similar to those included we omit them, in order to keep the length of the paper under control.For each pair (X i , Y i ), i = 1, ..., n, independently from one of these distributions, we transform the margins into unit Pareto using the empirical distribution functions. Then, we minimize the empirical density power divergence ∆ α (δ 1-ω , τ ) after replacing τ by 1, which can be the true value or a mis-specified one depending on the distribution, or by the external estimator τ k introduced in Beirlant et al. (2016):

Proposition 1 in 1 kn

 11 Beirlant et al. (2016), this estimator τ k (x, y) satisfies the condition √ kq ( τ -τ ) = O P (1) required by our Theorem 3. Moreover as advocated in Beirlant et al. (2016), we use for a and r the value 0.4 and K is taken as the power kernel K(t) := 6 t 5 1l {t∈(0,1)} . Concerning k, an extensive simulation study involving different percentages of contamination, different types of distributions and sample sizes, indicates that k = n/2 performs quite well in all the contexts.Thus, these values for a, r, K and k are those used throughout the paper. However, to avoid instabilities in the computation of our robust estimator of L due to the fact that τ k (x, y) can be too close to 0, we set τ k (x, y) at 0.3 if τ k (x, y) ∈ [0, 0.3] and at its definition[START_REF] Feuerverger | Estimating a tail exponent by modelling departure from a Pareto distribution[END_REF] otherwise.Similar proposals where second order rate parameters are bounded away from zero have already been used, for instance in[START_REF] Beirlant | Tail index estimation and an exponential regression model[END_REF] in the univariate framework or Beirlant et al.(2016) in the multivariate setting. Finally, since any stable tail dependence function satisfies max(x, y) ≤ L(x, y) ≤ x + y, this constraint has been imposed on all the estimators.

Figure 2 :Figure 2 ,

 22 Figure 2: Procedure of contamination

  in Dierckx et al., 2013). Recall that the first paper dealing with the estimation of the second order parameter τ in the multivariate framework is due to Fougères et al. (2015), and even in the univariate context where this topic has been extensively studied, no robust estimator exists. In the sequel, we pursue the simulations with τ fixed at the value 1 and also we keep α at the value 0.25.

.

  4,990 (x, y) proposed by Fougères et al. (2015): L m,0.4,990 (x, y) := L m (x, y) ∆ 990,0.4 (0.4x, 0.4y) -L m (0.4x, 0.4y) ∆ 990,0.4 (x, y) ∆ 990,0.4 (0.4x, 0.4y) -0.4 ∆ 990,0.4 (x, y) where ∆ 990,0.4 (x, y) := 0.4 -1 L 990 (0.4x, 0.4y) -L 990 (x, y); • the bias-corrected estimator proposed by Beirlant et al. (2016): L m,990 (x, y) := Ľm (x, y) -m 990 τ 990 (x,y) α990 (x, y) 990 (x,y) a -1 j L 990 (a j x, a j y) Note that in fact Fougères et al. (2015) have proposed two bias-corrected estimators for L(x, y)

) 20 weakly

 20 in D(0, c 1-ω b ]. Combining the delta method with the fact that

K 1 -ω = g( 1 , 1

 111 order 1 of c. Now recall that k n = q ← m n , i.e. with η = 1 and homogeneity of g.

Finally combining ( 23 Figure 3 :

 233 Figure3: Comparison between our robust and bias-corrected estimator L m (0.5, 0.5) with α = 0.25 (full orange line), 0.5 (dashed red line) and 1 (dotted red line). Here τ is estimated by[START_REF] Feuerverger | Estimating a tail exponent by modelling departure from a Pareto distribution[END_REF] and 0% of contamination. First row: Logistic; Second row: BPII(3); Third row: Cauchy; Fourth row: BPII(4) distributions. On the left: mean and on the right: MSE based on 500 replications.The horizontal line on the left is the true value of the parameter.

Figure 4 :

 4 Figure 4: Comparison between our robust and bias-corrected estimator L m (0.5, 0.5) with α = 0.25 (full orange line), 0.5 (dashed red line) and 1 (dotted red line). Here τ is estimated by[START_REF] Feuerverger | Estimating a tail exponent by modelling departure from a Pareto distribution[END_REF] and 5% of contamination. First row: Logistic; Second row: BPII(3); Third row: Cauchy; Fourth row: BPII(4) distributions. On the left: mean and on the right: MSE based on 500 replications.The horizontal line on the left is the true value of the parameter.

Figure 5 :

 5 Figure 5: Comparison between our robust and bias-corrected estimator L m (0.5, 0.5) with α = 0.25 (full orange line), 0.5 (dashed red line) and 1 (dotted red line). Here τ = 1 and 0% of contamination. First row: Logistic; Second row: BPII(3); Third row: Cauchy; Fourth row: BPII(4) distributions. On the left: mean and on the right: MSE based on 500 replications. The horizontal line on the left is the true value of the parameter.

  1], g is a continuous homogeneous function of order 0 and δ a function of constant sign in the neighbourhood of zero, with |δ| being a bivariate regularly varying function, that is, there exists a function ξ such that lim

t↓0 |δ|(tx, ty) |δ|(t, t) = ξ(x, y),

Acknowledgement

The authors thank the associate editor and the reviewers for their helpful comments which led to improvements of our first version of the paper. This work was supported by a research grant (VKR023480) from VILLUM FONDEN and an international project for scientific cooperation (PICS-6416).


Chao Lyu

Yann Capdeville

Liang

Lyu Chao

Capdeville Yann

Zhao Liang

Lyu & Al

Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation

come

Efficiency of the spectral element method with very high polynomial degree to solve the elastic wave equation

INTRODUCTION

Solving the acoustic and elastic wave equations numerically is a critical step for many research based on seismic data. It is especially important for most seismic imaging methods, from the exploration scale to the global Earth scale. Even though the Finite Difference (FD) Efficient high degree spectral element method still dominates in the seismic exploration community, the spectral element method (SEM) [START_REF] Maday | Spectral element methods for the incompressible navierstokes equations: State-of-the-art surveys on computational mechanics[END_REF][START_REF] Seriani | Spectral element method for acoustic wave simulation in heterogeneous media[END_REF][START_REF] Komatitsch | The spectral element method: an effective tool to simulate the seismic response of 2D and 3D geological structures[END_REF][START_REF] Komatitsch | Introduction to the spectral-element method for 3-D seismic wave propagation[END_REF][START_REF] Chaljub | Spectral element analysis in seismology[END_REF] has been gaining more and more popularity, especially in the academic community. It is often the chosen method for global or regional Earth scale seismic imaging developments based on adjoint methods and full waveform inversion method [START_REF] Capdeville | Towards global earth tomography using the spectral element method: a technique based on source stacking: Geophysical Efficient high degree spectral element[END_REF][START_REF] Tromp | Seismic tomography, adjoint methods, time reversal and banana-doughnut kernels[END_REF][START_REF] Fichtner | Full seismic waveform Efficient high degree spectral element tomography for upper-mantle structure in the australasian region using adjoint methods[END_REF][START_REF] Zhu | Structure of the European upper mantle revealed by adjoint tomography[END_REF][START_REF] Monteiller | Three-dimensional full waveform inversion of short-period teleseismic wavefields based upon the SEM-DSM hybrid method[END_REF][START_REF] Wang | The deep roots of the western pyrenees revealed by full waveform inversion of teleseismic p waves[END_REF][START_REF] Beller | Lithospheric architecture of the south-western alps revealed by multiparameter teleseismic full-waveform inversion[END_REF][START_REF] Trinh | Efficient time-domain 3D elastic and viscoelastic full-waveform inversion using a spectral-element method on flexible cartesian-based mesh[END_REF].

SEM is a finite element type method and, as such, is based on the weak form of the wave equation. Two versions of the SEM exist, one based on Chebyshev polynomials and one based on Legendre polynomials. In the following, we only use the Legendre version, which is the basis for many available programs. Compared to classical low-degree finite elements, the Legendre SEM is based on a tensorised high-degree polynomial approximation per element combined with a precise numerical quadrature associated with the so-called Gauss-Lobatto-Legendre (GLL) points. It has a spectral convergence with the element polynomial degree, thus leading to a low spatial dispersion. Moreover, its tensorial formulation leads to a diagonal mass matrix, which is a strong advantage for explicit time schemes. It can naturally and accurately handle free surface and material discontinuities. This capability to accurately model interface waves such as surface waves is one of the main reasons for its popularity.

Despite its quality, the SEM has one critical drawback: to warrant the accuracy of the method, each material discontinuity interface has to be explicitly meshed. Moreover, because this method is based on a tensorial formulation, the mesh needs to be based on quadrilateral (in 2-D) or hexahedron (in 3-D) elements. The meshing difficulty can strongly limit the applicability of the method, especially in 3-D. Note that some examples of SEM based on triangles exist but their efficiency is lower and with a limited polynomial degree range [START_REF] Komatitsch | Wave propagation in 2-D elastic media using a spectral element method with triangles and quadrangles[END_REF][START_REF] Mercerat | Triangular spectral element simulation of two-dimensional elastic wave propagation using unstructured triangular grids[END_REF][START_REF] Afanasiev | Modular and flexible spectral-element waveform modelling in two and three dimensions[END_REF]. Efficient high degree spectral element One of the key parameters when using SEM is the polynomial degree N used in each tensorial direction for each element. In seismology, N = 4 is often chosen [START_REF] Komatitsch | The spectral element method: an effective tool to simulate the seismic response of 2D and 3D geological structures[END_REF][START_REF] Komatitsch | Introduction to the spectral-element method for 3-D seismic wave propagation[END_REF] and sometimes N = 8 [START_REF] Chaljub | Solving elastodynamics in a solid heterogeneous 3-Sphere: a spectral element approximation on geometrically non-conforming grids[END_REF][START_REF] Capdeville | Coupling the spectral element method with a modal solution for elastic wave propgation in global earth models[END_REF]. Although it depends on the desired accuracy, the propagated distance and on how the minimum wavelength is estimated, the commonly admitted number of grid points per minimum wavelength (G) to obtain sufficient accuracy in a constant velocity medium is approximately G = 5 ∼ 6 for N = 4 and G = 4 ∼ 5 for N = 8 [START_REF] Priolo | A numerical investigation of Chebyshev spectral element method for acoustic wave propagation[END_REF][START_REF] De Basabe | Grid dispersion and stability criteria of some common finite-element methods for acoustic and elastic wave equations[END_REF][START_REF] Seriani | Dispersion analysis of spectral element methods for elastic wave propagation[END_REF]. A low G is obviously an advantage as it lowers the required computer memory and maybe the computing cost of a modeling. In particular, a low G can be important in the full waveform inversion (FWI) context. Indeed, FWI schemes operate on the forward and adjoint wavefield. To do so, some methods rely on the partial or compressed storage [START_REF] Komatitsch | Anelastic sensitivity kernels with parsimonious storage for adjoint tomography and full waveform inversion[END_REF][START_REF] Boehm | Wavefield compression for adjoint methods in full-waveform inversion[END_REF] or even full storage to avoid backward propagation [START_REF] Fichtner | Full seismic waveform Efficient high degree spectral element tomography for upper-mantle structure in the australasian region using adjoint methods[END_REF]. In 3-D, dividing G by 2 implies a storage requirement divided by 8, which is significant. In the 1-D case, it has been shown, for N = 60, G can be as low as about 2 ∼ 2.5 [START_REF] Priolo | A numerical investigation of Chebyshev spectral element method for acoustic wave propagation[END_REF]. Nevertheless, the most widely used degree N is still 4 not 60, which is for two reasons:

• Taking advantage of a large N and a low G imposes being able to use very large elements. For example, an element for N = 60 needs to be about 15 times larger than an element for N = 4 (assuming a constant G(N ) for the sake of simplicity).

In most realistic situations, this is not possible because of the necessity to mesh mechanical discontinuities and the domain geometry;

• The GLL points are not evenly spaced. They cluster toward the edges of the elements and this clustering is stronger with large N as shown in Figure 1 stability condition is governed by the minimum value of the ratio between the size of the grid cells and the P-wave velocity expressed as:

dt ≤ C dx min V P , ( 1
)
where C is the courant constant and V P is P-wave velocity. In the following, we refer to dt CFL as the largest stable time step of the Newmark scheme for a given test case.

As a consequence, the possibility to use very high degree spectral elements has not been considered as a valid option and has not been studied yet unlike for finite differences method [START_REF] Liang | Determining finite difference weights for the acoustic wave equation by a new dispersion-relationship-preserving method[END_REF]. Nevertheless, in the forward modeling context, the recent Efficient high degree spectral element introduction of the non-periodic homogenization [START_REF] Capdeville | 2-D non-periodic homogenization to upscale elastic media for P-SV waves[END_REF][START_REF] Cupillard | Non-periodic homogenization of 3-D elastic media for the seismic wave equation[END_REF] has solved many of the meshing difficulties and opens the door to the use of very high degree N . In the inverse problem context, models are generally smooth and using very high degree elements has always been an option.

The non-periodic homogenization is an asymptotic method designed to compute effective media in the case of deterministic media with no scale separation such as geological media. For a given maximum source frequency and a complex "true" medium, this method computes a smoothly varying effective medium in which the computed waveform is the same as a reference solution computed in the true medium, up to a controllable error. This is true for all type of waves, including surface, refracted, backscattered, etc, waves. In the forward modeling context, homogenization can be seen as a preprocessing step upscaling original media to the wavelength scale.

Regarding SEM, homogenization solves some aspects of the meshing issue: indeed, after homogenization, there is no more mechanical discontinuity to honor, except for the free surface. This makes the option of using very large elements with a high degree N possible, opening the door to a lower G than for the classical degree N = 4 and therefore allowing a lower memory. Nevertheless, about the computing cost, the fact that GLL points cluster toward the edges of the elements is still true and will go against the advantages of a lower G. The objective of this work is, therefore, to determine if using very high degree N can be an advantage, for both memory and computing time aspects for SEM and if so, to find if there is an optimal degree N .

The paper is organized as follows: we first give information about the different SEM codes, the criteria for the different test configurations and the methodology we use to address the paper's objective. We then perform our tests in 2-D and 3-D homogeneous media.

Finally, we perform the same test in simple heterogeneous media before concluding our work. For each program, we limit ourselves to the standard order 2 Newmark scheme for Efficient high degree spectral element the explicit time marching.

SEM CONSIDERATION AND NUMERICAL EXPERIMENTS

SETUP

SEM complexity consideration

As shown in the Appendix A, for a given accuracy, if the maximum time step allowed by the CFL condition given in equation 1 can be reached, the complexity C(N ), that is the number of operations of any SEM program, as a function of the degree N , is scaling as

C(N ) ∝ G d+1 (N )N 2 , ( 2
)
where d is the dimension of the problem. If the maximum time steps cannot be reached (typically for long time series for which the time step error can be large), the complexity scales as

C(N ) ∝ G d (N )N . (3)
In both cases, it shows that the decreasing G as a function of N is competing with the increasing N . Knowing that the decrease of G as a function of N is non-linear and cannot go below 2, it is already clear that there is a limit after which increasing N will not be efficient. But before that, an optimal N exists and needs to be determined numerically, which is the purpose of the next sections.

SEM codes

For a given source maximum frequency f max , geometry, elastic properties, and signal duration t max , there are two aspects in assessing whether using the very high degree in SEM Efficient high degree spectral element can be a good idea. First, we consider the amount of computer memory required to perform the simulation. This is directly related to the number of points per wavelength G and does not depend much on the particular code implementation of SEM. Second, we consider the computing time needed to perform the simulation. This is strongly related to code implementation. To mitigate a possible bias in our results, we use three different programs:

• SPEC2DY. This program has originally been written by G. Festa [START_REF] Festa | The newmark scheme as velocity-stress time-staggering: an efficient implementation for spectral element simulations of elastodynamics[END_REF] and has been widely modified over the years, but the philosophy and the core of the original program are still the same. In particular, its efficiency has not been seriously optimized.

• SEM3D: This program is a 3-D version of the above 2-D program. It has been originally written by G. Festa and E. Delavaud [START_REF] Delavaud | Simulation numérique de la propagation d'ondes en milieu géologique complexe: application à l'évaluation de la réponse sismique du bassin de caracas (venezuela)[END_REF]. In contrast to its 2-D version, its efficiency has then been optimized by a team of the CEA (Commissariat  l'nergie atomique, France);

• SPECFEM (2D and 3D). These two programs are from the popular SPECFEM program suites, both the 2-D and 3-D versions are well maintained and optimized [START_REF] Komatitsch | The spectral element method: an effective tool to simulate the seismic response of 2D and 3D geological structures[END_REF][START_REF] Komatitsch | Introduction to the spectral-element method for 3-D seismic wave propagation[END_REF].

Each of these three codes uses an order 2 Newmark time marching scheme.

Numerical experiments setup and models

In order to measure the numerical efficiency of SEM with degree N , we set up a series of 2-D and 3-D numerical experiments in homogeneous and heterogeneous elastic models. To make the efficiency measurements meaningful, we need to perform our experiments for a fixed accuracy. There are many ways to measure accuracy and we choose one that is close to many realistic situations: a rectangular domain Ω with a free surface on the top, shallow sources and shallow receivers with recorded waveforms in the time domain. Efficient high degree spectral element

We perform our tests in three different models: 2-D homogeneous, 3-D homogeneous and 2-D heterogeneous. We did not test the 3-D heterogeneous case because of large computing resources necessary to obtain an accurate reference solution in that case, but we believe that this does not detract from the main conclusion of this work. For the 2-D and 3-D homogeneous cases, we use V P = 3.4 km/s and V S = 2.0 km/s for the P and S wave velocities and ρ = 2000 kg/m 3 for the density. For the heterogeneous case, a monochromatic oscillatory heterogeneity is added on top of the homogeneous model. The detailed heterogeneity is defined in the 2-D heterogeneous experiment section below.

In the homogeneous case, the maximum frequency and S-wave velocity makes it possible to define the minimum wavelength of the propagating wavefield in the far-field as

λ min = V S f max . (4)
In the following, we measure each spatial distances as a function of λ min and time as a function of

t min = 1 f max . ( 5
)
The source is a vertical point force located at a 3 λ min distance below the free surface and its time wavelet is a Ricker function (second derivative of a Gaussian) of central frequency f 0 = 10 Hz with maximum frequency f max ≃ 3f 0 . This estimation of the maximum frequency is important because it determines the minimum wavelength and therefore G(N ).

This estimate can be changed and it will change the obtained G(N ) and it implies the absolute value of G must be taken with caution. Nevertheless, this estimate does not affect the relative results between different degrees.

To measure the error, we rely on two bins of receivers at two different epicentral distances, a short one (20λ min ) and a long one (200λ min ). Each bin is wide enough to contain at least one element, even for the largest degree N that we tested. The main reason for doing Efficient high degree spectral element so is that error is not constant within an element and depends on where exactly the receiver is located [START_REF] Moczo | 3-D finite-difference, finite-element, discontinuous-Galerkin and spectral-element schemes analysed for their accuracy with respect to P-wave to S-wave speed ratio[END_REF]. To average out this effect, we use many receivers within at least one element.

The boundary conditions are free normal stress conditions all around the domain to fully exclude the influence of absorbing boundaries. The chosen domain is wide enough to ensure that no reflecting waves from boundaries affect the results. The time duration of the signal is long enough to ensure that the full waveforms include P, S and the Rayleigh wave phases.

2-D homogeneous experiments

For the 2-D tests, Ω is a 900×450 λ 2 min rectangular domain. Each of the different epicentral 20 λ min and 200 λ min receiver bin has a size of 15 × 15 λ 2 min and contains 2601 receivers. Ω

and an example of energy snapshot of the wavefield at t = 270 t min is displayed in Figure 2.

The recorded signal at each receiver last for 60 t min for the 20 λ min distance receiver bin and 300 t min for the 200 λ min distance receiver bin.

3-D homogeneous experiment

For the 3-D test, we use a homogeneous media with the same elastic properties as in the 2-D homogeneous case above. Here, because of computing resources limitation, we only use a relatively small model and a short distance receiver bins (20 λ min ). For this test, Ω is 300 × 300 × 150 λ 3 min parallelogram. The receiver bin is a 15 3 λ 3 min cube, just below the surface, 20 λ min away from the source (see Figure 3), containing 1331 receivers. The 200 λ min receiver bin case has not been studied in 3D because the computing resources required to do so are beyond our capacity (indeed, to avoid domain border reflections, for the 200 λ min case an even larger domain is required). Efficient high degree spectral element In order to further analyze the efficiency of SEM as a function of the degree N in the heterogeneous 2-D case, we use the same geometry as for the 2-D homogeneous test and only focus on the long epicentral distance bin (200 λ min ). To keep our analysis simple, the heterogeneous mechanical properties have been chosen with the following form:

                       λ(x) = λ 0 f (x) µ(x) = µ 0 f (x) f (x) = 1 + a cos( 2π λ h k x • x) + cos( 2π λ h k z • x) k x = (cos(t 1 π/180 • ), sin(t 1 π/180 • )) k z = (cos(t 2 π/180 • ), sin(t 2 π/180 • )) (6)
where λ and µ are the Lamé elastic coefficients for the heterogeneous models, λ 0 and µ 0 the constant Lamé coefficients corresponding to V P = 3.4 km/s, V S = 2.0 km/s, ρ = 2000 kg/m 3 . The density is kept constant. We use a = 0.05, t 1 = 45 • , and t 2 = 135 • , corresponding to two orthogonal directions shown in Figure 4. In the following, only three values for λ h are tested, λ h = 2λ min , λ min and λ min /2. Note that a real model, once homogenized, has a continuous spectrum of heterogeneity, not just a single wavelength as here. In general, geological media lead to amplitude spectrum of heterogeneities that decreases as 1/λ. λ h = λ min /2 in equation 6 correspond to a strongly heterogeneous case, Efficient high degree spectral element

SEM mesh

We base our experiments on a trivial regular mesh made of n x × n z square elements of size L e in each direction. We keep n x = 2n z and the domain size constant. The element size L e can only be tuned by changing the number of horizontal elements n x . In each element, the polynomial expansion of degree N is used in each direction. The mesh is fully characterized by the (L e , N ) quantities. For a given mesh (L e , N ), we can compute dx min , Efficient high degree spectral element the minimum distance between two GLL points:

dx min (L e , N ) = L e 2 dξ N min . (7)
The number of grid points per minimum wavelength G in one direction is defined as follows:

G(L e , N ) = (N + 1) λ min L e . ( 8
)

Error measurement

There are several ways to measure the error of a given modeling with respect to a reference solution. One could for example compute the L 2 norm of the difference between computed and reference wavefield all over the domain for the final time step. Here, our choice is more closely related to situations faced in geophysics: receivers are located on the free surface or at shallow depths (in boreholes), with relatively long time series and for different epicentral distances. For a set of N r receivers, located in {x r , r ∈ {1..N r }}, the error E is computed as

E 2 = r=1,Nr tmax 0 (u -u ref ) 2 (x r , t) dt r=1,Nr tmax 0 (u ref ) 2 (x r , t) dt , ( 9
)
where u is the computed displacement including all the components, u ref is the reference solution.

The numerical error for a given model and receiver bin is jointly affected by the mesh design (L e , N ) and the time step dt:

E = E(L e , N, dt) .
(10) Efficient high degree spectral element

Reference Waveforms

The reference solution u ref should ideally be an analytical solution. Unfortunately, such analytical solutions are only available for simple models. In many cases, we need to rely on "converged" numerical solutions computed with a very small time step and spatially over-sampled numerical solution. Here, for the sake of completeness, we first show the equivalence of analytical and "converged" numerical solutions for 2-D and 3-D homogeneous half spaces. Afterward, we assume this equivalence holds for all our tests and all the reference waveforms are "converged" numerical solutions.

We first perform a 2-D comparison between a converged SEM solution and an analytical solution for the short distance receiver bin. The analytical solution was obtained using the program "EX2DDIR" (the source code can be found in www.spice-thn.org), which is based on the Caniard-de Hoop technique [START_REF] De Hoop | A modification of cagniard's method for solving seismic pulse problems[END_REF][START_REF] Johnson | Green's function for Lamb's problem[END_REF]. For the SEM solution, a 180×90 λ 2 min domain with a structural 100×50 spectral elements with degree 12 is designed. The resulting G value is about twice that recommended by [START_REF] Priolo | A numerical investigation of Chebyshev spectral element method for acoustic wave propagation[END_REF]. We finally use a very small time step, dt ≈ 1 250 dt CFL . The agreement between the two solutions is displayed in Figure 5 for a representative receiver. The total misfit, computed according to equation 9, is E = 1.5 × 10 -6 .

We then perform a similar test, but in 3-D. For the analytical solution, we use the program "CANHFS" (personal communication) which calculates the Green function and is also based on the Caniard-de Hoop technique. For the SEM solution, a 180×180×90 λ 3 min domain with a structural 100×100×50 spectral elements mesh with degree 12 is designed.

The resulting G value is once again about twice that recommended by [START_REF] Priolo | A numerical investigation of Chebyshev spectral element method for acoustic wave propagation[END_REF]. We finally use a very small time step, dt ≈ 1 100 dt CFL . The agreement between the two solutions is displayed in Figure 6 for a representative receiver. The total misfit is here, is E = 2.2 × 10 -3 . From these two tests, we conclude that we can replace the analytical solutions by con-Efficient high degree spectral element The residual (×500) is displayed in dotted line. The receiver is a typical receiver from the short distance bin (20 λ min distance). Efficient high degree spectral element verged numerical solution in our analysis. Indeed, our misfit threshold is 1%, and the misfit between converged and analytical solutions is much smaller than 1%. For the rest of the paper, we assume that we can use converged numerical solutions as reference solutions.

Finally, note the fact that the agreement between numerical and analytical solutions is not as good in 3-D than in 2-D is more related to our difficulty to accurately use the analytical solution, rather than related to a larger error in the 3-D numerical modeling. The 3-D misfit is nevertheless still much smaller than 1% and we did not investigate the code problems any further.

G AS A FUNCTION OF THE DEGREE N

Here, we study the minimum number of GLL points per wavelength G required to reach the 1% error threshold. To do so, we start by choosing a time step dt small enough so that the error due to the time marching scheme can be ignored. The error is then dominated by the spatial error and only depends on L e and N (E = E(L e , N )). Knowing that a fixed size domain is used, the elements size L e is controlled by the number of elements n z in the z direction (the number of elements in the x direction is tied to n z by the relation n x = 2n z ).

For given degree N , receiver bin (20λ min or 200λ min ) and model, we proceed as follow to determine G so that the error E = 1%:

• we start with a large n z so that the error is below 1%;

• we then gradually decrease n z until we find its values so that the error is just below and just above 1%;

• we compute G for each n z available (discrete) values and we finally obtain the G to reach exactly 1% using a quadratic interpolation.

An example of the procedure for N = 40 is given in Figure 7. Efficient high degree spectral element Following this procedure, we compute the 1% error G for a set of degree N (4, 8, 12, 16, 20, 24, 30 and 40), in our four 2-D models (homogeneous, heterogeneous with λ h = 2λ min , λ min and 0.5λ min ) and the 3-D homogeneous model for the two receiver bins. The results are given in Figure 8, Figure 9 and Table 1.

For the homogeneous case, it can be noted that, for N = 40, G ≃ 2.5 is reached in 2-D as well as in 3-D and for both short and long-distance bins. For the long-distance bin, it corresponds to twice fewer points per wavelength for N = 40 than N = 4. This implies that factor of 4 in 2-D and 8 in 3-D of computer memory can be saved using high degrees compared to low degrees. As expected, G is lower for short distance than for long, but this difference is significant only for low degree. For the heterogeneous case, Figure 9 displays a comparison of the G obtained in the three heterogeneous models versus the one obtained in the homogeneous case. It is worth noting that the sampling rule of thumbs used in homogenization [START_REF] Capdeville | Residual homogenization for elastic wave propagation in complex media[END_REF] is

G ε 0 ≃ G homo (1 + 1 2ε 0 ) , (11)
where G homo is G in the homogeneous case, G ε 0 is G in the heterogeneous case and ε 0 = λ h /λ min . Here we have ε 0 = 2, 1 and 0.5 , which leads to a G ε 0 /G homo ratio of 1.

OPTIMAL TIME STEP AS A FUNCTION OF THE DEGREE N

In this work, we only use the standard second-order explicit Newmark scheme. The error is therefore a quadratic function of dt. The time step is always subject to the CFL condition equation 1: dt ≤ dt CFL .

For a given degree N , we proceed in a similar way to the previous section to find the maximum time step dt c such that the error remains just below 1%. We first choose a large G (twice the optimal G) so that the spatial error is much smaller than 1%. In that case, the error is only dependent upon dt. We then find dt c by the method of trial and error and a quadratic interpolation. We finally check that the obtained dt c leads to an error indeed just below 1%. An example is given in Figure 10 for N = 4. However, this procedure is limited by the CFL condition. Indeed, for high degrees, the error cannot reach 1% even for dt = dt CFL . In that case, dt c cannot be determined and is set to dt CFL . For low degrees, it appears that dt c is almost independent of the degree and that it only depends on the receiver Efficient high degree spectral element It can be seen from Figure 10 that, for N = 4, dt = dt CFL leads to an error of 4.5%, which is large. In such a case, a higher order time scheme would be necessary to take advantage of a large time step.

From this section, we can conclude that the error in time and space are almost independent.

GLOBAL COMPUTING TIME AS A FUNCTION OF DEGREE N

We finally evaluate the computing time as a function of the degree N for the optimal sampling G and time step dt determined in the two previous sections. If the gain in memory For each case, we choose as reference the computing time obtained using SPECFEM2D or 3D programs with N = 4, normalized to one. In this section, the gains or losses in computing time are always in comparison to this reference. Each computing time measurement is obtained using the optimal G and dt, performing five runs and averaging the obtained elapsed computing time. We used 40, 280 and 100 computing cores for the 2D homogeneous, 3D homogeneous model and 2D heterogeneous models respectively. Efficient high degree spectral element

2-D homogeneous model case results

Figures 12 and 13 show the results for the homogeneous 2-D case, for the short and long distance bins respectively. For the short distance receiver bin, it can be seen that the computing time decreases with the increasing degree for both programs until degree 16 or 20 and then increases again. To its maximum, the high degree computing time gain compared to the low degree one is about a factor of 2. For the long distance bin, a similar pattern can be observed. Nevertheless, the computing time does not increase much after a degree 16. The pattern is slightly different for each program, but the overall result is similar:

3-D homogeneous model case results

For the 3-D case, we performed only the short distance bin case, because of computing resources limitation. The computing time as a function of N obtained for the two programs is displayed in Figure 14. First, it can be noted that there is no fundamental difference in computing time between the two programs. Next, the trend is very similar to the 2-D homogeneous short distance receiver case: first a decrease of the computing time with the degree and then an increase after a minimum. The minimum computing time is obtained in the 12-16 degree range.

For the 3-D long range case, we can only speculate that the results would be similar to the 2-D results. The short distance bin results are similar, and there is no specific reason to expect a difference for the long distance. It is nevertheless unfortunate that we can not check this point.

2-D heterogenous case results

We finally perform a test in the heterogeneous models described earlier. This test is only performed in 2-D for the long distance bin using the SPEC2DY program. We did not test the 3-D model because of computing resource limitation. We also did not test SPECFEM2D in that test because there is no simple way to input oscillating models in this program.

In Figure 15 the computing times as a function of the degree N for three 2-D heterogeneous models are shown (for λ h = 0.5λ min , λ min and 2λ min ) for the long distance bin. In Figure 13, it can be noted that the cases λ h = λ min and λ h = 2λ min display little differences compared to the homogeneous long distance case. The remaining differences are mostly due to the fact that the heterogeneous models have a smaller dt CFL (because of a denser mesh) compared to the long receiver homogeneous case. The case λ h = 0.5λ min is similar, but the minimum computing time occurs for lower degrees, 16-20, compared to 20-24 for the previous cases. Efficient high degree spectral element

DISCUSSION

In this work, to assess if using very high degree (N > 8) in SEM is of any interest, we had to make several choices. This first one is the error threshold that we chose to 1%. Obviously, one could have made a different choice and this would have affected the optimal G, dt and measured computing time results. For example, for the short distance bin, using a threshold to 5% moves the optimal G for N = 40 from 2.4 to 1.97. Nevertheless, it does not change either the observed trend or the conclusions of this work. The way of estimating the source maximum frequency f max has the same effect on G, but does not change the conclusions either. The second one is about the chosen programs to perform the computing time measurements. This aspect is problematic as different programs and different hardware can potentially significantly affect the computing time results. For example, we did not try the GPU spectral element implementations for which the behavior can be different. This work

shows that very high degree SEM can be interesting, but one should keep in mind that, for a significantly new or different SEM implementation, or significantly different hardware, a benchmark should be run to determine which degree is the best from the point of view of computing time. Nevertheless, the memory gain results would remain unchanged.

Using very high degree SEM implies using very large elements, much larger than for low degrees. In a forward modeling context, the interest of very high degree SEM therefore strongly relies on the homogenization technique. It is the key to remove mechanical discontinuities: it allows to release the mesh constraints and thus allows the use of large elements. In the inversion context, the models are smooth anyway and it is usually simple to use very large elements. In both cases, the elastic properties are not constant per pieces but continuous and oscillatory with space. In this work, we did not test any realistic homogenized model such as the Marmousi model in [START_REF] Capdeville | 2-D non-periodic homogenization to upscale elastic media for P-SV waves[END_REF]. Such models make it difficult to draw any general conclusion because results are then strongly dependent on the source and receiver locations. Instead, we tested single wavenumber heterogeneity Efficient high degree spectral element models in order to obtain a reference G as a function of the degree N for three different model roughness. Because realistic geological elastic models amplitude spectra decrease with the wavenumber, these periodic tests can be seen as what can be expected at worst.

The test model λ h = λ min /2 is an upper bound and a realistic geological model is expected to behave more like the λ h = λ min case. In the context of inversion, models are rather smooth and we expect they behave like the λ h = 2λ min case. Finally, let us mention that homogenization cannot remove free-surface or solid-fluid interface topographies, even if they can be homogenized [START_REF] Capdeville | A non-periodic two scale asymptotic method to take account of rough topographies for 2-D elastic wave propagation[END_REF]. This can be a complication, but it can be handled generally by deforming the elements in the vertical direction to match the topography. This comes at the price of a high degree element transformation (where a linear or quadratic element are classically used).

Regarding the number of points per wavelength G, our work confirms that it can be as low as 2.5 using N = 40 for both 2-D and 3-D homogeneous or smooth media. Note that our G is an average number meaning that, for example, it is possible to cast about 16 λ min within one degree 40 element direction. For rough heterogeneous media, G increases more or less following the rule of thumb of equation 11. Compared to a degree 4, a high degree G can be twice smaller. This implies that a factor up to 4 in 2-D and 8 in 3-D in memory can be gained by using a very high degree. This aspect can be particularly interesting in the adjoint inversion context where storing the wavefield on hard disks can be necessary either partially [START_REF] Komatitsch | Anelastic sensitivity kernels with parsimonious storage for adjoint tomography and full waveform inversion[END_REF] or fully to avoid backward propagation [START_REF] Fichtner | Full seismic waveform Efficient high degree spectral element tomography for upper-mantle structure in the australasian region using adjoint methods[END_REF].

It is well known that the GLL points cluster near the edge of the elements for high degrees. This fact a priori implies a dramatically small time step to respect the CFL stability condition. Nevertheless, this effect is not that dramatic because G also decreases for high degrees and this partly compensates the clustering of the GLL points. For example, even if the closest distance between two GLL points is 77 smaller in the reference element for N = 40 compared to N = 4, the time step is only 4.4 times smaller. Then, a low G Efficient high degree spectral element also implies a lower number of degrees of freedom and therefore less computation.

Finally, it is important to note that G(N ) strongly depends on the Poisson's ratio [START_REF] Seriani | Dispersion analysis of spectral element methods for elastic wave propagation[END_REF]. For a fixed degree N , G is degraded (larger) for Poisson's ratio close to 0.5 (in this work, we used a Poisson's ratio of 0.24). [START_REF] Seriani | Dispersion analysis of spectral element methods for elastic wave propagation[END_REF] also have shown that G is less degraded for large degrees than for low degrees. It implies that high degree element is even more interesting compared to low degrees for Poisson's ratio close to 0.5.

CONCLUSIONS

In the end, high-degree SEM is often more effective than lower degrees and when it is not, the losses in computing time is not dramatic whereas the gain in memory is significant. For example, the memory needed for a degree 20 is about a quarter that necessary for a degree 4 in 2-D and about one eighth in 3-D in our numerical experiments. The computation time with degree 12 to 24 can be up to twice faster than the classical degree 4. To conclude, considering using SEM with a significantly higher degree than usual is a good idea in some situations such as full waveform inversion. Efficient high degree spectral element complexity C(N ) scales as

C(N ) ∝ G d+1 (N )N 2 . (A-2)
This formula has a limit: it doesn't account for the fact that, in order to reach a given accuracy, it may happen that the maximum time steps allowed by the CFL cannot be reached.

In that case, SEM numerical complexity C(N ) scales as

C(N ) ∝ G d (N )N .
(A-3) Efficient high degree spectral element

Figure 1 :

 1 Figure 1: Minimum distance dξ N min between two GLL points for the 1-D reference element Λ = [-1, 1] as a function of the degree N with logarithmic scales.

Figure 2 :Figure 3 :

 23 Figure 2: The 2-D homogeneous domain used in this work. The source (black star) is located about 3λ min below the surface. The short and long distance receiver bins (black square) are represented. The kinetic energy snapshot for t = 270 t min and a 60×30 elements mesh, associated to a polynomial degree N = 40, are also displayed.

Figure 4 :

 4 Figure 4: Same as Figure 2 but for the 2-D heterogeneous media with λ h = λ min . The lower left zoom displays a 14.1λ min ×14.1λ min area of V S . The 100×50 elements mesh associated to the degree N = 40 is displayed.

Figure 5

 5 Figure5: 2-D homogeneous model waveform displacement comparison between the analytical solution (gray line), the spectral element simulated solution (dashed line) for the horizontal (top plot) and vertical components (bottom plot). The residual (×10 4 ) is displayed in dotted line. The receiver is a typical receiver from the short distance bin (20 λ min distance).

Figure 6 :

 6 Figure 6: 3-D homogeneous model waveform displacement comparison between the analytical solution (gray line) and the spectral element simulated solution (dashed line) for the horizontal (top plot) and vertical component (bottom plot). The residual (×500) is displayed in dotted line. The receiver is a typical receiver from the short distance bin (20 λ min distance).

Figure 7 :

 7 Figure 7: Error (see equation 9) as a function of G (dotted line) for the short distance receiver bin and for degree N = 40 and a very small dt. The triangles are measured values and the dashed line shows the quadratic interpolation estimation of G to obtain a 1% error. Here, G ≃ 2.35.

Figure 8 :

 8 Figure 8: G as a function of the degree N for the 2-D and 3-D homogeneous models and for the 20 and 200 λ min epicentral distance receiver bins.

Figure 9 :

 9 Figure 9: Heterogeneous models G ratio with the corresponding homogeneous model G as a function of the degree N . Three different heterogeneity roughness are displayed (λ h = {0.5, 1, 2}λ min ).

Figure 10 :

 10 Figure 10: Error (see equation 9) as a function of the time step dt (dotted line), for a large G (spatially over sampled), N =4 in the 2-D homogeneous model and for the long distance receiver bin (200λ min ). The 1% error threshold and its corresponding time step (here dt c = 4.7 × 10 -4 s) are displayed (dashed line).

Figure 11 :Figure 12 :

 1112 Figure11: Error as a function of the time step dt with the optimal G(N ) (in Figure7). Three different degrees are displayed (N = 4, 8 and 40). The optimal dt c = 4.7 × 10 -4 s is displayed (dashed line). For N = 40, the dt CFL is reached before the optimal dt c .

Figure 13 :

 13 Figure 13: Computing time as a function of the degree N for the 2-D homogeneous case, for the long distance bin (200 λ min ) and for the SPEC2DY and SPECFEM2D programs. Symbols with an extra black circle correspond to measurements done with dt = dt CFL .

Figure 14 :

 14 Figure 14: Computing time as a function of the degree N for the 3-D homogeneous case, for the short distance bin (20 λ min ) and for the SEM3D and SPECFEM3D programs. Symbols with an extra black circle correspond to measurements done with dt = dt CFL .

Figure 15 :

 15 Figure 15: Computing time as a function of the degree N for the 2-D heterogeneous cases for the long distance bin (200 λ min ) and for the SPEC2DY program. Symbols with an extra black circle correspond to measurements done with dt = dt CFL .

Table 1 :

 1 The number of N s, nz (number of elements in z direction), minimum dx ratio and maximum dt ratio in 2-D homogeneous half space case for the 200 λ min epicentral distances case. Efficient high degree spectral element

			Information for 200λ min distance
	N	4	8	12	16	20	24	30	40
	nz	505.9 186.2 111.7 79.9 61.9 50.3 39.5 28.9
	dx ref min,4 /dx ref min,N 1	3.45 7.40 12.85 19.82 28.28 43.81 77.20
	G N /G N =4	1	0.662 0.574 0.537 0.514 0.497 0.484 0.469
	dt N =4 /dt	1	1.270 1.634 2.028 2.426 2.812 3.419 4.413

Table Caption Table 1 :

 Caption1 The number of N s, nz (number of elements in z direction), minimum dx ratio and maximum dt ratio in 2-D homogeneous half space case for the 200 λ min epicentral distances case.

ACKNOWLEDGEMENTS

This research was supported by the National Natural Science Fund for Distinguished Young Scholars in China (Grant No.41625016,2017-2021) , the French National Research Agency (ANR-16-CE31-0022-01, 2016-2021) and the China Scholarship Council (File No. 201804910289).

We are very grateful to the computing center CCIPL in Nantes and the server "WUKONG" in Beijing. We thank the two reviewers, E. Chaljub and A. Fichtner for their comments which greatly helped to improve the manuscript. E. Chaljub is at the origin of the appendix of this article. We express our sincere gratitude to their authors for the SPECFEM2D and 3D

APPENDIX A COMPUTATIONAL COMPLEXITY ANALYSIS OF CODE-INDEPENDENT FEATURES FOR SEM

In this appendix, we estimate the theoretical scaling of the numerical complexity of a SEM simulation as a function of the element degree N . The overall numerical complexity of such a simulation is the product of the cost at the element level times the number of elements times the number of time steps. At the element level, the numerical complexity is dominated by the stiffness matrix multiplication, which scales as N d+1 . This is usually where the main optimization effort is mainly put into, following matrix-matrix multiplication strategy [START_REF] Deville | High-order methods for incompressible fluid flow[END_REF]. We have numerically checked that the calculation of the internal forces (the product between elemental stiffness matrix and displacement vector) indeed behaves in N d+1 for the SEM code used in this work. Regarding the number of elements N e , assuming a domain size of L d , we have L d = (N λ λ min ) d , where N λ is the number of the minimum wavelength. We also have λ min = G(N )∆ x where ∆ x the average distance between GLL points within each spectral element. The element size in one direction is L e = N ∆ x and therefore

The number of time steps is controlled by equation 1 and therefore by 1/dx min . Because of the GLL points cluster near the element edges, dx min scales as ∆ x /N (see Figure 1)

and therefore as (N G(N )) -1 . Gathering those estimates, one finds the SEM numerical Efficient high degree spectral element

Figure Captions

•