Nicholas Beck
email: nicholas.beck@hec.ca

Elena Di Bernardino

Mélina Mailhot

Semi-parametric Estimation of Multivariate Extreme Expectiles

Keywords: Dependence, extreme value theory, multivariate risk measures, optimization, semi-parametric estimation 2010 MSC: 60G70, 62H12, 90C53, 91G70

This paper focuses on semi-parametric estimation of multivariate expectiles for extreme levels of risk. Multivariate expectiles and their extremes have been the focus of plentiful research in recent years. In particular, it has been noted that due to the difficulty in estimating these values for elevated levels of risk, an alternative formulation of the underlying optimization problem would be necessary. However, in such a scenario, estimators have only been provided for the limiting cases of tail dependence: independence and comonotonicity. In this paper, we extend the estimation of multivariate extreme expectiles (MEEs) by providing a consistent estimation scheme for random vectors with any arbitrary dependence structure. Specifically, we show that if the upper tail dependence function, tail index, and tail ratio can be consistently estimated, then one would be able to accurately estimate MEEs. The finite-sample performance of this methodology is illustrated using both simulated and real data.

Introduction

Multivariate risk measures

Risk measures are viewed as a crucial tool in a multitude of fields relating to mathematics and statistics. Insurance, finance and economics must all use these measures in evaluating their prospective positions in the face of uncertain outcomes. Moreover, these measures are often required by regulatory organizations so that companies operate in good faith of the interests of stakeholders. To this end, it is easy to see why risk measures are constantly evolving. There have been numerous developments in this field, whether it be through establishing ideal properties [START_REF] Artzner | Coherent measures of risk[END_REF][START_REF] Ekeland | Comonotonic measures of multivariate risks[END_REF][START_REF] Riedel | Dynamic coherent risk measures[END_REF][START_REF] Schmidt | Non-parametric estimation of tail dependence[END_REF], extensions of univariate measures to higher dimension [START_REF] Cossette | Bivariate lower and upper orthant value-at-risk[END_REF][START_REF] Cossette | Vector-valued tail value-at-risk and capital allocation[END_REF][START_REF] Cousin | On multivariate extensions of value-at-risk[END_REF][START_REF] Embrechts | Bounds for functions of multivariate risks[END_REF], estimating these measures non-parametrically [START_REF] Beck | A consistent estimator to the orthant-based tail value-at-risk[END_REF][START_REF] Di Bernardino | Plug-in estimation of level sets in a non-compact setting with applications in multivariate risk theory[END_REF], or even the development of new measures [START_REF] Herrmann | Multivariate geometric expectiles[END_REF][START_REF] Herrmann | Multivariate geometric tail-and range-value-at-risk[END_REF][START_REF] Laniado | Multivariate extremality measure[END_REF][START_REF] Torres | A directional multivariate value at risk[END_REF].

Recently, the concept of elicitability has been developed for risk measures. Gneiting [START_REF] Gneiting | Making and evaluating point forecasts[END_REF] defines elicitability as the ability to express a risk measure in the form of an optimization problem. Specifically, a risk measure T α of a random variable (r.v.) X is considered elicitable if it takes the form

T α (X) = arg min x∈R E {S (x, α)} ,
for a risk level α ∈ (0, 1), where S is the score function associated to the risk measure T α . Elicitability is a desirable property as it allows for a natural backtesting methodology; see Bellini and Bignozzi [START_REF] Bellini | On elicitable risk measures[END_REF] for more information. Two of the most well-known elicitable risk measures are value-at-risk (VaR) and expectiles. From the usual definition in the univariate setting, VaR is the smallest such value a r.v. can take for which the probability of failing to exceed it is at least α. Thus, at level α, VaR is classically written VaR α (X) = inf{x ∈ R : F(x) ≥ α} but can be written in an elicitable form as VaR α (X) = arg min x∈R E{α(Xx) + + (1 -α)(Xx) -}, where x + = max{0, x} and x -= max{0, -x}. While elicitable, it is well documented that VaR is not a coherent risk measure in the sense of Artzner et al. [START_REF] Artzner | Coherent measures of risk[END_REF] as it does not respect the crucial sub-additivity property. The second measure, and the focus of this paper, is the α-level expectile. Expectiles e α (X) have been introduced by Newey and Powell [START_REF] Newey | Asymmetric least squares estimation and testing[END_REF] as the minimizers of an asymmetric quadratic loss. Indeed, at level α ∈ (0, 1) the expectile of a r.v. X is defined by

e α (X) = arg min x∈R E α(X -x) 2 + + (1 -α)(X -x) 2 -. (1)
Note that the uniqueness of the minimum in ( 1) is guaranteed by the strict convexity of the scoring function. When α = 1/2, it is well known that e α (X) = E[X], thus expectiles can be seen as an asymmetric generalization of the mean. The term "expectiles" has probably been suggested as a combination of "expectation" and "quantiles". For α < 1/2, expectiles are super-additive and therefore not coherent. However, counter to VaR, expectiles are coherent when α ≥ 0.5; see Ziegel [START_REF] Ziegel | Coherence and elicitability[END_REF]. As risk managers are often interested in elevated levels of risk α, this is quite advantageous.

Thus, the potential coherence of expectiles makes them an alternative risk measure to the VaR for practical applications. It should be noted that, prior to the recent discussions on expectiles, a third measure of risk had been discussed as a potential replacement for VaR in practical applications: expected shortfall (ES) [START_REF] Acerbi | Expected shortfall: A natural coherent alternative to value at risk[END_REF][START_REF] Yamai | On the validity of value-at-risk: Comparative analyses with expected shortfall[END_REF][START_REF] Yamai | Value-at-risk versus expected shortfall: A practical perspective[END_REF]. The expected shortfall at risk level α for a r.v. X is defined by

ES α (X) = 1 1 -α 1 α
VaR u (X)du, which can be interpreted as the expected value of X given that it surpasses VaR α (X). While expected shortfall is a coherent risk measure [START_REF] Acerbi | On the coherence of expected shortfall[END_REF], it is not elicitable. Given that expectiles can be both elicitable and coherent, they present themselves as a potential improvement over both VaR and ES. In fact, for α > 1/2, expectiles are the only risk measure that is coherent, elicitable and law-invariant [START_REF] Ziegel | Coherence and elicitability[END_REF].

Expectiles are uniquely identified by the first-order condition, α E {Xe α (X)} + = (1 -α) E {Xe α (X)} -.

The above equation can also be written as

1 -α α = E {X -e α (X)} + E {X -e α (X)} - .
This makes the economic interpretation of expectiles as risk measures clearer. Expectiles can be seen as the value of X that provides a profits/loss ratio of (1 -α)/α. A comparison between expectiles and more common financial risk measures, including VaR and ES can be found in Bellini and Di Bernardino [START_REF] Bellini | Risk management with expectiles[END_REF]. Note that both expectiles and VaR fall into the family of generalized quantiles, as described in Bellini et al. [START_REF] Bellini | Generalized quantiles as risk measures[END_REF], defined by

q α (X) = arg min x∈R (αE [Φ 1 {(X -x) + }] + (1 -α)E [Φ 2 {(X -x) -}]) ,
where Φ 1 , Φ 2 : [0, ∞) → [0, ∞) are strictly increasing convex functions satisfying Φ i (0) = 0 and Φ i (1) = 1 for i ∈ {1, 2}.

While univariate risk measures are well researched and developed, they limit one's ability for sound inference in real life situations as they often ignore any potential dependence that risks can have on one another. Ignoring such relationships can provide inaccurate inference and induce prohibitive losses for businesses. As such, our interest lies in exploring multivariate expectiles as these dependence structures can be incorporated directly into the measure.

First, the notion of elicitability is extended to the multivariate context. Here, for any d-dimensional random vector X ∈ R d an associated risk measure T α is elicitable if we can write it as

T α (X) = arg min x∈R d {S(x, α)} .
Now consider the multivariate extensions of expectiles. Herrmann et al. [START_REF] Herrmann | Multivariate geometric expectiles[END_REF] present a multivariate geometric definition of expectiles, which considers the possibility of selecting different levels of risk α for the components of the vector X. Maume-Deschamps et al. [START_REF] Maume-Deschamps | Multivariate extensions of expectiles risk measures[END_REF] define two notions of multivariate expectiles: L p -expectiles and Σ-expectiles.

The focus of this paper is on L p -expectiles. For any p ∈ [1, ∞], define the L p -expectile of a random vector X by e p α (X) = arg min

x∈R d E α (X -x) + 2 p + (1 -α) (X -x) - 2 p .
Specifically, this paper explores semi-parametric estimation of the L 1 -expectile, i.e., e α (X) := e 1 α (X) = arg min

x∈R d E          α         d i=1 |X i -x i | +         2 + (1 -α)         d i=1 |X i -x i | -         2          . ( 2
)
Analogously to the univariate case, the

L 1 -expectile in (2) is the unique solution in R d of 1 -α α = E[ (X -x) + 1 1{X k > x k }] E[ (X -x) -1 1{X k < x k }] , k ∈ {1, . . . , d}.
Thus, it can be interpreted as a ratio of expected positive scenarios over negative ones.

In terms of estimation, we provide a consistent estimator for e α (X) in (2) for elevated risk levels α ≈ 1 when the underlying dependence structure and marginal distributions are unknown, through the approximated optimization problem arg min

Θ∈R d L Λ (Θ) ,
for some loss function L and consistently estimated parameter set Λ. In Maume-Deschamps et al. [START_REF] Maume-Deschamps | Multivariate extensions of expectiles risk measures[END_REF], it was shown that multivariate expectiles could be consistently estimated using Robbins-Monro's stochastic optimization for moderate levels of α (see discussion in Section 4.3). However, for elevated levels of α their approach fails. As such, extra care must be given when estimating e α (X) at extreme levels of risk.

The outline of the paper is as follows: In the remainder of Section 1, multivariate extreme expectiles (MEEs) and some important related concepts are introduced. In Section 2, a general optimization framework is provided for MEEs. The main contribution of this paper is in Section 3, where the consistency of the proposed approximate loss function and its capability to accurately produce MEEs for extreme levels of risk α is provided (Theorem 1 and Corollary 3, respectively). In Section 4, the finite-sample performance of our procedure is illustrated with several simulation studies. In Section 5, possible improvements of our estimation procedure are discussed. Section 6 provides concluding remarks. Finally, in the supplementary material Section 7 we provide additional simulation, exploration, and discussion for our procedure (see Section 7.1), an application of our procedure for real data (see Section 7.2) and the gradient calculation of the considered loss function (see Section 7.3).

Multivariate Extreme Expectiles

The notion of extreme expectiles has been explored in the univariate case by Bellini and Bignozzi [START_REF] Bellini | On elicitable risk measures[END_REF], Mao et al. [START_REF] Mao | Asymptotic expansions of generalized quantiles and expectiles for extreme risks[END_REF], Bellini and Di Bernardino [START_REF] Bellini | Risk management with expectiles[END_REF], Daouia et al. [START_REF] Daouia | Estimation of tail risk based on extreme expectiles[END_REF][START_REF] Daouia | Extreme M-quantiles as risk measures: From L 1 to L p optimization[END_REF] and in the multivariate case by Maume-Deschamps et al. [START_REF] Maume-Deschamps | Multivariate extensions of expectiles risk measures[END_REF][START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF] and Herrmann et al. [START_REF] Herrmann | Multivariate geometric expectiles[END_REF], for example. In the latter works, the authors discuss situations for which expectiles can be solved for explicitly or for moderate risk levels α. We recall several asymptotic and multivariate concepts that will be crucial to this work. The first of which is the notion of multivariate regular variation (MRV). Recall in one dimension a function f is considered to be regularly varying at a ∈ {0, +∞} with index θ ∈ R,

denoted f ∈ RV θ (a), if lim t→a f (tx) f (t) = x θ
for all x > 0. Extending this definition to random vectors, one has the following equivalent definitions:

Definition 1 (MRV equivalent definitions). Let X be a random vector on R d , the following definitions are equivalent:

• The vector X has regularly varying tail of index θ.

• There exists for all x > 0 a finite measure µ on the unit sphere

S d-1 , a normalizing function b : (0, ∞) → (0, ∞) such that lim t→+∞ Pr X > xb(t), X X ∈ • = x -θ µ(•).
The measure µ depends on the chosen norm, it is called the spectral measure of X.

• There exists a finite measure µ on the unit sphere S d-1 , a slowly varying function L, and a positive real θ > 0 such that

lim x→+∞ x θ L(x) Pr X > x, X X ∈ B = µ(B)
for all B ∈ B(S d-1 ) with µ(∂B) = 0, where B represents the Borel σ-algebra and ∂B the boundary of the set B.

For more information on the concept of regular variation and how it relates to MEEs, see Section 1 of Maume-Deschamps et al. [START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF]. Drawing further from Maume-Deschamps et al. [START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF], the focus here is on marginal distributions with asymptotically equivalent tails. This can be stated in the following assumption.

Assumption 1 (Equivalent Tails). Assume that X has equivalent regularly varying marginal tails, which means that both (i) and (ii) below are satisfied:

(i) F 1 ∈ RV -θ (+∞) with θ > 0;
(ii) The tails of X 1 , . . . , X d are equivalent. That is for all i ∈ {2, . . . , d}, there is a positive constant c i such that

lim x→+∞ F i (x) F 1 (x) = c i .
Conditions (i) and (ii) imply that all margins are RV -θ (∞). Denote c i the tail ratio. Obviously c 1 = 1. Extreme asymptotic relationships within a random vector X on R d are further quantified by the upper tail dependence function defined by

λ(x) = lim t→0 + Pr {F 1 (X 1 ) ≥ 1 -tx 1 , . . . , F d (X d ) ≥ 1 -tx d } t = lim t→0 + C(tx 1 , . . . , tx d ) t , ( 3
)
where C is the survival copula associated to X (see, for instance, Equation (1.4) or (1.5) in Maume-Deschamps et al. [START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF]). Pairwise upper tail dependencies will be of particular importance in this paper, so we denote by λ ik (x i , x k ), the upper tail dependence function for the random pair (X i , X k ). It is also beneficial to define the stable tail dependence function (stdf)

L(x) = lim t→0 + 1 t Pr{F(X 1 ) > 1 -tx 1 or . . . or F(X d ) > 1 -tx d }
and the bivariate relationship L(x, y) = x + y -λ(x, y). The following proposition from Maume-Deschamps et al. [START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF] connects all these concepts to MEEs by providing an alternative system of equations when α → 1.

Proposition 1 (Maume-Deschamps et al. [START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF], Proposition 3.3). Assume that X has regular varying multivariate distribution in the sense of Definition 1 and that Assumption 1 holds. Consider the L 1 -expectile e α (X) = (e i α (X)) i∈{1,...,d} . Then any limit vector Θ = (η, β 2 , . . . , β d ) of (1 -α)/F 1 {e 1 α (X)}, e 2 α (X)/e 1 α (X), . . . , e d α (X)/e 1 α (X) for α → 1 satisfies the following system of equations

1 θ -1 -η β θ k c k = - d i=1,i k        ∞ β i β k λ ik c i c k t -θ , 1 dt -η β θ-1 k c k β i        , k ∈ {1, . . . , d}. (4)
In particular, Maume-Deschamps et al. [START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF] show that the system solutions are

η ⊥ = 1 (θ -1) 1 + d j=2 c 1/(θ-1) j , β ⊥ k = c 1/(θ-1) k , (5)
for asymptotically independent random vectors and

η + = 1 (θ -1) , β + k = c 1/θ k , (6)
for comonotonic random vectors, for all k ∈ {1, . . . , d}.

Denote Θ ⊥ = (η ⊥ , β ⊥ 2 , . . . , β ⊥ d ) and Θ + = (η + , β + 2 , . . . , β + d ).
In the rest of this paper we are interested in using Proposition 1 towards estimating Θ and using this result to calculate MEEs. Given these definitions, it is easy to see that estimating e α (X) for extreme risk levels α will require one to first consistently estimate θ, c = (c 2 , . . . , c d ) and λ ik .

Optimization problem for MEE

As mentioned in Section 1, for moderate levels of α, effective calculation of e α (X) is possible through Robbins-Monro's stochastic optimization. However, such a method does not perform well for α ≈ 1 (see Section 4.3 for further details). An alternative methodology will be required. In Definition 2 below, we use Proposition 1 as inspiration to define an alternative optimization problem for MEEs. Definition 2. Consider a d-dimensional random vector X with regularly varying multivariate distribution in the sense of Definition 1 with margins possessing equivalent tails as given in Assumption 1. Let β * = (β 1 , . . . , β d ) with β 1 = 1 and Θ = (η, β 2 , . . . , β d ). Let c 1 = 1 and Λ = (θ, c 2 , . . . , c d , λ(•)), where c i is the tail ratio for i ∈ {2, . . . , d}, θ is the tail index and λ is the upper tail dependence function associated to the random vector X in (3). Define the loss function

L Λ (Θ) := 1 2 F Λ (Θ) 2 2 , (7)
where

F Λ (Θ) = F (1) Λ (Θ), . . . , F (d) Λ (Θ) = g (1) Λ (Θ) + f (1) Λ (Θ), . . . , g (d) Λ (Θ) + f (d) Λ (Θ) ,
with, for all k ∈ {1, . . . , d},

g (k) Λ (Θ) = 1 θ -1 -η β θ k c k and f (k) Λ (Θ) = d i k        ∞ β i β k λ ik c i c k t -θ , 1 dt -η β θ-1 k c k β i       
.

Define an optimal vector Θ * , obtained by optimizing the loss function L Λ in [START_REF] Bellini | Generalized quantiles as risk measures[END_REF], i.e.,

Θ * = arg min Θ L Λ (Θ). ( 8
)
Corollary 1 establishes that the optimization problem in (8) can be used to calculate MEEs.

Corollary 1. Consider a d-dimensional random vector X with regularly varying multivariate distribution in the sense of Definition 1. Consider the L 1 -expectile e α (X) = (e i α (X)) i∈{1,...,d} as in [START_REF] Acerbi | On the coherence of expected shortfall[END_REF]. Under Assumption 1 and the assumption that the vector (1 -α)/F 1 {e 1 α (X)}, e 2 α (X)/e 1 α (X), . . . , e d α (X)/e 1 α (X) has a unique limit point Θ = (η, β 2 , . . . , β d ) for α → 1, then Θ satisfies the optimization problem (8) in Definition 2. Moreover, the optimization problems (2) and (8) are equivalent for α → 1, i.e., lim α→1 e i α (X) e i α (X)

= 1, for i ∈ {1, . . . , d}, e α (X) := VaR α (X 1 )η 1/θ (1, β 2 , . . . , β d ) .

Proof. From Proposition 1 and the uniqueness of the limit Θ = (η, β 2 , . . . , β d ) we get that Θ satisfies the optimization problem in [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF]. Furthermore, from Proposition 5.1 in Maume-Deschamps et al. [START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF], we know that e α (X) ∼ VaR α (X 1 )η 1/θ (1, β 2 , . . . , β d ) , for α → 1. Hence, the result is proved.

To solve the optimization problem (8) in Definition 2, the Broyden-Fletcher-Goldfarb-Shanno (BFGS) descent algorithm will be used in the present paper. The algorithm is a member of the family of quasi-Newton optimization methods and allows one to avoid calculating second derivatives, a feat which can drastically improve computation time. The BFGS algorithm is given in Algorithm 1.

Remark 1. For convergence results of the BFGS algorithm the interested reader is referred to Nocedal and Wright [START_REF] Nocedal | Numerical Optimization[END_REF]. It is important to note that, in the general case, the true optimization problem [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF] does not satisfy the convexity assumption required for global convergence of the BFGS algorithm. Nevertheless, the choice to use BFGS is motivated by it's simplicity of implementation and further validated by the satisfactory numerical results provided in Section 4. More sophisticated algorithms which extend BFGS would indeed be possible, including L-BFGS-B [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF] and modified BFGS [START_REF] Li | A modified BFGS method and its global convergence in nonconvex minimization[END_REF][START_REF] Li | On the global convergence of the BFGS method for nonconvex unconstrained optimization problems[END_REF], the latter of which allows one to relax the convexity assumption and still achieve global convergence. These methods are briefly discussed in Section 5.

Algorithm 1 BFGS quasi-Newton optimization algorithm (Step 0) Put counter k := 0 and choose initial values Θ 0 ∈ R d , H 0 ∈ R d×d initial approximation to the inverse of the Hessian matrix of L Λ , σ ∈ (0, 1/2), ρ ∈ (σ, 1), and ≥ 0.

(Step 1) Let L Λ as in Definition 2. If ∇L Λ Θ k ≤ : STOP. (Step 2) Calculate the direction d k = -H k ∇L Λ Θ k . ( Step 3
) Determine the step size t k > 0 such that the Wolfe conditions are satisfied, i.e.,

L Λ Θ k + t k d k ≤ L Λ Θ k + σt k ∇L Λ Θ k , ∇L Λ Θ k + t k d k d k ≥ ρ∇L Λ Θ k d k .
(Step 4) Update the following:

• Θ k+1 := Θ k + t k d k • y k := ∇L Λ Θ k+1 -∇L Λ Θ k , • s k := Θ k+1 -Θ k • H k+1 := I -ρ k s k y k H k I -ρ k s k y k + ρ k s k s k , where ρ k = 1/y k s k . ( Step

Approximated Optimization Problem

Since the underlying distribution is not known in practice, crucial information, including the values of λ, θ, and c i for i ∈ {2, . . . , d}, is not available. As such, direct application of the BFGS algorithm for the optimization problem ( 8) is not possible. Instead, an approximated version of (8) will be solved, i.e., one can focus on the approximate optimum arg min

Θ∈R d L Λ(Θ) (9)
for some vector of estimators Λ = ( θ, ĉ2 , . . . , ĉd , λ). Specifically, convergence of the estimated optimum can be shown in the following way. The first step will be to show that Λ Pr -→ n→∞ Λ, and subsequently that

L Λ(Θ) Pr -→ n→∞ L Λ (Θ) and ∇L Λ(Θ) Pr -→ n→∞ ∇L Λ (Θ)
, where Pr -→ n→∞ denotes convergence in probability. Then, we will prove the consistency of every iteration of the BFGS algorithm (see Theorem 1)

Θ k Pr -→ n→∞ Θ k , k ∈ {1, 2, . . .}.
Finally, by using a two-step procedure (see Algorithm 2), we can study the behaviour of the proposed method (see Corollary 3).

Loss function

Establishing the consistency of every iteration of the BFGS algorithm applied to the approximated loss function comes down to establishing the convergence of the approximate loss function L Λ and its gradient ∇L Λ. In the end, this will require convergence arguments for estimators λik , θ, and ĉi for i ∈ {2, . . . , d}, of λ ik , θ, and c i for i ∈ {2, . . . , d}, respectively. Moreover, the convergence

∞ β i β k λik c i c k t -θ , 1 dt Pr -→ n→∞ ∞ β i β k λ ik c i c k t -θ , 1 dt
is also required. For n observations of the random vector X, denote the random sample X i = (X i1 , . . . , X id ), i ∈ {1, . . . , n}. Next, denote X j:1,n ≤ X j:2,n ≤ • • • ≤ X j:n,n the order statistics associated to the realizations of the j th random variable of X for j ∈ {1, . . . , d}. First, consider the estimated tail index θ via Hill's estimator Hill [START_REF] Hill | A simple general approach to inference about the tail of a distribution[END_REF],

γ = 1 θ θ i=1 ln(X 1:n-i+1,n ) -ln(X 1:n-θ ,n ), θ = 1 γ , (10)
for an intermediate integer sequence

, i.e., θ = θ (n) such that θ (n) = o(n) and θ (n) -→ n→∞ ∞. The tail ratio, c i , can be estimated by ĉi = X i:n-i +1,n X 1:n-i +1,n - θ , i ∈ {2, . . . , d}, (11)
where i = i (n) is an intermediate integer sequence and θ is as defined in [START_REF] Chaouch | Design-based estimation for geometric quantiles with application to outlier detection[END_REF]. The consistency of θ and ĉi is

Upper tail dependence function

One group of estimators for the upper tail dependence function take the form

λ(x) = n λ C λ n x ,
where C is some estimator of the survival copula C for the random vector X and λ = λ (n) is an intermediate integer sequence. Thus, to estimate the upper tail dependence function effectively, a consistent estimator of the copula is required. In this paper, the empirical beta copula was chosen. For a random sample X 1 , . . . , X n with corresponding ranks denoted by R i, j = n k=1 1{X k j ≤ X i j }, i ∈ {1, . . . , n} and j ∈ {1, . . . , d}, the empirical beta copula is defined as

C Beta n (u) = 1 n n i=1 d j=1 F R i, j ,n (u i ),
where u ∈ [0, 1] d and F r,n (u) is the cumulative distribution function (cdf) of U r:n , the r th order statistic of n independent uniform random variables on [0, 1]. This also corresponds to U r:n ∼ Beta(r, n + 1r). In two dimensions this gives

C Beta n (u, v) = 1 n n i=1 F R i,1 ,n (u)F R i,2 ,n (v).
Moreover, as C Beta n is a genuine copula, using the bivariate relationship

C(u, v) = u + v -1 + C(1 -u, 1 -v) gives pairwise upper tail dependence function λik Beta (x i , x k ) = n λ          λ n x i + λ n x k -1 + 1 n n j=1 F R j,i ,n 1 - λ n x i F R j,k ,n 1 - λ n x k          . ( 12
)
Proposition 2 establishes the consistency of the upper tail dependence function based on the empirical beta copula.

Proposition 2. For the estimator of the upper tail dependence function based on the empirical beta copula, as defined in [START_REF] Cossette | Vector-valued tail value-at-risk and capital allocation[END_REF], it holds that λik Beta is a consistent estimator for λ ik , i.e.,

λik Beta (x i , x k ) Pr -→ n→∞ λ ik (x i , x k ).
Proof. The proof is straightforward from the work of Kiriliouk et al. [START_REF] Kiriliouk | An estimator of the stable tail dependence function based on the empirical beta copula[END_REF]. They define the empirical process

B Beta n, λ = λ (L Beta n, λ -L),
where

L Beta n, λ (x) = n 1 -C Beta n (1 -λ x/n) / λ
is the stdf defined under the empirical beta copula. In particular, they established the weak convergence (as defined in van der Vaart and Wellner [START_REF] Van Der Vaart | Weak Convergence and Empirical Processes[END_REF])

B Beta n, λ B on ∞ ([0, 1]
), where B is some stochastic process with continuous trajectories. Using this result, and the bivariate relationship that for any stdf λ(x, y) = x + y -L(x, y), the empirical process for upper tail dependence function based on the empirical beta copula will have the same (up to a sign) asymptotic distribution. This implies that

λik Beta (x i , x k ) Pr -→ n→∞ λ ik (x i , x k ),
establishing the result.

100

In particular, Proposition 2 gives n (0.5, 1.5) (red), and the empirical beta λBeta (0.5, 1.5) (green), for varying values of λ with sample size ranging from 100 to 10000. The true value, λ 23 (0.5, 1.5) (black), is also included. The margins are Pareto type I. Specifically,

λik Beta c i c k t -θ , 1 Pr -→ n→∞ λ ik c i c k t -θ , 1 .
X i ∼ P{2, 1.25(1 + i)}, i ∈ {1, 2, 3}.
In Fig. 1, an illustration of the performance of λBeta is provided using independent Pareto type I distributions. Specifically, we consider

X i ∼ P{2, 1.25(1 + i)}, i ∈ {1, 2, 3} with cdf F i (x) = 1 -{1.25(1 + i)/x} 2 .
For comparison, the pairwise upper tail dependence function based on the standard empirical copula

λik n (x i , x k ) = 1 λ n j=1 1 F in (X ji ) > 1 - λ n x i , F jn (X jk ) > 1 - λ n x k ,
where X ji is the j th realization of the r.v. X i , is also used. As can be seen, λBeta provides a far more accurate estimate of λ and doesn't exhibit any of the oscillating behavior that is seen when considering λn . The performance is also better at smaller sample sizes when a larger limiting subsequence is considered.

Integral of the upper tail dependence function

The final piece of the estimation procedure is to approximate

∞ β i β k λ ik c i c k t -θ , 1 dt. ( 13
)
To our knowledge, a consistent estimator for the integral of the upper tail dependence function has not been established 105 in prior literature. Again, the underlying upper tail dependence function is estimated using the empirical beta copula.

The result is established in the following proposition.

Proposition 3. Consider the estimator of λ as proposed in [START_REF] Cossette | Vector-valued tail value-at-risk and capital allocation[END_REF], then

∞ β i β k λik Beta c i c k t -θ , 1 dt Pr -→ n→∞ ∞ β i β k λ ik c i c k t -θ , 1 dt.
Proof. The proof boils down to an application of the Dominated Convergence Theorem. First, the convergence of the integrand follows from Proposition 2. Second, as the empirical beta copula is a proper copula, it is easy to see that its upper tail dependence function satisfies the bounds

λ ik,⊥ c i c k t -θ , 1 ≤ λik Beta c i c k t -θ , 1 ≤ λ ik,+ c i c k t -θ , 1 and hence 0 ≤ λik Beta c i c k t -θ , 1 ≤ min c i c k t -θ , 1 .
See for instance Proposition 2.1 in Kiriliouk et al. [START_REF] Kiriliouk | An estimator of the stable tail dependence function based on the empirical beta copula[END_REF]. Third, it is shown in Maume-Deschamps et al. [START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF] (see the proof of Lemma 3.4) that

∞ β i β k λ ik,+ c i c k t -θ , 1 dt = c i c k β i β k -θ+1                β k β i c k c i -1 θ -1        + c k c i β k β i -θ + 1 θ -1         1 +        β i β k c k c i -1 θ -1        +         -θ+1          < ∞.
Thus, by dominated convergence, we get the result.

In Fig. 2, an illustration of the finite sample performance of the estimator is provided. In particular, for each sample size and subsequence combination, the calculation of

∞ β 2 β 3 λ 23 Beta c 2 c 3 t -θ
, 1 dt was replicated 500 times. As can be seen 110 for samples of size n > 400 and subsequences n q for q > 0.4, the estimator performs well.

Once again, it is compared to an equivalent quantity, using the standard empirical copula. In fact, under the empirical copula, a closed form estimator for (13) exists, as outlined in Proposition 4 below.

Proposition 4. Under the empirical copula c 3 t -θ , 1 dt for various sample sizes and subsequences. The variable q refers to the power used for our subsequence, i.e., λ = n q , q ∈ {0.1, 0.2, . . . , 0.9} with sample size n ∈ {100, 200, . . . , 1000, 2500, 5000}. They are compared to the true value

C n (u) = 1 n n i=1 d j=1 1 R i, j n ≤ u j , it holds that ∞ β i β k λik n c i c k t -θ , 1 dt = 1 λ n j=1         n -R i, j + 1 λ • c k c i -1 θ - β i β k         + • 1 X k j ≥ X k:n-λ ,n .
∞ β 2 β 3 λ 23 c 2 c 3 t -θ
, 1 dt under Pareto type I margins with a comonotonic dependence structure (green horizontal line). Specifically,

X i ∼ P{2, 1.25(1 + i)}, i ∈ {1, 2, 3}. Proof. One can write ∞ β i β k λik n c i c k t -θ , 1 dt = ∞ β i β k 1 λ n j=1 1        1 -F n,i (X i j ) < λ c i c k t -θ n , 1 -F n,k (X k j ) < λ n        dt = ∞ β i β k 1 λ n j=1 1        n -R i, j + 1 n < λ c i c k t -θ n        1 X k j ≥ X k:n-λ ,n dt = 1 λ n j=1 ∞ β i β k 1        n -R i, j + 1 λ • c k c i -1 θ > t        1 X k j ≥ X k:n-λ ,n dt = 1 λ n j=1 n-R i, j +1 λ • c k c i -1 θ β i β k 1        n -R i, j + 1 λ • c k c i -1 θ ≥ β i β k        1 X k j ≥ X k:n-λ ,n dt = 1 λ n j=1         n -R i, j + 1 λ • c k c i -1 θ - β i β k         1        n -R i, j + 1 λ • c k c i -1 θ ≥ β i β k        1 X k j ≥ X k:n-λ ,n = 1 λ n j=1         n -R i, j + 1 λ • c k c i -1 θ - β i β k         + 1 X k j ≥ X k:n-λ ,n .
Fig. 3 illustrates the performance of

∞ β 2 /β 3 λn (c 2 /c 3 )t -θ , 1 dt as in Proposition 4.
Comparing to Fig. 2, one sees that 115 the empirical beta copula far surpasses the standard empirical copula in terms of accuracy.

While the convergence of each individual estimator has been established, they must now be combined to establish the convergence of our estimation procedure. Using the consistency of ( 10), [START_REF] Cossette | Bivariate lower and upper orthant value-at-risk[END_REF], and (12) yields the following corollary Corollary 2. Taking Λ = ( θ, ĉ2 , . . . , ĉd , λik Beta ) as defined in [START_REF] Chaouch | Design-based estimation for geometric quantiles with application to outlier detection[END_REF], [START_REF] Cossette | Bivariate lower and upper orthant value-at-risk[END_REF], and (12), one has

∞ β i β k λik Beta ĉi ĉk t -θ, 1 dt Pr -→ n→∞ ∞ β i β k λ ik c i c k t -θ , 1 dt.
Proof. For simplicity, we establish a similar result for λik Beta , i.e.,

λik Beta ĉi ĉk t -θ, 1 Pr -→ n→∞ λ ik c i c k t -θ , 1 .
The proof makes use of the continuous mapping theorem (denoted CMT in the following); see, e.g., van der Vaart [START_REF] Van Der | Asymptotic Statistics[END_REF], Theorem 2.3. Beginning with the first statement, using several iterations of the triangle inequality yields c 3 t -θ , 1 dt as in Proposition 4 for various sample sizes and subsequences. The variable q refers to the power used for our subsequence, i.e., λ = n q , q ∈ {0.1, 0.2, . . . , 0.9} with sample size n ∈ {100, 200, . . . , 1000, 2500, 5000}. They are compared to the true value

∞ β 2 β 3 λ 23 c 2 c 3 t -θ
, 1 dt under Pareto type I margins with a comonotonic dependence structure (green horizontal line). Specifically,

X i ∼ P{2, 1.25(1 + i)}, i ∈ {1, 2, 3}. λik Beta ĉi ĉk t -θ, 1 -λ ik c i c k t -θ , 1 = λik Beta ĉi ĉk t -θ, 1 -λ ik ĉi ĉk t -θ, 1 + λ ik ĉi ĉk t -θ, 1 -λ ik c i c k t -θ , 1 = λik Beta ĉi ĉk t -θ, 1 -λ ik ĉi ĉk t -θ, 1 + λ ik ĉi ĉk t -θ, 1 -λ ik c i c k t -θ, 1 + λ ik c i c k t -θ, 1 -λ ik c i c k t -θ , 1 = λik Beta ĉi ĉk t -θ, 1 -λ ik ĉi ĉk t -θ, 1 + λ ik ĉi ĉk t -θ, 1 -λ ik c i ĉk t -θ, 1 + λ ik c i ĉk t -θ, 1 -λ ik c i c k t -θ, 1 + λ ik c i c k t -θ, 1 -λ ik c i c k t -θ , 1 ≤ λik Beta ĉi ĉk t -θ, 1 -λ ik ĉi ĉk t -θ, 1 (i) + λ ik ĉi ĉk t -θ, 1 -λ ik c i ĉk t -θ, 1 (ii)
+ λ ik c i ĉk t -θ, 1 -λ ik c i c k t -θ, 1 (ii)
+ λ ik c i c k t -θ, 1 -λ ik c i c k t -θ , 1 (iii) Pr -→ n→∞ 0,
where (i) by the convergence of λik Beta (Proposition 2), (ii) by the CMT on λ ik and the consistency of ĉi , i ∈ {1, . . . , d} (Maume-Deschamps et al. [START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF], Proposition 4.2), and (iii) by the CMT on λ ik and the consistency of θ (de Haan and Ferreira [START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF], Theorem 3.2.2).

The proof of the main statement would follow in the same way, or as a result of the dominated convergence theorem, similar to Proposition 3.

For an illustration of Corollary 2, simulations are presented in Fig. 4. It is clear that adding extra uncertainty through calculation of θ and ĉi , i ∈ {2, . . . , d}, decreases the accuracy of the full estimator. The variability of these quantities is well documented and this would naturally require an increased number of observations to estimate with confidence. Moreover, the selection of the subsequences θ , λ and i , i ∈ {2, . . . , d} is crucial. For instance, the impact of the choice of λ is illustrated in Figs. 234. In a real data study, the subsequences could be calibrated by some stability criterion. More precisely, in Section 7.2, we will use the selection procedure of intermediate sequences described, for instance, by Cai et al. [START_REF] Cai | Estimation of the marginal expected shortfall: the mean when a related variable is extreme[END_REF] and Di Bernardino and Prieur [START_REF] Di Bernardino | Estimation of the multivariate conditional tail expectation for extreme risk levels: Illustration on environmental data sets[END_REF]. We will consider a range of possibilities for intermediate sequences for which the associated estimators showed some stability. Then, to gain robustness, we average the estimates corresponding to the selected ranges.

Consistency

In Theorem 1 below, we first prove that the approximated loss function (approximated gradient function) in the optimization problem [START_REF] Cai | Estimation of the marginal expected shortfall: the mean when a related variable is extreme[END_REF] converges to the loss function in (8) (to the gradient function). Furthermore, every iteration of the BFGS algorithm applied to the approximated optimization problem is consistent for the true optimization problem. Subsequently, we show that, under the assumption that BFGS captures the global minimum of the true optimization problem [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF], the approximate optimization problem (9) can be used to solve for MEEs (see Algorithm 2 and Corollary 3 below). ĉ3 t -θ, 1 dt for various sample sizes and subsequences. The variable q refers to the power used for our subsequence, i.e., λ = n q , q ∈ {0.1, 0.2, . . . , 0.9} with sample size n ∈ {100, 200, . . . , 1000, 2500, 5000}. They are compared to the true value

∞ β 2 β 3 λ 23 c 2 c 3 t -θ
, 1 dt under Pareto type I margins with a comonotonic dependence structure (green horizontal line). Specifically,

X i ∼ P{2, 1.25(1 + i)}, i ∈ {1, 2, 3}.
Moreover, given identical starting values Θ 0 , H 0 , σ ∈ (0, 1/2), ρ ∈ (σ, 1) and ≥ 0, for any step k, it holds that

Θ k Pr -→ n→∞ Θ k . ( 15
)
Proof. Noting the structure of L Λ and ∇L Λ (Θ) (detailed in 7.3), the consistency results in ( 14) are a direct result of the consistency of θ, ĉi , i ∈ {2, . . . , d}, Propositions 2, 3 and Slutzky's Theorem. For [START_REF] Daouia | On kernel smoothing for extremal quantile regression[END_REF], we use induction. First, fix Θ 0 , H 0 , σ ∈ (0, 1/2), ρ ∈ (σ, 1) and ≥ 0. Recalling the iteration step of BFGS algorithm, we have

d0 = -H 0 ∇L Λ(Θ 0 ) Pr -→ n→∞ -H 0 ∇L Λ (Θ 0 ) = d 0 .
Similarly, it is clear that t0

Pr -→ n→∞ t 0 . It follows that

Θ 1 = Θ 0 + t0 d0 Pr -→ n→∞ Θ 0 + t 0 d 0 = Θ 1 ,
Θ k+1 = Θ k + tk dk = Θ k + tk H -1 k ∇L Λ Θ k Pr -→ n→∞ Θ k + t k H -1 k ∇L Λ Θ k = Θ k+1 .
Thus, for any step k ∈ {1, 2, . . .}, Θ k Pr -→ n→∞ Θ k by induction.

We now proceed by using an iterated two-step procedure. We first provide an adequate estimate of the true loss function (and its gradient) and then proceed with the optimization procedure. This is formalized in Algorithm 2.

Algorithm 2 Two-steps estimation procedure for MEEs (Step 1) Taking the limit n → ∞. Establish the consistency of Λ and subsequently L Λ and ∇L Λ (see Theorem 1).

(Step 2) Taking the limit k → ∞.

Optimize the consistently approximated problem from Step 1 using the BFGS algorithm.

145

The resulting asymptotic behavior of our estimation scheme using Algorithm 2 is established in Corollary 3. It is important to note the limit order in the non-exchangeable iterated limit [START_REF] Daouia | Estimation of tail risk based on extreme expectiles[END_REF].

Corollary 3. Under the assumptions of Corollaries 1, 2 and Theorem 1, and further assuming that the BFGS algorithm solves for the global minimum of (8), it holds that

lim k→∞ (plim n→∞ Θ k ) = Θ * . ( 16
)
Proof. Under the assumption that BFGS solves for the global minimum of (8), i.e., Θ k -→ k→∞ Θ * , we have

lim k→∞ lim n→∞ Pr Θ k -Θ * ≥ = lim k→∞ lim n→∞ Pr Θ k -Θ k + Θ k -Θ * ≥ ≤ lim k→∞ lim n→∞ Pr Θ k -Θ k + Θ k -Θ * ≥ ≤ lim k→∞ lim n→∞ Pr Θ k -Θ k ≥ /2
Step 1

+ Pr Θ k -Θ * ≥ /2
Step 2 = 0, where firstly (Step 1) converges to zero as n → ∞ by Theorem 1 and subsequently (Step 2), which does not depend on n, converges to 0 as k → ∞ by assumption.

Numerical study

In this section, several numerical examples are provided demonstrating the effectiveness of our estimation procedure. In Section 4.1, we focus on the limiting cases of tail dependence: asymptotic independence and comonotonicity. As per Maume-Deschamps et al. [START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF], closed form solutions exist for these tail dependence structures so they will provide a good measuring stick for our methodology. In Section 4.2, we provide examples using non-trivial tail dependence structures, i.e., neither asymptotic independence nor comonotonicity, where closed-form solutions for the system of equations in (4) exist. As such, we can further illustrate the performance of our estimates in this intermediate situation.

Results in each case are provided for both Θ and êα (X). In the latter case, using Proposition 5.1 from Maume-Deschamps et al. [START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF], the final calculation of MEE estimates can be done using êα (X) = VaR α (X 1 ) η1/ θ(1, β2 , . . . , βd ) ,

where

VaR α (X 1 ) = X 1:n-k(n)+1,n k(n) (1 -α)n 1/ θ , (17)
called Weissman's estimator [START_REF] Weissman | Estimation of parameters and large quantiles based on the k largest observations[END_REF], for an intermediate integer sequence k(n) is better suited than standard estimation techniques, using the empirical cdf, for large risk levels α ≈ 1. While both Θ are êα (X) are provided, results for Θ will more effectively demonstrate the capabilities of our procedure as they do not rely on the performance of ( 17).

The simulations were conducted for sample sizes n ∈ {50, 100, 250, 500, 1000, 2500} and repeated 500 times for each sample size. The subsequences for θ, c i for i ∈ {2, . . . , d} and λ were θ = n 0.75 , i = n 0.75 for i ∈ {2, . . . , d} and λ = n 0.5 respectively. The risk level was set to α = 1 -1/5000 = 0.9998 to reflect the largest sample size considered. True system solutions Θ can be found in the last line of each table. First, we consider the limiting cases of tail dependence: asymptotic independence and comonotonicity.

Asymptotic independence and comonotonicity

For the limiting cases of tail dependence, we consider three and five dimensional random vectors with Pareto type I marginals. Specifically, X i ∼ P(b i , a) where a = θ = 3.5 and b i = 1.25(i + 1), i ∈ {1, . . . , d}, d ∈ {3, 5}. For Pareto margins with cdfs F i , i ∈ {1, . . . , 5} and these parameter settings, the tail ratios are

c i = lim x→∞ F i (x) F 1 (x) = b i b 1 a =                    1.000, i = 1, 4.134, i = 2, 11.314, i = 3, 24.705, i = 4, 46.765, i = 5.

Three-dimensional case

The summarized results of the simulations are presented in Table 1 and boxplots can be found in Fig. 5. In the threedimensional case, it can be seen that the estimation procedure proposed in this paper performs quite well. At lower samples, the results are only slightly biased, but with a fair amount of standard error. As the sample size increases however, the algorithm is able to successfully estimate the true value of Θ and subsequently the true value of e α (X)

for α = 0.9998 with minimal standard errors. It should be noted that for the comonotonic simulations the slight bias and small standard errors would mean that the true value would likely not be contained in any reasonable confidence interval, though the differences become smaller and smaller as the sample size increases. Fig. 5: Boxplots for the estimated solution vector for simulations using a three-dimensional random vector with Pareto type I margins for various sample sizes n. Specifically, X i ∼ P{3.5, 1.25(1 + i)}, i ∈ {1, 2, 3} with θ = i = n 0.75 and λ = n 0.50 . Results are for independent margins (left) and comonotonic margins (right). True values Θ ⊥ and Θ + (dashed green lines) provided for comparison.

Five-dimensional case

The five-dimensional case is summarized in Table 2. It can be seen that the algorithm scales well to the increased dimension. This could be explained in part by the fact we are only ever dealing with pairwise dependencies through the pairwise upper tail dependence function λ ik . The results presented focus on the independent five-dimensional case, though results for comonotonic random vectors (not included here) were very similar.

Pareto-Clayton model

In this section we present a simulation study using a non-trivial tail dependence structure. Specifically, we consider a random vector X with Pareto type I margins, as defined previously, whose dependence structure is defined by a survival Clayton copula. The interested reader is referred to Cuberos et al. [START_REF] Cuberos | Copulas checker-type approximations: Application to quantiles estimation of sums of dependent random variables[END_REF], Maume-Deschamps et al. [START_REF] Maume-Deschamps | Impact of dependence on some multivariate risk indicators[END_REF], Yeh [START_REF] Yeh | The frailty and the Archimedean structure of the general multivariate Pareto distributions[END_REF]. Recall that for a pair of standard uniform rv's (U, V) with a survival Clayton copula, we get

C(u, v; θ C ) = max u -θ C + v -θ C -1, 0 -1/θ C ,
where the parameter θ C ∈ [-1, ∞) \ {0} determines the strength of the dependence in the upper tail. In particular, the pairwise upper tail dependence function can now be defined as

λ(x, y; θ C ) = x -θ C + y -θ C -1/θ C .
For a general development of MEEs for Archimedean copulas, see Maume-Deschamps et al. [START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF]. Notice that for this Pareto-Clayton model, it follows that For simplicity, and to facilitate closed form solutions to the system of equations in Proposition 1, we consider the case where θ C = θ -1 . Then, the integral of λ ik can be written as

λ c i c k t -θ , 1; θ C =        c i c k -θ C t -θ θ C + 1        -1/θ C .
∞ β i β k λ c i c k t -θ , 1; θ -1 dt = 1 θ -1 c i c k        β i β k + c i c k 1/θ        -θ+1
, giving a system of equations

1 θ -1 -η β θ k c k = - i k         1 θ -1 c i c k        β i β k + c i c k 1/θ        -θ+1 -η β θ-1 k c k β i         , k ∈ {1, . . . , d}.
We now consider two cases: (i) d = 2 , θ = 2, c 2 > 0, and (ii) d ≥ 2, θ > 1, c i = 1, i ∈ {1, . . . , d}. In case (i) and (ii), the closed-form MEEs system solutions are given by

Θ =        1 + c 2 (c 3/4 2 + √ c 2 ) -1 c 3/4 2 + 1 , c 3/4 2        and Θ = (d -1)2 -θ+1 + 1 d(θ -1) , 1, . . . , 1 ,
respectively. In the latter case, this solution makes heuristic sense. Having equivalent margins with a unit tail ratio would imply that they are identically distributed. As such, one would expect that the components of the MEEs are identical and solely based on the dependence parameter θ. In Table 3, we present the results of our simulation for each case. Here we see that performance is again satisfactory. In both cases our estimation procedure is able to accurately capture the true value of MEEs. An additional study of MEEs for the Pareto-Clayton model can be found in the Supplementary material Section (see Section 7.1). In particular, we study their behavior for varying θ C .

(b) Case (ii), with X i ∼ P(2, 2.5), i ∈ {1, 2, 3}. This implies d = 3, θ C = 1/2, and c i = 1, i ∈ {1, 2, 3} n η β2 β3 ê1 α (X) ê2 α (X) ê3 α (

Comparison with alternative methods

To further substantiate the performance of our estimation procedure, we compare the results in Table 1 to simulations using the fully parametric estimators

η⊥ = 1 ( θ -1) 1 + d j=2 ĉ 1 θ-1 j , β⊥ k = ĉ 1 θ-1 k and η+ = 1 θ -1 , β+ k = ĉ 1 θ k (18)
for k ∈ {2, . . . , d}, based on ( 5) and ( 6), respectively, from Maume-Deschamps et al. [START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF]. Using these fully parametric estimators for the system solutions, estimators for MEEs can be written for independent and comonotonic random vectors, respectively, where VaR α (X 1 ) is the Weissman's estimator [START_REF] Weissman | Estimation of parameters and large quantiles based on the k largest observations[END_REF] in [START_REF] Daouia | Extreme M-quantiles as risk measures: From L 1 to L p optimization[END_REF]. Analogously to Table 1, in Table 4 we consider median parameter estimates from the simulation study of independent and comonotonic Pareto margins in three dimensions. Standard deviations for the estimates are presented in parentheses. Here, as in Section 4.1,

X i ∼ P{3.5, 1.25(1 + i)}, i ∈ {1, 2, 3}.
Unsurprisingly, the fully parametric estimators perform better in particular for smaller sample sizes and have smaller standard errors. However, our semi-parametric method (see Table 1) does not require any pre-imposed dependence assumption and still provides satisfactory results. Note that the standard errors being zero for β2 and β3 for the comonotonic case are a result of the estimators structure and the distributional assumptions, i.e.,

β+ i = ĉ 1 θ i = X i:n-i +1,n X i:n-i +1,n = F -1 i (U n-i +1,n ) F -1 1 (U n-i +1,n ) = b i b 1
for i ∈ {2, 3}, where U n-i +1,n is the ni + 1, nth order statistic from the standard uniform sample used to generate the data. For comonotonic random vectors, each realization of the copula is generated using the same uniform realization for each component, i.e., X ji = F -1 j (U i ), j ∈ {1, . . . , d}, i ∈ {1, . . . , n}.

In the full non-parametric framework, Maume-Deschamps et al. [START_REF] Maume-Deschamps | Multivariate extensions of expectiles risk measures[END_REF] use Robbins-Monro's stochastic optimization to solve the system of equations

α d i=1 E [(X i -x i ) + 1{X k > x k }] = (1 -α) d i=1 E [(X i -x i ) -1{X k < x k }] , k ∈ {1, . . . , d}.
For a comparison with our results under the same independent Pareto model, the interested reader is referred to Fig. 7 in Maume-Deschamps et al. [START_REF] Maume-Deschamps | Multivariate extensions of expectiles risk measures[END_REF] (page 40). It is easy to see that for elevated levels of α, the non-parametric approach solved through Robbins-Monro's stochastic optimization performs quite poorly.

Table 4: Medians of results from MEE simulation study for independent (top) and comonotonic (bottom) random vectors using the full parametric estimators for Θ and e α in [START_REF] Deheuvels | Almost sure convergence of the Hill estimator[END_REF] with Pareto type I margins and various sample sizes n. Specifically, n ∈ {50, 100, 250, 500, 1000, 2500, 5000} and X i ∼ P{3.5, 1.25(1 + i)}, i ∈ {1, 2, 3} with θ = i = n 0.75 , λ = n 0.50 . Standard deviations for the estimates are presented in parentheses.

(a) Independent Margins.

n η β2 β3 ê1 α (X) ê2 α (X) ê3 α (

Discussion

As mentioned in Section 2, in the general case the loss function in ( 7) could be non-convex, and thus, might not satisfy the assumptions required for global convergence of the BFGS algorithm. Nevertheless, we chose this method for its relative simplicity and ease of implementation, as well as the satisfactory numerical results it has provided for our problem. To build a more rigorous estimation procedure, one could consider implementing methods which build upon the foundation of the original BFGS algorithm. For instance, Li and Fukushima [START_REF] Li | A modified BFGS method and its global convergence in nonconvex minimization[END_REF][START_REF] Li | On the global convergence of the BFGS method for nonconvex unconstrained optimization problems[END_REF] discuss modifications to the BFGS algorithm which allow one to relax the assumption of convexity and still ensure global convergence. Specifically, they discuss the importance of Lipschitz continuity and using a damping method on the iterates of the approximate inverse Hessian matrix. Equivalently, Yuan et al. [START_REF] Yuan | Global convergence of BFGS and PRP methods under a modified weak Wolfe-Powell line search[END_REF] discuss modifications to the Wolfe conditions in Algorithm 1 which ensures global convergence of BFGS for general functions. Such methods could further improve our results. However, while BFGS is already included in many R packages, the above methods are not readily available. Moreover, manually programming them and doing so efficiently would require a specific work which is beyond the scope of the present work.

Additionally, while the BFGS algorithm works well for our problem, it was designed for unconstrained optimization problems. Clearly for extreme multivariate expectiles it is required that Θ > 0 componentwise. A further improvement to this problem would be incorporating a box-constrained methodology for solving [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF] and, more importantly, [START_REF] Cai | Estimation of the marginal expected shortfall: the mean when a related variable is extreme[END_REF]. One such choice would be box-constrained BFGS algorithm (or BFGS-B). This algorithm would allow us to incorporate properties of MEEs. For instance, as comonotonicity and independence represent the limiting cases of dependence in the upper tail, all expectiles will be bounded above and below by comonotonic and independent expectiles. Moreover, from Lemmas 2.4 and 2.5 Maume-Deschamps et al. [START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF] we know that for asymptotically independent and comonotonic random vectors the solutions to (8) are Θ ⊥ and Θ + defined by ( 5) and ( 6), respectively.

As these situations represent the extreme cases of MEEs, this allows one to bound any solution Θ * of ( 8) to the d-

dimensional box Θ L 1 , Θ U 1 × • • • × Θ L d , Θ U d where Θ L i = min Θ ⊥ i , Θ + i and Θ U i = max Θ ⊥ i , Θ + i for i ∈ {1, .
. . , d}. However, because the estimation is done without prior knowledge of the underlying distribution, the box-constraints in the algorithm must also be estimated. To this end, the bounds of each box will be slightly relaxed. In other words, the algorithm would look for a solution in

Θ L 1 , Θ U 1 × • • • × Θ L d , Θ U d ,
where

Θ L i = (1 -γ L ) min Θ ⊥ i , Θ + i and Θ U i = (1 + γ U ) max Θ ⊥ i , Θ + i for some γ L , γ U > 0.
Another change that could be made to the general BFGS algorithm is to incorporate limited memory storage of the inverse hessian H k . This would prove most beneficial when the dimension of the problem is large and we would want to limit the number of operations on large inverse Hessian matrices. Combining these two modifications gives the limited-memory box-constrained BFGS algorithm (L-BFGS-B). Preliminary simulations using L-BFGS-B algorithm applied to [START_REF] Cai | Estimation of the marginal expected shortfall: the mean when a related variable is extreme[END_REF] were conducted and the results were encouraging. Indeed, incorporating this methodology could provide further benefits in terms of both accuracy and computational efficiency. For more information on BFGS, L-BFGS or L-BFGS-B algorithms see Nocedal and Wright [START_REF] Nocedal | Numerical Optimization[END_REF], Liu and Nocedal [START_REF] Liu | On the limited memory BFGS method for large scale optimization[END_REF] and Byrd et al. [START_REF] Byrd | A limited memory algorithm for bound constrained optimization[END_REF]. Incorporating these changes along with the modifications described in Li and Fukushima [START_REF] Li | A modified BFGS method and its global convergence in nonconvex minimization[END_REF][START_REF] Li | On the global convergence of the BFGS method for nonconvex unconstrained optimization problems[END_REF] and Yuan et al. [START_REF] Yuan | Global convergence of BFGS and PRP methods under a modified weak Wolfe-Powell line search[END_REF] would also be interesting.

Conclusions

In this paper, we have presented an effective semi-parametric method for accurately estimating multivariate expectiles for extreme levels of risk α ≈ 1. Using the empirical beta copula, taking advantage of the fact that it is a genuine copula, we are able to establish the consistency of the approximate loss function in (9) when using empirical estimators for the tail index, tail ratio and upper tail dependence function. We have also proposed a new consistent estimator for the integral of the upper tail dependence function. Numerical simulations demonstrate the satisfactory finite-sample performance of our estimation scheme when compared to the true MEEs stemming from a variety of dependence structures, where closed form solutions exist. Furthermore, the set-up provided by the system of equations in Proposition 1 allows us to mitigate the curse of dimensionality as we are only ever required to consider pairwise dependencies. As such, the number of operations required for each iteration of the optimization algorithm will be d 2 .

In the future, it would be interesting to compare the performance of the algorithm for a variety of estimators. For instance, if one were to consider the checkerboard copula instead of the empirical beta copula. Finally, while finite-sample results are encouraging, the estimation procedure could be further streamlined by considering an algorithm which could take advantage of additional properties of our optimization problem. For example, L-BFGS-B algorithm could improve estimation by bounding solutions to (9) by the MEEs under the limiting cases of upper tail dependence.

Furthermore the present work could open interesting perspectives in functional statistics field. Indeed, the univariate well known conditional expectile/quantile model (see, e.g., Daouia et al. [START_REF] Daouia | On kernel smoothing for extremal quantile regression[END_REF], Usseglio-Carleve [START_REF] Usseglio-Carleve | Estimation of conditional extreme risk measures from heavy-tailed elliptical random vectors[END_REF], Girard et al. [START_REF] Girard | Nonparametric extreme conditional expectile estimation[END_REF]) can be adapted to our semi-parametric estimation scheme for new conditional multidimensional L 1 -expectile from heavy-tailed distributions when functional covariate information is available. The interested reader is also referred for instance to Gardes et al. [START_REF] Gardes | Functional nonparametric estimation of conditional extreme quantiles[END_REF], Gardes and Stupfler [START_REF] Gardes | An integrated functional Weissman estimator for conditional extreme quantiles[END_REF] where functional non-parametric estimators for conditional extreme quantiles are proposed.

More precisely, let (X i , Z i ), i = 1, . . . , n, be n independent copies of a random pair (X, Z) ∈ R d × E where (E, d) is a not necessarily finite-dimensional Polish space endowed with a semi-metric d. For instance, E can be the standard p-dimensional space R p , a space of continuous functions over a compact metric space, or a Lebesgue space L p (R), to name a few. From Equation (2), a future development of this work would be to consider the functional conditional multidimensional L 1 -expectile extension given by e 1 α (X, z) := arg min

x∈R d E          α         d i=1 |X i -x i | +         2 + (1 -α)         d i=1 |X i -x i | -         2 Z = z         
.

Then the study in the present paper could be adapted to estimate extreme e 1 α (X, z) by using the extrapolation technique when α → 1 together with a nonparametric estimation of the conditional tail copula (see, e.g., Gardes and Girard [START_REF] Gardes | Nonparametric estimation of the conditional tail copula[END_REF]). The main objective of this semi-parametric approach for extreme functional e 1 α (X, z) would be to balance the trade-off between the high sensitivity to dimension of the nonparametric models and the relative lack of flexibility of the completely parametric models (see, e.g., Goia and Vieu [START_REF] Goia | An introduction to recent advances in high/infinite dimensional statistics[END_REF]).

Supplementary material

Pareto-Clayton model: additional studies

Here we present the results of several extra simulations demonstrating the impact of dependence parameter θ C on our results. Note that in general θ C ∈ [-1, +∞) \ {0} where the cases θ C = -1, θ C → 0 and θ C → ∞ correspond to countermonotonic, independent and comonotonic dependence structure, respectively. As such, we test our methodology by considering very small and large values of θ C so as to see if we can recapture these limiting results presented in Section 4.1. Furthermore, we consider a range of values in between the limiting cases to study the behavior of MEEs as the dependence becomes stronger.

Once more, we consider X i ∼ P(b i , a) where a = θ = 3.5 and b i = 1.25(i + 1), i ∈ {1, 2, 3}. Recall that for these settings the system solutions for independent and comonotonic random variables are Θ ⊥ = (0.0740, 1.764, 2.639), Θ + = (0.400, 1.500, 2.000), respectively (see Section 4.1). For all simulations in this section, we consider a survival Clayton copula with dependence parameter θ C ∈ {0.01, 0.1, 0.25, 0.5, 0.75, 1, 2.5, 5, 7.5, 10, 15, 20, 100}. Finally, θ = i = n 0.75 and λ = n 0.50 . In Table 5, it can be seen that our estimation procedure accurately estimates the true system solutions for Pareto-Clayton models which are arbitrarily close to independence and comonotonicity. Moreover, in Fig. 6, we plot the boxplots for β2 from our simulations. We also provide the true values under independence and comonotonicity for comparison (green and red dotted lines in Fig. 6, respectively). Here, we are able to extract more insight into the behavior of MEEs as the dependence parameter shifts. For the Pareto-Clayton model, we see that the parameter estimates seem to decrease monotonically from those representing independence to those representing comonotonicity.

n η β2 β3 ê1 α (X) ê2 α (X) ê3 α (

Data study

Here we present an application of MEEs to a real data set. In particular, we consider the Pima Indians Diabetes Database (https://www.kaggle.com/uciml/pima-indians-diabetes-database). This data set has been used in previous studies of multivariate risk measures, see, e.g., Girard and Stupfler [START_REF] Girard | Extreme geometric quantiles in a multivariate regular variation framework[END_REF], Chaouch and Goga [START_REF] Chaouch | Design-based estimation for geometric quantiles with application to outlier detection[END_REF]. The data set consists of 9 variables (Pregnancies, Glucose, Diastolic Blood Pressure (DBP), Skin Thickness, Insulin Levels, body mass index (BMI), Diabetes Pedigree Function, Age and Outcome) for n = 768 individuals from the Pima Indian population. The study was conducted in an attempt to understand the prevalence of diabetes in this population.

We consider the two-dimensional case for pairs (X i1 , X i2 ) where X i1 is the BMI and X i2 is the DBP for individual i ∈ {1, . . . , 768}.

The first step is to estimate θ and ĉ2 . Using the selection procedure described, for instance, by Cai et al. [START_REF] Cai | Estimation of the marginal expected shortfall: the mean when a related variable is extreme[END_REF], we consider a range of possibilities for intermediate sequences θ and 2 for which the estimators showed some stability.

Then, to gain robustness, we average the estimates corresponding to the selected ranges. The interested reader is also referred to Di Bernardino and Prieur [START_REF] Di Bernardino | Estimation of the multivariate conditional tail expectation for extreme risk levels: Illustration on environmental data sets[END_REF]. For θ this turned out to be values of θ ∈ [100, 150]. The corresponding Hill plots are presented in Fig. 7. Considering these values the tail index was estimated to be θ = 9.126. The same procedure is followed for ĉ2 and the resulting estimated value is ĉ2 = 1103.046. In Fig. 8, the resulting expectiles are plotted. Recalling the first order conditions of L 1 expectiles, e α (X) is the value of x such that

d i=1 E [(X i -x i ) + 1{X k > x k }] d i=1 E [(X i -x i ) -1{X k < x k }] = 1 -α α , k ∈ {1, 2}.
This could be interpreted as the value of x which fixes the ratio of participation in positive scenarios to negative scenarios at (1 -α)/α for all variables. Consider the MEE at level α = 1 -1/768 = 0.9987, we have ê0.9987 (X) = (45.433, 106.493). As was noted in similar studies of this data set, e.g., Chaouch and Goga [START_REF] Chaouch | Design-based estimation for geometric quantiles with application to outlier detection[END_REF], an individual with a BMI above 31.7 kg/m 2 would qualify as severely overweight and a DBP value above 95 mm Hg would be at risk of hypertension. Both these factors are considered to impact the prevalence of diabetes. First, note that ∇L Λ (•) is given by

∇L Λ (•) = ∇         1 2 d k=1 g (k) Λ (•) + f (k) Λ (•) 2         =                                    d k=1 g (k) Λ (•) + f (k) Λ (•) × ∂ ∂η g (k) Λ (•) + f (k) Λ (•) d k=1 g (k) Λ (•) + f (k) Λ (•) × ∂ ∂β 2 g (k) Λ (•) + f (k) Λ (•) . . . d k=1 g (k) Λ (•) + f (k) Λ (•) × ∂ ∂β d g (k) Λ (•) + f (k) Λ (•)                                    . (19)
Then, the marginal derivatives in [START_REF] Di Bernardino | Plug-in estimation of level sets in a non-compact setting with applications in multivariate risk theory[END_REF] can be written explicitly evaluated

∂ ∂Θ j g (k) Λ (Θ) =              - β θ k c k , j = 1,
-ηθ

β θ-1 k c k , j = k 1, 0, otherwise, ∂ ∂Θ j f (k) Λ (Θ) =                  - β θ-1 k c k i k β i , j = 1, i k λ ik ci c k β i β k -θ , 1 β i β 2 k -η (θ-1)β θ-2 k c k β i , j = k 1, -λ jk c j c k β j β k -θ , 1 1 β k -η β θ-1 k c k , otherwise,
j, k ∈ {1, . . . , d}.

5 )

 5 Set k ← k + 1 and go to (Step 1).

  Limiting subsequence λ = n 0.25 Limiting subsequence λ = n 0.40 Limiting subsequence λ = n 0.75

Fig. 1 :

 1 Fig. 1: Comparison between estimates of the upper tail dependence function based the empirical copula, λ23n (0.5, 1.5) (red), and the empirical beta λBeta (0.5, 1.5) (green), for varying values of λ with sample size ranging from 100 to 10000. The true value, λ 23 (0.5, 1.5) (black), is also included. The margins are Pareto type I. Specifically, X i ∼ P{2, 1.25(1 + i)}, i ∈ {1, 2, 3}.

  Fig. 2: Performance of ∞ β 2 β 3 λ23 Beta

  Fig. 3: Performance of ∞ β 2 β 3 λ23 n c 2

Theorem 1 .

 1 Fig. 4: Performance of ∞ β 2 β 3 λ23 Beta

( a )

 a Simulation results for Pareto-Clayton model with θ C = 0.01. The final row Θ represents the known solutions under independence.

Fig. 6 :

 6 Fig. 6: Simulation results for the estimate of β2 under the Pareto-Clayton model with varying dependence parameter θ C and Pareto type I margins. Specifically, X i ∼ P{3.5, 1.25(1 + i)}, i ∈ {1, 2, 3} with θ = i = n 0.75 and λ = n 0.50 . Here the sample size is fixed at n = 5000 and θ C ∈ {0.01, 0.1, 0.25, 0.5, 0.75, 1, 2.5, 5, 7.5, 10, 15, 20}. Dotted lines provide true values for asymptotic independence (green) and comonotonicity (red) with β ⊥ 2 = 1.764 and β + 2 = 1.5, respectively.

Fig. 7 :

 7 Fig. 7: Hill Plots for estimating θ using BMI measurements. We consider a range of possibilities for intermediate sequences θ for which the estimator θ shows some stability. In left panel we consider a wide range where θ ∈ [10, 250] and the right panel shows the chosen estimation range where θ ∈ [100, 150].

Fig. 8 :

 8 Fig.8: Extreme expectiles (ê 1 α (X), ê2 α (X)) for pairs of BMI and DBP in the Pima Indians Diabetes Study for several risk values α ∈ [0.9, 1) (black dots). The extreme expectiles for α = 1 -1/n ≈ 0.9987,(45.433, 106.493), are represented by red triangles.

7. 3 .

 3 Gradient ∇L Λ (•) in Algorithm 1Recall that β * = (β 1 , . . . , β d ) with β 1 = 1 and Θ = (η, β 2 , . . . , β d ). Then, with our notation, Θ j = η, j = 1 β j , j ∈ {2, . . . , d} .

Table 1 :

 1 Medians of results from MEE simulation study for Pareto type I margins in three dimensions under independence (top) and comonotonicity (bottom) for various sample sizes n. Specifically, n ∈ {50, 100, 250, 500, 1000, 2500, 5000} and X i ∼ P{3.5, 1.25(1 + i)}, i ∈ {1, 2, 3} with θ = i = n 0.75 and λ = n 0.50 . Standard deviations for the estimates are presented in parentheses.

	(a) Independent Margins

Table 2 :

 2 Medians of results from MEE simulation study for Pareto type I margins in five dimensions linked with the independence copula for various sample sizes n. Specifically, n ∈ {50, 100, 250, 500, 1000, 2500, 5000} and X i ∼ P{3.5, 1.25(1 + i)}, i ∈ {1, . . . , 5} with θ = i = n 0.75 , λ = n 0.50 . Standard deviations for the estimates are presented in parentheses.

			(a) Parameter estimates.
	n	η	β2	β3	β4	β5
	50	0.040 (0.058) 1.779 (0.773) 2.672 (1.159) 3.641 (1.405) 4.662 (1.887)
	100	0.037 (0.033) 1.760 (0.404) 2.622 (0.602) 3.553 (0.851) 4.565 (1.136)
	250	0.032 (0.017) 1.779 (0.288) 2.674 (0.410) 3.667 (0.547) 4.783 (0.757)
	500	0.032 (0.012) 1.783 (0.205) 2.668 (0.300) 3.661 (0.410) 4.690 (0.545)
	1000 0.031 (0.008) 1.782 (0.142) 2.653 (0.206) 3.661 (0.286) 4.723 (0.386)
	2500 0.031 (0.005) 1.769 (0.109) 2.653 (0.148) 3.638 (0.210) 4.695 (0.282)
	5000 0.030 (0.004) 1.763 (0.073) 2.645 (0.108) 3.618 (0.154) 4.667 (0.205)
	Θ	0.029	1.764	2.639	3.607	4.655
			(b) Expectiles estimates.
	n	ê1 α (X)	ê2 α (X)	ê3 α (X)	ê4 α (X)	ê5 α (X)
	50	12.301 (7.676) 21.755 (22.48) 33.493 (31.78) 46.092 (43.99) 58.649 (53.44)
	100	11.546 (4.261) 19.608 (8.535) 29.990 (12.93) 41.265 (19.67) 53.211 (26.50)
	250	10.968 (2.510) 19.893 (5.057) 29.810 (7.971) 41.149 (11.05) 52.947 (15.67)
	500	10.761 (1.870) 19.108 (3.689) 28.967 (5.726) 38.747 (8.207) 50.317 (11.57)
	1000	10.653 (1.261) 19.185 (2.324) 28.377 (3.698) 38.982 (5.352) 50.274 (7.274)
	2500	10.660 (0.784) 18.776 (1.591) 28.107 (2.491) 38.583 (3.622) 49.900 (4.956)
	5000	10.581 (0.628) 18.653 (1.108) 27.999 (1.823) 38.329 (2.632) 49.508 (3.671)
	e 0.9998 (X)	10.390	18.330	27.420	37.475	48.372

Table 3 :

 3 Medians of results from MEE simulation study for Pareto type I margins linked with a survival Clayton copula for different dimensions d, dependence levels θ C , and various sample sizes n. Specifically, n ∈ {50, 100, 250, 500, 1000, 2500, 5000} with θ = i = n 0.75 and λ = n 0.50 . Standard deviations for the estimates are presented in parentheses. Case (i), with X i ∼ P{2, 1.25(1 + i)}, i ∈ {1, 2}. This implies d = 2, θ C = 1/2 and c 2 = 1.25.

	n	η	β2	ê1 α (X)	ê2 α (X)
	50	0.468 (0.802) 1.873 (9.347) 140.525 (1241.8) 304.536 (5744.1)
	100 0.480 (0.645) 1.872 (1.175) 125.333 (678.19) 276.406 (826.40)
	250 0.508 (0.399) 1.895 (1.249) 120.870 (189.65) 251.578 (332.57)
	500 0.520 (0.313) 1.927 (0.776) 123.610 (99.534) 253.192 (179.82)
	1000 0.580 (0.255) 1.848 (0.620) 126.894 (59.084) 247.724 (93.983)
	2500 0.574 (0.213) 1.883 (0.517) 128.624 (42.478) 252.158 (64.426)
	5000 0.575 (0.184) 1.857 (0.491) 129.376 (32.306) 250.832 (49.495)
	Θ	0.590	1.837	135.798	249.477

(a)

Table 5 :

 5 Medians of results from MEE simulation study for Pareto type I margins using a Clayton survival copula with dependence parameter θ

C = 0.01 (top) and θ C = 100 (bottom) for various sample sizes. Specifically, n ∈ {50, 100, 250, 500, 1000, 2500, 5000} and X i ∼ P{3.5, 1.25(1 + i)}, i ∈ {1, 2, 3} with θ = i = n 0.75 and λ = n 0.50 . Standard deviations for the estimates are presented in parentheses.

  .6 43.6 42.8 41.3 40.1 39.3 38.6 38.0 37.5 36.8 36.4 35.7 35.4 34.9 .6 39.4 39.3 39.2 39.0 38.9 38.7 38.5 38.4 38.3 38.1 38.0

	alpha (CI, p =0.95)	4 6 8 10 12 14 16 18	10 22 34 46 58 70 82 94 107 122 137 152 167 182 197 212 227 242 49.6 45Order Statistics	alpha (CI, p =0.95)	8 9 10 11	100 39.8 39Order Statistics 104 108 112 116 120 124 128	132	136	140	144	148
			(a) θ ∈ [10, 250]			(b)

θ ∈ [100, 150]

established, for instance, in Deheuvels et al.[START_REF] Deheuvels | Almost sure convergence of the Hill estimator[END_REF] or de Haan and Ferreira[START_REF] De Haan | Extreme Value Theory: An Introduction[END_REF] and Maume-Deschamps et al.[START_REF] Maume-Deschamps | Extremes for multivariate expectiles[END_REF], respectively. This leaves the upper tail dependence function as the sole remaining piece.

Acknowledgments:

The authors express their gratitude to the Associate Editor and the three anonymous Referees for their valuable comments which have helped to greatly improve the quality of this manuscript. The authors would also like to thank Professor Christian Genest for providing several insightful comments during the earlier stages of this research. This work has been partially supported by the SIMONS Foundation, Society of Actuaries James C. Hickman Scholar Program and the Centre de recherche mathématiques de Montréal through the Simons-CRM Program. This work has been supported by the French government, through the 3IA Côte d'Azur Investments in the Future project managed by the National Research Agency (ANR) with the reference number ANR-19-P3IA-0002.