Identifying groups of variables with the potential of being large simultaneously
Maël Chiapino Anne Sabourin Johan Segers











Maël Chiapino Anne Sabourin LTCI, Télécom ParisTech, Université Paris-Saclay 
 46, rue Barrault, 75013 Paris, France.

              mael.chiapino@telecom-paristech.fr , anne.sabourin@telecom-paristech.fr           

           Johan Segers Université catholique de Louvain, Institut de Statistique, Biostatistique et Sciences Actuarielles, 
Voie du Roman Pays 20, B-1348 Louvain-la-Neuve, Belgium.

          johan.segers@uclouvain.be






Identifying groups of variables with the potential of being large simultaneously
  
  



    Maël Chiapino Anne Sabourin
Johan Segers

        

    Received: date / Accepted: date
=========================================================================================









  Identifying groups of variables that may be large simultaneously amounts to finding out which joint tail dependence coefficients of a multivariate distribution are positive. The asymptotic distribution of a vector of nonparametric, rank-based estimators of these coefficients justifies a stopping criterion in an algorithm that searches the collection of all possible groups of variables in a systematic way, from smaller groups to larger ones. The issue that the tolerance level in the stopping criterion should depend on the size of the groups is circumvented by the use of a conditional tail dependence coefficient. Alternatively, such stopping criteria can be based on limit distributions of rank-based estimators of the coefficient of tail dependence, quantifying the speed of decay of joint survival functions. Numerical experiments indicate that the algorithm's effectiveness for detecting tail-dependent groups of variables is highest when paired with a criterion based on a Hill-type estimator of the coefficient of tail dependence.

  
  MSC 62G32 MSC 62H15 MSC 62H05 MSC 62H20








































§ INTRODUCTION


  

A question that often arises when monitoring several variables is which groups of variables are prone to be large simultaneously. In food risk management, for instance, the variables under consideration may be the concentrations of different contaminants in blood samples of consumers. In environmental applications, one may be interested in several physical variables such as wind speed and precipitation recorded at several locations, with the purpose of setting off a regional warning when several of these variables exceed a high threshold. In the context of semi-supervised anomaly detection, when the training sample is mostly made of normal instances, identifying the groups of variables which are likely to be large together allows to label certain new instances as abnormal.

The latter use case is the motivation behind the DAMEX algorithm <cit.>. In a regular variation framework, identifying those groups among d variables that may be large simultaneously amounts to identifying the support of the exponent measure. The algorithm returns the list of groups of features α⊂{1,…, d} such that the mass of the empirical exponent measure on certain cones exceeds a user-defined threshold. However, when the empirical version of the exponent measure is scattered over a large number of such cones, the DAMEX algorithm does not discover a clear-cut structure. <cit.> encounter this difficulty for extreme streamflow data recorded at several locations of the French river system.

To overcome this issue, the same authors come up with the CLEF (CLustering Extreme Features) algorithm. Instead of partitioning the sample space, CLEF considers nested regions corresponding to increasing subsets of components. A group of variables is enlarged until there is no longer enough evidence that all features in it may be large together. In this respect, CLEF resembles the Apriori algorithm <cit.>, which is a data-mining tool for discovering maximal sets of items among d available items that are frequently  bought together by consumers. Apriori considers increasing itemsets that are made to grow until their frequency falls below a user-defined threshold. In CLEF, the stopping criterion concerns the relative frequency of simultaneous occurrences of large values of all components in a considered subset compared to the frequency of simultaneous occurrences of larges values of all but one component in this subset. <cit.> find the method to work well on real and simulated data but do not investigate the asymptotic properties of the statistic underlying the stopping criterion.

Our contributions are three-fold. First, we investigate the asymptotic behavior of the statistic underlying CLEF. In this way, the informal stopping criterion can be turned into a proper hypothesis test with controllable level. A second issue concerns the specification of the null hypothesis in the CLEF stopping criterion. Originally, a certain conditional tail dependence coefficient, κ_α, related to a given group of variables α⊂{1, …, d} is supposed to be above a strictly positive, user-defined and therefore somewhat arbitrary threshold. We propose instead to base the stopping criterion on the hypothesis that a multivariate version of the coefficient of <cit.> and <cit.> is equal to one. The test is based on the limit distributions of multivariate extensions of nonparametric estimators in <cit.> and <cit.>. Third, we conduct a numerical experiment to compare the finite-sample performance of the DAMEX algorithm and the CLEF algorithm with the various stopping criteria. We find that overall, the multivariate extension of the Hill-type estimator in <cit.> yields the most reliable procedure to detect maximal groups of asymptotically dependent variables.

Section <ref> casts the problem in the language of regular variation and introduces the tail dependence coefficients upon which the CLEF stopping criteria will be based. Necessary background on empirical tail dependence functions and processes is reviewed in Section <ref>, including a new result for the empirical joint tail function. In Section <ref>, we derive the asymptotic distribution of the statistic used in CLEF and turn the heuristic stopping criterion implemented in <cit.> into a statistical test with asymptotically controllable level. Two alternative tests based on the asymptotic distributions of estimators of the Ledford–Tawn–Ramos coefficient of tail dependence are constructed in Sections <ref> and <ref>. We report the results of our simulation experiments in  Section <ref>. Section <ref> concludes. Proofs are gathered in Appendix <ref> while the pseudo-code for the CLEF algorithm and variations is provided in Appendix <ref>. 
























§ REGULAR VARIATION AND TAIL DEPENDENCE COEFFICIENTS



Bold letters denote vectors and binary operations between vectors are understood componentwise. 

The indicator function of a set A is denoted by _A. For t∈∪{∞}, we let t_α denote the constant vector of (∪{∞})^α with all coordinates equal to t. In the special case α = {1, …, d}, the index α is usually omitted for brevity when clear from the context: for instance, = _{1,…,d} = (0, …, 0) ∈^d.

Let X = (X_1, …, X_d) be a random vector in ^d with cumulative distribution function F, whose margins F_1, …, F_d are continuous. We assume that the transformed vector V = (V_1, …, V_d) with V_j = 1/{1 - F_j(X_j)} for all j ∈{1,…,d} is regularly varying on the cone [0,∞]^d ∖{} with (nonzero) limit or exponent measure μ. This means that μ is finite on Borel sets of [0,∞]^d ∖{} bounded away from the origin and that 

    lim_t →∞ t  [  V∈  t A] = μ(A),

for all Borel sets A ⊂ [0,∞]^d ∖{} such that ∉∂ A and μ(∂ A)=0. The measure μ is homogeneous, i.e., μ(s  · ) = s^-1μ( · ) for all 0 < s < ∞, and therefore assigns no mass to hyperplanes parallel to the coordinate axes. As a consequence, (<ref>) applies to finite and infinite rectangles that are bounded away from the origin and whose sides are parallel to the coordinate axes. The measure μ characterizes the extremal dependence structure of X. The reader is referred to <cit.> for an introduction to regular variation.  

Let ∅α⊂. Particular instances of (<ref>) include the extremal coefficient λ_α <cit.> and the joint tail coefficient 
ρ_α: 

    λ_α   =
      lim_t →∞  t  [ ∃ j ∈α : V_j >  t ]
      =
      μ ( {u∈ [0, ∞)^d |∃ j ∈α : u_j > 1 } ),
      
    ρ_α   =
      lim_t →∞  t  [ ∀ j ∈α : V_j >  t ]
      =
      μ ( {u∈ [0, ∞)^d |∀ j ∈α : u_j > 1 } ).


In the bivariate case, |α|=2, and with our choice of Pareto margins, we have ρ_α = lim_t→∞(V_α_1>t | V_α_2>t), the upper tail dependence coefficient denoted by χ in <cit.>.

Our general objective is to propose statistically sound procedures to recover maximal subgroups α of components that are likely to be concomitantly large. Our aim can thus be phrased as recovering the maximal subsets α⊂ such that ρ_α>0.

Since ρ_α≤ρ_β as soon as α⊃β, any positive tolerance level with which we would like to compare an estimate of ρ_α should depend on α and in particular be decreasing as a function of the cardinality |α|. To circumvent this issue, <cit.> consider for α such that |α| ≥ 2 the conditional tail dependence coefficient

    κ_α =
      lim_t→∞[
      ∀ j ∈α : V_j > t 
       | ∑_j∈α{V_j > t}≥ |α| -1    
      ],

which is the limiting conditional probability that all variables in α exceed a large threshold given that all but at most one already do. In contrast to ρ_α, the coefficient κ_α has no particular reason to decrease as a function of |α|. Note that ρ_α = μ(Γ_α) while κ_α = μ(Γ_α) / μ(Δ_α) = ρ_α / μ(Δ_α) where Γ_α = {x∈ [0, ∞)^d |∀ j ∈α : x_j > 1 } and Δ_α = { x ∈ [0, ∞)^d |∑_j ∈α_{x_j≥ 1}≥ |α|-1}, provided |α|≥ 2. If μ(Δ_α) = 0, then μ(Γ_β) = 0 for all β⊂α with |β| = |α| - 1; in that case, we define κ_α = 0.







In the CLEF algorithm <cit.>, the criterion to decide whether ρ_α>0 or not is that κ_α≥ C, where C is a user-defined tolerance level, κ_α = μ(Γ_α)/μ(Δ_α), and μ is the empirical exponent measure in (<ref>) below. The level C can be chosen independently of α. Still, its choice is somewhat arbitrary, and in particular, the user has no control of false positives. In Section <ref>, we will provide the asymptotic distribution of κ_α and propose a test statistic with a guaranteed asymptotic level.

If ρ_α=0 (or κ_α=0), the limiting distributions of the statistics √(k) (ρ_α - ρ_α) and √(k) (κ_α - κ_α) are degenerate at zero.  We therefore have no control on the asymptotic levels of tests based on those statistics under H_0 : κ_0 = 0. This is why will have to define a CLEF stopping criterion in terms of a test of H_0 : κ_α≥κ_min versus H_1 : κ_α < κ_min instead, in terms of a user-defined level κ_min > 0. The choice of κ_min is somewhat arbitrary; in the simulation experiments (Section <ref>), we choose κ_min = 0.08. 

In Sections <ref> and <ref>, we consider alternative CLEF stopping criteria based on estimators of the coefficient of tail dependence η_α∈ (0,1]. For bivariate distributions, the coefficient has been introduced by <cit.> and extended by <cit.> in order to model situations in between asymptotic dependence (ρ_{1,2} > 0) and full independence of X_1 and X_2. <cit.> and <cit.> proposed and studied a multivariate extension of η_α for |α| ≥ 3. The model assumption is that there exist η_α∈(0,1] and a slowly varying function ℒ_α such that

    [ ∀ j ∈α : V_j > t ] = t^-1/η_αℒ_α(t).



Suppose that the limit ρ_α in (<ref>) exists and that (<ref>) holds. Then ρ_α > 0 implies η_α = 1. The converse is true as well, provided lim inf_t→∞ℒ_α(t) > 0. Modulo this side condition, which we will take for granted, the null hypothesis ρ_α>0 corresponds to the simple hypothesis η_α = 1. 

We will test the null hypothesis η_α = 1 via multivariate extensions of nonparametric estimators of η_α in <cit.> and <cit.>. The null limit of the test statistic is non-degenerate, so that the asymptotic level of the test can be controlled, with no need to introduce an additional tolerance parameter κ_min. The estimators that we will study are related to the Pickands estimator and the Hill estimator for the extreme value index of T_α = min_j∈α V_j, respectively. The maximum likelihood estimator, also considered in <cit.>, is less suitable to our context due to its relative computational complexity, since the test is destined to be performed on a large number of subsets of {1,…,d}. See also the review <cit.> and the references therein.



The DAMEX algorithm <cit.> is designed to recover the family ℳ of non-empty subsets α of  with the property that

    μ(
        { x ∈ [0, ∞)^d  |  x_∞≥ 1 ; ∀ j ∈α , x_j > 0  and ∀ j
        ∉α,  x_j = 0 })
      > 0.

In contrast, our focus is on 𝕄 = {α|ρ_α > 0} = {α|κ_α > 0 }. Still,  the maximal elements of 𝕄 for the inclusion order are also the maximal elements of ℳ <cit.>. The two problems of finding the maximal elements of 𝕄 or ℳ are thus equivalent.






§ EMPIRICAL TAIL DEPENDENCE FUNCTIONS AND PROCESSES



To find the asymptotic distribution of nonparametric estimators of the
various dependence coefficients, we rely on empirical tail
processes. Let the random vector X∼ F be as in
Section <ref>; in particular, assume regular
variation as in (<ref>) with exponent measure μ. Let
Λ be the push-forward measure of μ on
[0, ∞]^d ∖{∞} induced by the
transformation x↦ 1/x = (1/x_1, …, 1/x_d),
i.e.,
Λ( · ) = μ( {x∈ [0, ∞]^d ∖{0}| 1/x∈ · } ).

For ∅α⊂, consider the stable tail dependence function ℓ_α : [0, ∞)^α→ [0, ∞) 
and the joint tail dependence function r_α : [0, ∞]^α∖{∞_α}→ [0, ∞) given by

    ℓ_α( x) 
         = lim_t → 0 t^-1 [ ∃ j ∈α : F_j(X_j) > 1 -  t x_j]
      = Λ( {y|∃ j ∈α : y_j < x_j } ),
      
    
      r_α( x) 
         = lim_t → 0 t^-1 [ ∀ j ∈α : F_j(X_j) > 1 - t x_j] 
      = Λ( {y|∀ j ∈α : y_j < x_j } ).

From (<ref>) and (<ref>), clearly λ_α = ℓ_α(1_α) and ρ_α = r_α(1_α). For brevity, we write ℓ = ℓ_ and r = r_. Note that ℓ_α( x ) = ℓ( xe_α ) for x∈ [0, ∞)^α, where e_α∈{0, 1}^d has components e_α, j = _α(j). Similarly, r_α( x ) = r( xι_α ) for x∈ [0, ∞]^α∖{∞_α}, where _α∈{1,∞}^d denotes the vector such that _α,j = 1 if j ∈α and _α,j = +∞ otherwise. By the inclusion–exclusion formula, for x∈ [0, ∞)^α, 
writing x_β = (x_j)_j ∈β, we have

    r_α( x) 
         = 
      ∑_∅β⊂α
      (-1)^|β|+1ℓ_β ( x_β),    ℓ_α( x) 
         = 
      ∑_∅≠β⊂α
      (-1)^|β|+1 r_β ( x_β).

Let X_i = (X_i,1, …, X_i,d), for i ∈{1,…,n}, be an independent random sample from F, having continuous margins and satisfying (<ref>). Let k = k(n)→∞ as n→∞, while k(n) = (n). Following for instance <cit.>, we rely on ranks to obtain an approximately Pareto-distributed sample V_i = (V_i,1, …, V_i,d). Let F_j(x) = n^-1∑_i=1^n _{X_i,j < x} be the (left-continuous) empirical distribution function of component j ∈{1,…,d} and put V_i,j = 1/{1 - F_j(X_i,j)} = n/(n + 1 - R_i,j), where R_i,j 
is the rank of X_i,j among X_1,j, …, X_n,j. The empirical counterparts to μ and Λ are

    μ( · ) 
           = 
        1/k∑_i=1^n δ_ (k/n) V_i( · ),    Λ( · )
           =
        1/k∑_i=1^n δ_ (n/k)/ V_i( · ),

respectively, with δ_a the Dirac measure at the point a. Replacing Λ by Λ in the definition of ℓ_α and r_α produces the empirical tail dependence function

    ℓ_α( x )
         = k^-1∑_i=1^n {∃ j ∈α : n + 1 - R_i,j≤⌊ k x_j ⌋}
       = k^-1∑_i=1^n {∃ j ∈α : X_i,j≥ X_(n-⌊ k x_j ⌋ +1),j}

and the empirical joint tail function

    r_α( x )
         = k^-1∑_i=1^n {∀ j ∈α : n + 1 - R_i,j≤⌊ k x_j ⌋}
       = k^-1∑_i=1^n {∀ j ∈α : X_i,j≥ X_(n-⌊ k x_j ⌋ +1),j},

where X_(1),j≤…≤ X_(n),j are the ascending order statistics of X_1,j, …, X_n,j and ⌊ · ⌋ is the floor function. The identities (<ref>) hold for ℓ_α and r_α as well.

<cit.> find the weak limit of the empirical process √(k) ( ℓ - ℓ ) on [0, T]^d for any T > 0. We leverage their theorem to show a similar result for √(k)( r_α - r_α ), jointly in α. The following conditions stem from the cited article.
[Uniform tail convergence]

There exists γ > 0 such that, uniformly in x ∈ [0,1]^d with ∑_j=1^d x_j = 1, we have 

    t^-1[ ∃ j = 1, …, d :   F_j(X_j) > tx_j ] - ℓ( x) 
      = (t^γ),      t →∞.



[Moderate k]

The sequence k = k(n) satisfies k = (n^2γ/(1+2γ)) as n →∞, with γ > 0 as in Condition <ref>.


[Smoothness]

For all j ∈{1,…,d}, the partial derivative ∂_j ℓ = ∂ℓ / ∂ x_j exists and is continuous on the set { x ∈ [0,∞)^d | x_j > 0}.


Since ℓ is convex, it is continuously differentiable Lebesgue almost everywhere <cit.>. Condition <ref> is satisfied for many popular max-stable models (logistic, asymmetric logistic, Brown–Resnick) but fails for max-linear models. Under Condition <ref>, the partial derivative ∂_j r_α = ∂ r_α / ∂ x_j (j ∈α) exists and is continuous on {x∈ [0, ∞)^α| x_j > 0 } and satisfies ∂_j r_α( x ) = ∑_β : j ∈β⊂α (-1)^|β|+1∂_j ℓ_β(x_β), where x_β = (x_s)_s ∈β.

<cit.> and <cit.> consider a centered Gaussian process W indexed by the Borel sets of [0,∞]^d∖{∞} bounded away from ∞ with covariance function 

    [W(A)   W(B)] = Λ(A ∩ B ).

Note that W(∅) = 0 almost surely. For ∅≠α⊂{1,…, d} and x∈ [0, ∞)^α, write

    W_α( x ) = W({y∈ [0, ∞]^d |∀ j ∈α : y_j < x_j }).

We consider weak convergence as in <cit.>; notation . We work in the metric space ℓ^∞(S) of bounded, real functions f on an arbitrary set S, the metric being the one induced by the supremum norm, f _∞ = sup_x ∈ S| f(x) |; the double use of the symbol ℓ should not give rise to any confusion. The proof of the following proposition and of other results in the paper is deferred to Appendix <ref>.



Let X_i = (X_i,1, …, X_i,d), for i ∈{1,…,n}, be an independent random sample from F, having continuous margins and satisfying (<ref>). Let k = k(n)→∞ as n→∞, while k(n) = (n). If Conditions <ref>, <ref> and <ref> hold, then, for T > 0, in the product space ∏_∅α⊂ℓ^∞( [0, T]^α ), we have, as n →∞, the weak convergence

    √(k){r_α( x) - r_α( x)}
      W_α( x) - ∑_j∈α∂_j r_α( x)   W_{j}(x_j)
      =
      Z_α( x ).

  



























§ ESTIMATING THE CONDITIONAL TAIL DEPENDENCE COEFFICIENT



This section investigates the asymptotic distribution of the empirical conditional dependence coefficient κ_α based on the empirical exponent measure μ. This is achieved by re-writing  κ_α  as a function of the empirical joint tail coefficients  ρ_α, the distribution of which follows from Proposition <ref>. We also propose consistent estimators of the asymptotic variance of κ_α. Combining the two yields a test for the null hypothesis κ_α≥κ_min where κ_min∈(0,1) is a tolerance level fixed by the user, to be seen as the minimal limiting conditional probability that all components in a random vector exceed a threshold, given that all of them but at most one already do.  

Let ∅α⊂ and recall the sets Γ_α = {x∈ [0, ∞)^d |∀ j ∈α : x_j > 1 } and, provided α has at least two elements, Δ_α = { x ∈ [0, ∞)^d |∑_j ∈α_{x_j≥ 1}≥ |α|-1}. Write α∖ j = α∖{j} for j ∈α. Since Δ_α is the disjoint union of the sets Γ_α∖ j∖Γ_α and Γ_α, where j ∈α, we find, for every Borel measure ν, the equality

    ν( Δ_α ) 
      = ∑_j ∈αν(Γ_α∖ j) - (|α| -1)  ν(Γ_α).

Recall ρ_α = μ(Γ_α) and κ_α = μ(Γ_α) / μ(Δ_α) in (<ref>). By (<ref>) applied to ν = μ, we have

    κ_α 
      = 
      ρ_α/∑_j ∈αρ_α∖ j - (|α| - 1)ρ_α.


Recall the joint tail function r_α and its nonparametric estimator r_α in (<ref>) and (<ref>), respectively. Since ρ_α = r_α( 1_α ), we define the estimators ρ_α = μ( Γ_α ) = r_α( 1_α ) and, provided |α|≥ 2,

    κ_α
      =
      μ( Γ_α )/μ( Δ_α )
      =
      ρ_α/∑_j ∈αρ_α∖ j - (|α| - 1) ρ_α.

The asymptotic distribution of the vector of empirical joint tail coefficients follows immediately from Proposition <ref>. Write ρ̇_α, j = ∂_j r_α(1_α). 


 
In the setting of Proposition <ref>, we have, jointly in ∅α⊂, the weak convergence

    √(k_n)( ρ_α - ρ_α)
       Z_α( 1_α ) = G_α,
           n →∞.

The limit distribution is centered Gaussian with covariance matrix

    [ G_α G_α' ]
      = ρ_α∪α' -
      ∑_j ∈αρ̇_j,αρ_α' ∪{j} - ∑_j' ∈α'ρ̇_j',α'ρ_α∪{j'}
      +∑_j ∈α∑_j' ∈α'ρ̇_j,α ρ̇_j',α' ρ_{j, j'}.



The asymptotic distribution of κ_α follows from the one of (ρ_β)_β via the delta method. The asymptotic variance involves the partial derivative ∂_j κ_α = ∂κ_α / ∂ x_j of the function

    κ_α( x) 
       = r_α( x)/∑_j ∈α r_α∖ j( x_α∖ j) - (|α| -1) r_α( x)

for x ∈ [0,∞)^α. Note that κ_α( 1_α) = κ_α. Write κ̇_j,α = ∂_j κ_α( 1_α ).
  


In the setting of Corollary <ref>, we have, as n →∞ and jointly in α⊂{1,…,d} such that |α| ≥ 2 and μ(Δ_α) > 0, the weak convergence
  
    √(k)( κ_α - κ_α)
        μ(Δ_α)^-2{( ∑_j∈αρ_α∖ j) G_α
          - ρ_α∑_j∈α G_α∖ j}.

For a fixed such α, the limit distribution is 𝒩(0, σ_κ,α^2) with
  
    σ^2_κ, α = (1 - κ_α)κ_α{μ(Δ_α)^-1 - ∑_j∈ακ̇_j,α}
        + ∑_i∈α∑_j∈ακ̇_i,ακ̇_j,αρ_{i,j}
    
        + κ_α∑_ j ∈ακ̇_j,α{ 1 - μ(Δ_α)^-1ρ_α∖ j}.



Following <cit.>, the asymptotic variance σ^2_κ,α in (<ref>) can be estimated consistently by estimating the partial derivatives κ̇_i,α via finite differencing applied to the empirical version of κ_α(x) in (<ref>) obtained by replacing r_α and r_α∖ j by r_α and r_α∖ j, respectively:

    κ_α( x) 
      = 
      ∑_i=1^n {∀ j ∈α : X_i,j ≥ X_(n - ⌊ kx_j ⌋ +1 ), j}/∑_i=1^n {∃ m∈α: ∀ j ∈α∖ m : X_i,j ≥  X_(n - ⌊ kx_j ⌋ +1 ), j}

Define

    κ̇_j, α, n = 1/2 k^-1/4{κ_α(
      1_α + k^-1/4 e_j) - κ_α( 1_α - k^-1/4 e_j ) },

with e_j the canonical unit vector of ^α pointing in direction j ∈α, and put

    σ^2_κ, α = (1 -
      κ_α)κ_α{μ(Δ_α)^-1 - ∑_j∈ακ̇_j,α, n}
      + 
      ∑_i,j∈ακ̇_i,α, nκ̇_j,α,nρ_{i,j}
    
      + 
      κ_α∑_ j ∈ακ̇_j,α,n{ 1 - μ(Δ_α)^-1ρ_α∖ j}.


  

Under the conditions of Proposition <ref>, we have
  σ^2_κ,α = σ^2_κ,α + _(1) as n →∞,
  so that √(k)(κ_α - κ_α) / σ_κ,α𝒩(0,1), provided σ^2_κ,α > 0.


The proof relies on the weak convergence of the empirical process √(k){κ_α( · ) - κ_α( · ) } on [0,T]^α for any T>0. This property follows in turn from Proposition <ref> and the functional delta method.

We consider a tolerance level κ_min∈ (0, 1) under which the tail dependence between components j∈α is deemed negligible compared to the one between components j∈β⊊α.  In other words, we aim at testing H_0: κ_α≥κ_min. Since κ_α = ρ_α / μ(Δ_α), the null hypothesis is that ρ_α is greater than some level depending on α. Let 0<δ<1 be a (small) probability, and consider the test

    τ_α,n 
      = 
      {κ_α < κ_min + q_δ k^-1/2σ_κ,α}

where q_δ is the δ-quantile of the standard normal distribution. By Proposition <ref>, if σ_κ,α > 0, the test in (<ref>) has asymptotic level δ for H_0 against H_1: κ_α< κ_min.

If ρ_α = 0, then, in Proposition <ref>, we have √(k)( ρ_α - ρ_α ) = _(1) as n →∞: indeed, on the one hand, we have √(k)( ρ_α - ρ_α ) = √(k)ρ_α≥ 0, and on the other hand, its limit distribution is centered Gaussian. Likewise, we have √(k)( κ_α - κ_α ) = _(1) as n →∞ in Proposition <ref> if κ_α = 0. As a consequence, under the simple hypothesis H_0 : ρ_α = 0, the asymptotic level of a test based on the asymptotic distribution of √(k)(ρ_α - ρ_α) or √(k)(κ_α - κ_α) cannot be controlled. This is why the test in (<ref>) concerns the null hypothesis H_0 : κ_α≥κ_min for some κ_min > 0 instead. Alternatively, we propose tests based on estimators of the coefficient of tail dependence η_α in (<ref>). In Sections <ref> and <ref>, we consider two such estimators, extending the ones of <cit.> and <cit.>, respectively, to the multivariate setting.







§ COEFFICIENT OF TAIL DEPENDENCE: PENG'S ESTIMATOR



For bivariate distributions, Peng's <cit.> estimator of the coefficient of tail dependence η = η_{1,2} is based on the property that the curve t↦ (log t, log[V_1>t, V_2>t]) has an affine asymptote with slope - 1/η. A similar idea motivates Pickands' <cit.> estimator for the extreme value index. Estimating the ordinate of the curve at t = n/k and t = n/(2k) allows to estimate that slope. Under a second-order regular variation condition, <cit.> shows that his estimator is asymptotically normal, both if η = 1 and if η < 1. In the former case, the asymptotic variance depends on the tail dependence function and its partial derivatives, which are unknown but may be estimated consistently, thus leading to tests whose asymptotic levels can be controlled.

Let α⊂ have at least two elements. Recall the empirical joint tail function r_α in (<ref>). We define the multivariate extension of Peng's <cit.> estimator of η_α in (<ref>) as 

    η_α^P 
      = log(2) /
      log{r_α( 2_α) / r_α( 1_α) }.

The asymptotic normality of η_α^P follows from Proposition <ref> and
the delta method.



In the setting of Proposition <ref>, we have, as n →∞ and jointly in α⊂ such that |α|≥ 2 and ρ_α > 0, the weak convergence
  
    √(k) (η_α^P - 1) -1/2ρ_αlog 2{ Z_α( 2_α) - 2 Z_α( 1_α) }.

  The right-hand side is a 𝒩(0, σ_α,P^2) random variable with variance
  
    σ_α,P^2 = 1/2(ρ_αlog 2)^2[
          ρ_α - 4 ρ_α^2 +  2 ∑_j∈αρ̇_j,α 
          r_α( 2_α∧_j) 
    
          + ∑_j∈α∑_j'∈αρ̇_j,αρ̇_j',α{ 3 ρ_{j, j'} - 2
          r_{j, j' }(2,1) }],

  where ρ_{j, j'} = r_{j, j'}(2, 1) = 1 if j = j' and
  where _j ∈{1,∞}^α is the vector which all coordinates equal to 1 except for the j-th one which equals ∞, so that     ( 2_α∧_j)_m = 1 if  m ∈α∖ j and ( 2_α∧_j)_m = 2 if m = j.



By extending the proof of <cit.>, it is also possible to obtain asymptotic normality of η_α^P in the case ρ_α=0 and η_α < 1 in (<ref>). This would require a   multivariate extension of the second-order regular variation condition in <cit.> in the style of Condition <ref> below. For the application as a stopping criterion in the CLEF algorithm, we are only interested in the asymptotic distribution of η_α^P under the hypothesis ρ_α>0, so we do not pursue this idea any further. 



As in Proposition <ref>, the asymptotic variance σ_α,P^2 in (<ref>) involves unknown quantities, all of which we can estimate consistently. For α⊂{1,…,d} and j∈α, define

    ρ̇_j,α,n 
      = 
      1/2 k^-1/4{r_α( 1_α + k^-1/4 e_j) 
        - 
        r_α( 1_α - k^-1/4 e_j) 
      },

where e_j is the canonical unit vector in ^α pointing in dimension j. Define

    σ^2_α,P 
      = 1/2(ρ_αlog 2)^2[ 
        ρ_α 
        + 
        ∑_j∈αρ̇_j,α, n{ 
    	- 4 ρ_α 
    	+ 2 r_α( 2_α∧_j) 
          }
    
        + 
        ∑_j∈α∑_j'∈αρ̇_j,α,nρ̇_j',α,n{ 3 ρ_{j, j'} - 2 r_{j, j'}(2,1) }].




In the setting of Proposition <ref>, we have σ_α, P^2 = σ_α, P^2 + _(1) as n →∞, where α⊂ is such that |α|≥ 2 and ρ_α > 0. If σ_α, P^2 > 0, then √(k) (η_α^P - 1) / σ_α, P𝒩(0,1) as n →∞.


The proof parallels the one of Proposition <ref> and is omitted for brevity. The main step is to verify that ρ̇_j,α, n = ρ̇_j,α + _(1) as n →∞, which follows from Proposition <ref>. 

To test the hypothesis H_0 : ρ_α > 0 at significance level δ∈ (0, 1), we propose

    τ_α,η^P,n 
      = 
      {η_α^P 
        < 1 - q_1-δ k^-1/2σ_α,P},

where q_1-δ is the (1-δ)-quantile of the standard normal distribution. In the setting of Proposition <ref>, the test in (<ref>) has asymptotic level δ for H_0 against H_1: η_α < 1.





§ COEFFICIENT OF TAIL DEPENDENCE: HILL ESTIMATOR




The coefficient of tail dependence η_α in (<ref>) is the tail index of the random variable T_α = min_j ∈α V_j: the function t ↦[T_α > t] is regularly varying at infinity with index -1/η_α. A tractable alternative to Peng's estimator for η_α is a Hill-type estimator as in <cit.>. Replacing the unobservable Pareto variables V_i,j by the rank-based versions V_i,j = n / (n + 1 - R_ij) in Section <ref> yields an approximate sample

    T_i, α = min_j∈αV_i,j,
           i = 1, …, n,

from the distribution of T_α. Let T_(1), α≤…≤T_(n), α denote the order statistics of T_1,α, …, T_n,α. The Hill estimator for η_α is defined as

    η_α^H 
      = 1/k∑_i=1^k 
      logT_(n-i+1), α/T_(n-k), α.

Under the second-order regular variation conditions stated below, the asymptotic normality of η_α^H follows from <cit.>. The results in the cited reference cover the bivariate case only. In this section, we verify that they remain valid in any dimension d≥ 2, and we provide the general expression for  the asymptotic variance. Put E_α = [0, ∞]^α∖{∞_α}.






For each α⊂{1,…,d} with |α|≥ 2, there exist functions c_α, c_1,α: E_α→ [0, ∞)   such that c_1,α is neither constant nor a multiple of c_α, and there exists q_1,α: (0, ∞) → (0, ∞), with q_1,α(t) → 0 as t→ 0, such that, for all x ∈ E_α, we have

    lim_t → 0{[ ∀ j ∈α : 1 - F_j(X_j) ≤ tx_j]/[ ∀ j ∈α : 1 - F_j(X_j) ≤ t ] 
          - c_α( x) 
        }/ 
        q_1, α(t) 
        = c_1,α( x).




Under Condition <ref>, the function q_α(t) = [ ∀ j ∈α : 1 - F_j(X_j) ≤ t ] is regularly varying at 0 with some index 1/η_α. Condition <ref> implies that the first-order condition (<ref>) holds with the same index 1/η_α.  In addition, c_α( 1_α) = 1 and c_α is homogeneous of order 1/η_α, i.e., c_α(t  x ) =
t^1/η_αc_α( x) for t>0, see <cit.>. Under the regular variation assumption (<ref>), we have  ρ_α = lim_t → 0 q_α(t)/t, so that, under Condition <ref>, ρ_α > 0 implies η_α = 1, as in <cit.> for the bivariate case. Finally, if ρ_α>0, then c_α( x )= r_α( x) /r_α( 1_α) = r_α( x) / ρ_α. Note that in <cit.>, our ρ_α is denoted by l for α = {1,2}. 

The asymptotic variance of the Hill estimator (<ref>) involves a Gaussian process whose distribution depends on whether ρ_α=0 or ρ_α > 0. As in <cit.>, introduce a centered Gaussian process W_1 on E_α with covariance function [ W_1( x)   W_1( y) ]  = c_α(  x∧ y) for x, y ∈ E_α. Recall the stochastic process Z_α in (<ref>) and the random variable G_α = Z_α(1_α) in (<ref>). 



Let X_i = (X_i,1, …, X_i,d), for i ∈{1,…,n}, be an independent random sample from F, having continuous margins and satisfying (<ref>). Let k = k(n)→∞ as n→∞, while k(n) = (n). If Conditions <ref>, <ref>, <ref>, and <ref> hold, then, as n →∞,
  
    √(k)( η_α^H - η_α)
        𝒩(0,σ_α,H^2),

with σ_α,H^2 = η_α^2 {W̃( 1_α) }, where W̃( x) = W_1( x) if ρ_α = 0 and W̃(x) = ρ_α^-1/2 Z_α( x) if ρ_α > 0. In particular, if ρ_α > 0, we have

    σ_α,H^2
      =
      ρ_α^-1(G_α)
        = 1 - 2 ρ_α
          
        +
        ρ_α^-1∑_j ∈α∑_j' ∈αρ̇_j,αρ̇_j',αρ_{j,j'}.



The proof of Proposition <ref> is based on the arguments developed in the proofs of <cit.>, <cit.>, and <cit.>, which we gather in Appendix <ref>.   
























Again, the unknown terms in (<ref>) may be replaced by their empirical counterparts, leading to an asymptotically consistent test. Recall ρ̇_j,α,n in (<ref>) and define
  
    σ^2_α, H
        = 1 - 2 ρ_α
          
        +
        ρ_α^-1∑_j∈α∑_j'∈αρ̇_j,α,nρ̇_j',α, nρ_{j,j'}.

  
The proof of the consistency of the variance estimator follows the same lines as the proofs of Propositions <ref> and <ref> and is omitted.



Under the conditions of Proposition <ref>, if ρ_α > 0, we have σ_α,H^2 = σ_α,H^2 + _(1) as n →∞ and thus √(k) (η_α^H - 1) / σ_α, P𝒩(0,1), provided σ_α,H^2 > 0.


We may exploit Corollary <ref> to test H_0 : ρ_α > 0 in the same way as we did by using Peng's estimator in (<ref>): at significance level δ∈ (0, 1), the null hypothesis is rejected in favour of H_1: η_α < 1 when η_α^H < 1 - q_1-δ k^-1/2σ_α,H.


The condition σ_α,H^2 > 0 in Corollary <ref> is satisfied whenever 0<ρ_α<1. Indeed, in (<ref>), we have ρ_{j,j'}≥ρ_α and ρ̇_j,αρ̇_j',α≥ 0, whence σ_α,H^2 ≥ 1 - 2ρ_α + ∑_(j,j')∈α^2ρ̇_j,αρ̇_j',α = 1 - 2 ρ_α + ρ_α^2 = (1 - ρ_α)^2. 





§ SIMULATION STUDY





Our aim is to compare the finite sample performance of the various tests proposed in Sections <ref>, <ref> and <ref> within the framework of the CLEF algorithm, the pseudo-code of which is given in Appendix <ref>. Three variants of the CLEF algorithm are obtained by varying the criterion according to which a subset α is declared as tail-dependent:
 
 
 
  κ_α > κ_min - q_δσ_κ,α / √(k) for CLEF-asymptotic;



 η_α,P > 1 - q_δσ_α,P / √(k) for CLEF-Peng; and



 η_α,H > 1 - q_δσ_α,H / √(k) for CLEF-Hill.

The original CLEF criterion was κ_α > C for some constant C chosen by the user.
For completeness, the output of the DAMEX algorithm <cit.> is included in the comparison. 
 
 













 
 
 
 
 
 

In practice, the dependence tests based on the tail dependence coefficient should not be carried out to the letter when the test statistic is not defined or when its estimated variance is infinite. Thus, in our experiments, CLEF-Peng and CLEF-Hill are modified so as to take into account additional, common-sense stopping criteria. A subset α will not be part of the list returned by the algorithms under the following conditions: 

 

 
 
 
 
 


 
 
  * Concerning  CLEF-Hill, when ρ_α = 0, that is, no extreme record  impacts  all coordinates in α, the estimated variance of the Hill estimator of η_α is infinite. Therefore, ρ_α = 0 is considered as a stopping criterion in  CLEF-Hill.
 
  * Concerning CLEF-Peng, when r_α( 2_α) = r_α( 1_α), the Peng estimator (<ref>) is ill-defined. Such a case arises when there are very few points in the joint tail within the subspace generated by α.


When the estimated derivatives ρ̇_j,α,n  are close to zero, and when ρ_α≪ 1,  the estimated variance σ_α,P^2 in (<ref>) becomes large, preventing rejection of the null hypothesis. 




To prevent these issues, each of the conditions ρ_α< 0.05  and r_α( 2_α) = r_α( 1_α) are declared as a stopping criterion in CLEF-Peng.














 
 


 




  
Experimental setting.
CLEF <cit.> is designed to face situations where  DAMEX <cit.> fails to exhibit a clear-cut dependence  structure. A major issue reported in <cit.> for certain hydrological data is the high variability of the groups of features for which large values occur simultaneously. Because of this, the empirical exponent measure μ assigns low mass to any sub-region partitioning the sample space, see Remark <ref>. The empirical finding motivating the latter work is that the various subsets α involved in simultaneous extreme records could nevertheless be clustered, meaning that many of them have a significant intersection, whereas many symmetric differences comprise just a single or at most a few features.

A natural assumption in this context is that a `true' list of dependent subsets ℳ = {α_1,…, α_K} exists such that μ(𝒞_α)>0 for α∈ℳ and that noisy features are involved in each extreme event. Observed large records then concern groups of the kind α' = α∪{j}, where α∈ℳ and j∈∖α.
  
In our experiments, datasets are generated as follows: The  dimension is fixed to d=100. A family of `true' dependent subsets ℳ= {α_1,…,α_K} of cardinality K = 80 is randomly chosen: the subset sizes |α| follow a truncated geometric distribution, with a maximum subset size set to 8. For simplicity, we forbid nested subsets, so α_j⊄α_k whenever j  k. The maximal elements of 𝕄 = {α⊂|ρ_α>0} are then precisely the elements of ℳ, as explained in  Remark <ref>. Finally, two different subsets may have at most two features in common. 

Once the dependence structure ℳ has been fixed, the data X_1, …, X_n are sampled independently from  d-dimensional asymmetric logistic distributions <cit.>, using Algorithm 2.2 in
<cit.>. The underlying `true' distribution function is 

    G( x) = exp[- ∑_m=1^K{∑_j∈α_m (| A(j)| x_j)^-1/w_α_m}^w_α_m],

where A(j) = {α∈ℳ| j ∈α} and w_α_m is a dependence parameter which is set to 0.1 in our simulations. Actually, to mimic the noisy situation described above, each point X_i is simulated according to a slightly different version, G_i, of G. For each i = 1, …, n and k = 1, …, K, we randomly select an additional `noisy feature' j_i,k∈{1, …, d}∖α_k and set α_i,k' = α_k ∪{ j_i,k}. Then ℳ_i' = {α_i,1', …, α_i,K' } is the collection of `noisy subsets' for X_i and G_i(x) is as in (<ref>) with 𝒜(j) replaced by 𝒜'_i(j) = {α' ∈ℳ'_i | j ∈α'}.


























































  
Results.
We generate datasets of size n=5e4 and n=1e5. For each sample size, 50 independent datasets are simulated according to the
procedure summarized in the preceding paragraph. We compare the average performance of the three proposed versions of CLEF, together with the original CLEF and DAMEX algorithms,  for different choices of k
 

and confidence level δ.  






































Tables <ref>
and <ref> gather the results 

for a  confidence level δ equal to 0.001 and 0.0001, respectively.
In both tables, 

the results obtained with the original version of CLEF
and DAMEX are included in the comparison with an identical choice of tuning parameters, so that the last two lines of the two tables are the same. 




In CLEF,  the threshold C was chosen  by trial and error  in the interval (0, κ_min), namely C=0.05. 


Imposing that C < κ_min is intended to reproduce the effect of the variance term upon the stopping criterion  in CLEF-asymptotic.   
In DAMEX, the 80 subsets with highest empirical mass are retained and  the subspace thickening parameter ϵ is set to the default  value of 0.1, following the guidelines of the authors. 




Each algorithm produces a list, 𝕄, of groups of features α∈{1, …, d}. This list is to be compared with the one of K = 80 `true' subsets ℳ. The performance of each algorithm is measured in terms of two criteria: the number of `true' subsets α∈ℳ that appear in 𝕄 (third column of Tables <ref> and <ref>); the number of `errors', that is, the subsets α∈𝕄 that do not belong to ℳ. These can be understood as `false positives'. Among these errors, we make the distinction between those which are  respectively proper subsets (fourth column of Tables <ref> and <ref>) or proper supersets (fifth column) of some true β∈ℳ, and the other errors (sixth column).

CLEF-Peng obtains  the best overall scores for both values of δ, but as explained above, a special treatment is reserved for the case ρ_α≤ 0.05, and this threshold constitutes an arbitrary tuning parameter, which can impact the performance significantly. On the other hand, CLEF-Hill does not require any other adjustment than for the special case ρ_α=0 and performs nearly as well as CLEF-Peng with δ=0.0001 and k/n=0.005. In addition, CLEF-Hill outperforms all the other methods. In particular, CLEF-asymptotic is globally less accurate than CLEF-Peng and CLEF-Hill. This  reflects the fact that the  null hypothesis in this algorithm  involves an arbitrary κ_min>0 fixed by the user. Our own choice κ_min = 0.08 was fixed by trial and error, which is straightforward with synthetic data and could also be achieved by cross-validation in a real use case.  Finally, as expected, DAMEX obtains very low scores, because it is not designed to handle the addition of noisy features, as explained earlier.



















































































































































































§ CONCLUSION



In this work, we propose three variants of the CLEF algorithm <cit.>, replacing the heuristic criterion in the original version with a formal test for asymptotic dependence, and this for all possible subsets of features among {1,…,d}. As in the original CLEF implementation, only a small proportion of all 2^d-1 subsets has to be examined, while the computational complexity for each such subset is low. Experimental results indicate that the CLEF algorithm is most effective when based on a test constructed from an extension of the Hill estimator <cit.> of the multivariate coefficient of tail dependence.

The procedure we propose is nonparametric and rank-based. Parametric approaches, based for instance on the nested asymmetric logistic distribution <cit.>, could have a greater sensitivity, at the cost of increased model risk and greater computational complexity. We have also assumed that the observations are serially independent; in the contrary case, the asymptotic variances of the various estimator need to be estimated by some form of bootstrap, which, in high dimensions, poses important theoretical and computational challenges; see <cit.> for the bivariate and serially independent case.







§ PROOFS




For ∅α⊂ and x∈ [0, ∞)^α, put 

    L_α(x) 
         = 
      {y∈ [0, ∞]^d |∃ j ∈α : y_j < x_j }, 
    
      R_α(x) 
         = 
      {y∈ [0, ∞]^d |∀ j ∈α : y_j < x_j }.

If α =, then just write L rather than L_. Note that L_α(x) = L(xe_α) with e_α = (_α(j))_j=1^d and that L_{j}(x_j) = R_{j}(x_j) and thus W(L_{j}(x_j)) = W_{j}(x_j). <cit.> show that, in the space ℓ^∞([0, T]^d) and under Conditions <ref>, <ref> and <ref>, we have weak convergence 

    √(k){ℓ(x) - ℓ(x) } 
      W(L(x)) - ∑_j=1^d ∂ℓ_j(x) W_{j}(x_j)

as n →∞. Here, we have taken a version of the Gaussian process W such that the trajectories x↦ W(L(x)) are continuous almost surely.

As in (<ref>), we have, for ∅α⊂ and x∈ [0, ∞)^α, the identity

    r_α(x) 
      = ∑_∅β⊂α (-1)^|β|+1ℓ( x_βe_β )

where x_β = (x_j)_j ∈β. Hence, we can view the vector (√(k)( r_α - r_α ))_∅α⊂ as the result of the application to √(k)(ℓ-ℓ) of a bounded linear map from the space ℓ^∞([0, T]^d) to the product space ∏_∅α∈ℓ^∞([0, T]^α). By the continuous mapping theorem, we obtain, in the latter space, the weak convergence

    √(k){r_α( x) - r_α( x)}∑_∅β⊂α 
      (-1)^|β|+1{
      W(L_β(x_β))
      -
      ∑_j=1^d ∂_j ℓ_β(x_β) W_{j}(x_j _β(j))
      }.

Here we used ℓ(x_βe_β) = ℓ_β(x_β).

The set-indexed process W satisfies the remarkable property that W(A ∪ B) = W(A) + W(B) almost surely whenever A and B are disjoint Borel sets of [0, ∞]^d ∖{∞} that are bounded away from ∞: indeed, (<ref>) implies [{W(A ∪ B) - W(A) - W(B)}^2] = 0. It follows that the trajectories of W obey the inclusion-exclusion formula, so that, for ∅α⊂ and x∈ [0, ∞)^α, we have, almost surely,

    ∑_∅β⊂α 
      (-1)^|β|+1 
      W(L_β(x_β))
         =
      ∑_∅β⊂α 
      (-1)^|β|+1 
      W( ⋃_j∈β R_{j}(x_j) ) 
       =
      W( ⋂_j∈α R_{j}(x_j) )
      =
      W(R_α(x))
      =
      W_α(x).

We can make this hold true almost surely jointly for all such α and x: first, consider points x with rational coordinates only and then consider a version of W by extending W_α to points x with general coordinates via continuity. Similarly, since W_{j}(0) = W(∅) = 0 almost surely, we have

    ∑_∅β⊂α (-1)^|β|+1∑_j=1^d ∂_j ℓ_β(x_β) W_{j}(x_j _β(j))
         =
      ∑_j ∈α∑_β : j ∈β⊂α∂_j ℓ_β(x_β) W_{j}(x_j) 
       =
      ∑_j ∈α∂ r_j(x) W_{j}(x_j).

We have thus shown weak convergence as stated in (<ref>).



  The weak convergence statement (<ref>) is a special case of (<ref>): set x = 1_α. The covariance formula (<ref>) follows from the fact that
  

    [W_α(1_α) W_α'(1_α')]
         =
      Λ( {y∈ [0, ∞]^d |∀ i ∈α∪α' : y_i < 1 } ) 
       =
      μ( {u∈ [0, ∞)^d |∀ i ∈α∪α' : u_i > 1 } )
      =
      ρ_α∪α';

the first equality follows from (<ref>) and the last one from (<ref>). We obtain (<ref>) by expanding G_α = Z_α(1_α) using (<ref>) and working out [ G_α G_α' ] with the above identity.


 
Let α = {α_1, …, α_S }⊂{1, …, d} with S = |α|≥ 2 and such that μ(Δ_α) > 0. In view of (<ref>), we have κ_α = g_α( θ_α ) and κ_α = g_α( θ_α ) where θ_α = ( ρ_α, ρ_α∖α_1, …, ρ_α∖α_S ), θ_α = (ρ_α, ρ_α∖α_1, …, ρ_α∖α_S), and

    g_α( x_0, x_1, …, x_S) = x_0/∑_j=1^S x_j - (S-1) x_0,      x∈[0,∞)^1+S.

Let ∇ g_α(x) denote the gradient vector of g_α evaluated x and let ⟨ ·  ,  · ⟩ denote the scalar product in Euclidean space. Proposition <ref> combined with the delta method as in <cit.> gives, as n →∞,

    √(k) ( κ_α - κ_α )
      =
      √(k){ g_α( θ_α ) - g_α( θ_α ) }   =
      ⟨∇ g_α(θ_α),  √(k) (θ_α - θ_α) 
      ⟩
      +
      _(1)
      
       ⟨∇ g_α(θ_α),   
        (G_α,G_α∖α_1,…,G_α∖α_S)
      ⟩,

the weak convergence holding jointly in α by Slutsky's lemma and Proposition <ref>. The partial derivatives of g_α are

    ∂ g/∂ x_0( x ) 
         = ∑_j=1^S  x_j /{∑_j=1^S  x_j - (S-1) x_0 }^2, 
    ∂ g/∂ x_j( x )
         = - x_0/{∑_j=1^S  x_j - (S-1) x_0 }^2,
             j = 1, …, S.

Evaluating these at x = θ_α and using ∑_j ∈αρ_α∖ j - (S-1) ρ_α = μ(Δ_α) as in (<ref>) and (<ref>), we find that

    ⟨∇ g_α(θ_α),   
        (G_α,G_α∖α_1,…,G_α∖α_S)
      ⟩
      =
      μ(Δ_α)^-2{( ∑_j∈αρ_α∖ j) G_α
        - ρ_α∑_j∈α G_α∖ j},

in accordance to the right-hand side in (<ref>).

To calculate the asymptotic variance σ_κ,α^2, we introduce a few abbreviations: we write R_β = R_β(1_β) and W_β^∩ = W_β(1_β) = W(R_β) for ∅β∈ and we put W_j = W_{j}(1) for j = 1, …, d, so that G_α = W_α^∩ - ∑_j∈αρ̇_j,α W_j. We find

    H_α   = 
      ( ∑_i∈αρ_α∖ i) G_α 
      - 
      ρ_α∑_i∈α G_α∖ i
       =  
      (∑_i∈αρ_α∖ i)
      ( W^∩_α - ∑_j∈αρ̇_j,αW_j) 
      - 
      ρ_α∑_i∈α( W^∩_α∖ i - ∑_j∈α∖ iρ̇_j,α∖ i W_j).

From the proof of Proposition <ref>, recall that W(A ∪ B) = W(A) + W(B) almost surely for disjoint Borel sets A and B of [0, ∞]^d ∖{∞} bounded away from ∞; moreover, for such A and B, the variables W(A) and W(B) are uncorrelated. Since R_α∖ i is the disjoint union of R_α and R_α∖ i∖ R_α, we have therefore W_α∖ i^∩ = W_α^∩ + W(R_α∖ i∖ R_α) almost surely. In addition, ∑_i ∈αρ_α∖ i = μ(Δ_α) + (S-1)ρ_α by (<ref>) applied to ν = μ. As a consequence,

    H_α
      = 
      {μ(Δ_α) - ρ_α}
      W_α^∩
      - ρ_α∑_j ∈α W(R_α∖ j∖ R_α)
      +  ∑_j∈α K_α, j W_j

where

    K_α,j
      =
      ρ_α(∑_i∈α∖ jρ̇_j,α∖ i) 
      -
      (∑_i∈αρ_α∖ i)ρ̇_j,α,
           j ∈α.

The S+1 variables W_α^∩ = W(R_α) and W(R_α∖ j∖ R_α), j ∈α, are all uncorrelated, since they involve evaluating W at disjoint sets; W_j = W(R_{j}) is uncorrelated with W(R_α∖ j∖ R_α), for the same reason. Moreover, [ W_α^∩ W_j ] = Λ(R_α∩ R_{j}) = Λ(R_α) = ρ_α and similarly [ W(R_α∖ i∖ R_α) W_j ] = Λ(R_α∖ i∖ R_α) = ρ_α∖ i - ρ_α if i,j∈α and i  j. Hence


    (H_α)
      = 
      {μ(Δ_α) - ρ_α}^2 ρ_α  
      +
      ρ_α^2 ∑_j ∈α (ρ_α∖ j - ρ_α)
      + 
      ∑_i,j ∈α K_α,i K_α,jρ_{i,j}
    
      + 
      {μ(Δ_α) - ρ_α}ρ_α∑_j∈αK_α,j
      - 
      ρ_α∑_j∈αK_α,j∑_i∈α∖ j (ρ_α∖ i - ρ_α).

As ∑_j ∈α (ρ_α∖ j - ρ_α) = μ(Δ_α)- ρ_α and ∑_i ∈α∖ j(ρ_α∖ i - ρ_α) = μ(Δ_α) - ρ_α,j, we get


    (H_α)
      = {μ(Δ_α) - ρ_α}ρ_α{μ(Δ_α) + ∑_j∈αK_α,j}
        + ∑_i,j∈αK_α,iK_α,jρ_{i,j}
    
        -ρ_α∑_j ∈αK_α,j{μ(Δ_α) - ρ_α∖ j}.

Recall κ_α(x) in (<ref>). We have

    ∂/∂ x_j(
        1/κ_α(x))_x = 1_α   =
      ∂/∂ x_j(
        ∑_i ∈α r_α∖ i(x_α∖ i)/r_α(x))
      
       =
      ρ_α^-2(
        ρ_α∑_i ∈α∖ jρ̇_j, α∖ i
        -
        ρ̇_j, α∑_i ∈αρ_α∖ i)
      =
      ρ_α^-2 K_α, j.

It follows that κ̇_j,α = - ρ_α^-2 K_α, j / (1/κ_α)^2 = - K_α,j / μ(Δ_α)^2.
     
     
     
     
     
     
     
     
     
     
     
     
     
     











     
     
     
By (<ref>), we find that σ_κ,α^2 = μ(Δ_α)^-4( H_α ) is equal to the right-hand side of (<ref>).









     


  
  
  
  
  
  
  
  
  
We only need to prove that σ^2_κ,α = σ_κ, α^2 + _(1) as n →∞. In view of the expressions (<ref>) and (<ref>) for σ_κ, α^2 and σ_κ,α, it is enough to show that κ̇_j,α, n = κ̇_α,j + _(1), with κ̇_j,α, n in (<ref>); indeed, Corollary <ref> already gives consistency of μ(Δ_α) and ρ_β. Now since 2^-1 k^1/4{κ_α(
1_α + k^-1/4 e_j) -  κ_α( 1_α - k^-1/4 e_j ) }→κ̇_α,j as n →∞, a sufficient condition is that for some ϵ>0, 

    sup_[1-ϵ,2+ϵ]^α k^1/4|κ_α( x) - κ_α( x)|
      = _(1),      n →∞.

In turn, (<ref>) follows from weak convergence of k^1/2( κ_α - κ_α ) as n →∞ in the space ℓ^∞([1-, 1+]^α). In light of the expressions of κ_α and κ_α in terms of the (empirical) joint tail dependence functions r_β and r_β, respectively, weak convergence of k^1/2( κ_α - κ_α ) follows from Proposition <ref> and the functional delta method <cit.>. The calculations are similar to the ones for the Euclidean case in the proof of Proposition <ref>; an extra point to be noted is that if α is such that μ(Δ_α) > 0, then the denominator in the definition of κ_α(x) in (<ref>) is positive for all x in a neighbourhood of 1_α.














































Proposition <ref> implies, as n →∞, the weak convergence


    ( 
        √(k){r_α(2_α) - r_α(2_α},
        √(k){r_α(1_α) - r_α(1_α})
      (
        Z_α( 2_α), Z_α( 1_α)
      ).

Now η_α^P = g(r_α( 2_α), r_α ( 1_α)) and η_α = 1 = g(r_α( 2_α),  r_α ( 1_α)) = g(2ρ_α,ρ_α), with g(x,y) = log(2)/ log (x/y); note that the function r_α is homogeneous. Since the gradient of g is ∇ g(x,y) = log(2) (log(x/y))^-2 ( - x^-1, y^-1 ), the delta method gives

    √(k) (η^P - 1)      ⟨∇ g(2ρ_α, ρ_α) , ( Z_α( 2_α),  Z_α( 1_α) ) 
      ⟩
       =
      1/ρ_αlog 2 ⟨ ( -1/2, 1   ) ,  ( Z_α( 2_α), Z_α( 1_α)) ⟩
       = 
      -1/2 ρ_αlog 2 { Z_α( 2_α) - 2 Z_α( 1_α) }.

The first part of the assertion follows. As for the variance,

    (  Z_α( 2_α) - 2Z_α( 1_α) ) 
         = 
      ( Z_α( 2_α) )
      + 
      4( Z_α( 1_α) )
      -
      4 (Z_α( 2_α), Z_α( 1_α)),

  
  
  
  
  
  
  

The function r_α is homogeneous of order 1, so that  ∂_jr_α is constant along rays, that is, the function 0 < t ↦∂_j r_α(t   x) is constant. Moreover, the measure Λ is homogeneous of order 1 too. In view of (<ref>) and (<ref>), it follows that ( Z_α(t x) ) = t ( Z_α( x ) ) for t > 0; in particular ( Z_α(2_α) = 2 ( Z_α(1_α ). Further, ρ_α = (d r_α(t, …, t) / dt)_t=1 = ∑_j ∈αρ̇_j,α and thus 

    (Z_α( 1_α))
            =  ρ_α - 2∑_j∈αρ̇_j,αρ_α  +  ∑_j∈α∑_j'∈αρ̇_j,αρ̇_j',αρ_{j,j'}
        = ρ_α - 2 ρ_α^2 +  ∑_j∈α∑_j'∈αρ̇_j,αρ̇_j',αρ_{j,j'}.

The covariance term is

    (Z_α( 2_α), Z_α( 1_α))
        =ρ_α 
          - ∑_j∈αρ̇_j,αρ_α
          - ∑_j∈αρ̇_j,α
          r_α( 2_α∧_j) 
    
          + ∑_j∈α∑_j'∈αρ̇_j,αρ̇_j',αr_{j,j'}(2,1),

with 2_α∧_j as explained in the statement of the proposition. Since ∑_j∈αρ̇_j,α = ρ_α, we can simplify and find

    (  Z_α( 2_α) - 2Z_α( 1_α) ) 
         = 
      6 ( Z_α(1_α) ) - 4 (Z_α( 2_α), Z_α( 1_α)) 
       = 
      2 ρ_α   - 8 ρ_α^2 +  4∑_j∈αρ̇_j,α r_α( 2_α∧_j)
      
            
      + 
      ∑_j∈α∑_j'∈αρ̇_j,αρ̇_j',α[ 6 ρ_{j, j'} - 4 r_{j, j' }(2,1)].

Divide the right-hand side by (2 ρ_αlog 2)^2 to obtain (<ref>).




   To alleviate
  notations, ∅α⊂{1,…, d} is fixed
  and the subscript α is omitted throughout the proof.
  Introduce the tail empirical process Q_n(t) = T_(n - ⌊ kt ⌋) for 0 < t < n/k.
  The key is to represent the Hill estimator as a statistical tail functional <cit.> of Q_n, i.e., η^H = Θ(Q_n), where Θ is the map defined for any measurable function
  z: (0,1]→ as Θ(z) = ∫_0^1  log^+ { z(t) / z(1) } t when the integral is finite and Θ(z) = 0 otherwise. Let z_η: t∈(0,1]↦ t^-η denote the quantile function of a standard Pareto distribution with index 1/η; it holds that Θ(z_η) = η. The map Θ is scale invariant, i.e., Θ(t z) = Θ(z), t>0.
  
  The proof consists of three steps:
  
  
  * Introduce a function space D_η,h allowing
    to control Q_n(t) and z_η(t) as t→ 0. In this space and up to rescaling, 
    Q_n - z_η converges weakly
    to a Gaussian process.
  
  * Show that the map Θ is Hadamard
    differentiable at z_η tangentially to some well chosen
    subspace of D_η,h.
  
  * Apply the functional delta method to show that η^H = Θ(Q_n) is asymptotically normal and compute its asymptotic variance via the Hadamard derivative of Θ.
  

  

  
Step 1.
  Let ϵ>0 and h(t) = t^1/2+ ϵ,  t∈[0,1]. Then
  h∈ℋ, where
  
    ℋ = {z: [0,1]→| z  continuous, lim_t→ 0 z(t) t^-1/2 (loglog(1/t))^1/2 = 0
        }.
 Introduce the function space
  
    D_η, h = {z: [0,1]→|lim_t→ 0 t^η h(t)
        z(t) = 0 ;  t↦ t^ηh(t) z(t) ∈ D[0,1] },

  where D[0,1] is the space of càdlàg functions. Notice that
  z_η∈ D_η,h. Equip D_η, h with the seminorm z_η,h = sup_t∈(0,1]|t^η h(t) z(t)|. Let m = ⌈ n q^← (k/n)⌉, with ⌈ · ⌉ the ceil function, so that k/m→ρ;
  for self-consistency of the present paper, the roles of
  k and m are reversed compared to the notation in <cit.>. From <cit.>, we have, for all t_0>0, 
  in the space D_η,h, the weak convergence
  
    √(k)( m/nQ_n - z_η) ( η t^-(η +1)W̅(t))_t ∈ [0,t_0]

  where W̅ (t) = W̃( t_α), and W̃ is defined as in the
  statement of Proposition <ref>. Indeed, the
  process W̅ in the statement from <cit.>
  has same distribution as W_1( t_α) in
  the case ρ=0; recall that our ρ is denoted by l in <cit.>. Put U_i,j = 1 - F_j(X_i,j), and let U_(1),j≤…≤ U_(d),j be the order statistics of U_1,j, …, U_n,j. In the case ρ>0,
  W̅ equals in distribution W_dra ( t_α) where
  W_dra appears in Lemma 6.1 in the cited reference as the
  limit in distribution (for α= {1,2}), for
  x∈ E_α, of
  
    Δ_n,k,m ( x) 
           =  √(k)[  1/k∑_i=1^n {∀ j ∈α : U_i,j≤ U_(⌊ m x_j ⌋), j} - c( x)]  
    
        
        
    
    
           = 
        √(m/k)_→ρ^-1/2√(m)[ 
          1/m∑_i=1^n {∀ j ∈α : U_i,j≤ U_(⌊ m x_j ⌋),j}_r_n( x) with k replaced by m -
          r( x)k/mρ_→ 1].
 
  From Proposition <ref> and Slutsky's Lemma, we have Δ_n,k,mρ^-1/2 Z_α in ℓ^∞([0, 1]^α). Therefore, W_dra = ρ^-1/2 Z_α, as claimed.

  

  
Step 2.
  The right-hand side of
  (<ref>) 
  belongs to
  𝒞_h,η = { z ∈ D_η,h|z is continuous}.  To apply the functional delta-method   <cit.>, we must verify that
  the restriction of Θ to D̅_η, h is
  Hadamard-differentiable tangentially to 𝒞_η,h, with
  derivative Θ', where D̅_η, h is a subspace of
  D_η,h such that (Q_n ∈D̅_η,h)→ 1 as
  n→∞; see the remark following Condition 3 in <cit.>. Then it will follow
  from the scale invariance of Θ, the identities Θ(Q_n) = η^H and Θ(z_η) = η, and the weak convergence in (<ref>) that
  
    √(k)( η^H- η)
        =
        √(k)( Θ(m/nQ_n) - Θ(z_η) )
        Θ'[( η t^-(η +1)W̅(t))_t ∈ [0,1]]

  as n →∞. From <cit.>, the
  restriction of Θ to D̅_η,h, the subset of
  functions on D_η,h which are positive and non increasing, is
  indeed Hadamard differentiable; letting ν denote the measure ν (t) = t^η t + ϵ_1(t), with ϵ_1 a point mass at 1, the derivative is
  
    Θ'(z) = ∫_0^1 t^η z(t)  t - y(1) = ∫_[0,1] z(t)
        ν (t).

  
  

  
Step 3.
  The weak limit in (<ref>) is thus equal to ∫_[0,1]η t^-(η +1)W̅(t) ν(t).  From <cit.>, the latter random variable is centered Gaussian with variance
  
    σ^2 = ∬_[0,1]^2η^2 (st)^-(η +1)(W̅(s), W̅(t))ν(s)ν(t).

  
  
  By definition of ν and by symmetry of the covariance,

  
    σ^2 / η^2  = 2  ∫_s=0^1∫_t=0^s (st)^-1(W̅(s), W̅(t))  t  s_A
    
                             - 2 ∫_s=0^1 (W̅(s), W̅(1)) s^-1 s_B + (W̅(1)).

  For any s∈(0,1),
  
    ∫_t=0^s (W̅(s), W̅(t))(st)^-1 t
            =∫_u=0^1 (W̅(s), W̅(us))(su)^-1 u 
       = ∫_u=0^1 (W̅(1), W̅(u))(u)^-1 u  
        = B.

  The penultimate equality follows from (W̅(λ s),W̅( λ t)) = λ(W̅(s),W̅( t))
  for   λ>0 and s,t∈(0,1]. Therefore A=B and
  σ^2 = η^2 (W̅(1)), as required.

  

  
  




§ CLEF ALGORITHM AND VARIANTS



The CLEF algorithm is described at length in <cit.>. For completeness, its pseudo-code is provided below. The underlying idea is to iteratively construct pairs, triplets, quadruplets… of features that are declared `dependent' whenever κ_α≥ C for some user-defined tolerance level C > 0. Varying this criterion produces three variants of the original algorithm, namely CLEF-Asymptotic, CLEF-Peng, and CLEF-Hill. The pruning stage of the algorithm is the same for all three variants.








  This work was supported by a public grant as part of the
  Investissement d’avenir project, reference ANR-11-LABX-0056-LMH,
  LabEx LMH.




29



urlstyle

[Agrawal et al(1994)Agrawal, Srikant et al]agrawal1994fast
Agrawal R, Srikant R, et al (1994) Fast algorithms for mining association
  rules. In: Proc. 20th int. conf. very large data bases, VLDB, vol 1215, pp
  487–499

[Bacro and Toulemonde(2013)]bacro2013measuring
Bacro JN, Toulemonde G (2013) Measuring and modelling multivariate and spatial
  dependence of extremes. Journal de la Société Française de
  Statistique 154(2):139–155

[Bücher and Dette(2013)]bucher+d:2013
Bücher A, Dette H (2013) Multiplier bootstrap of tail copulas with
  applications. Bernoulli 19(5A):1655–1687

[Chiapino and Sabourin(2016)]chiapinofeature
Chiapino M, Sabourin A (2016) Feature clustering for extreme events analysis,
  with application to extreme stream-flow data. In: ECML-PKDD 2016, workshop
  NFmcp2016

[Coles et al(1999)Coles, Heffernan, and Tawn]coles1999dependence
Coles S, Heffernan J, Tawn J (1999) Dependence measures for extreme value
  analyses. Extremes 2(4):339–365

[De Haan and Zhou(2011)]HaanZhou11residual
De Haan L, Zhou C (2011) Extreme residual dependence for random vectors and
  processes. Advances in Applied Probability 43(01):217–242

[Draisma et al(2001)Draisma, Drees, Ferreira, and
  de Haan]draisma2001tail
Draisma G, Drees H, Ferreira A, de Haan L (2001) Tail dependence in
  independence. Eurandom preprint

[Draisma et al(2004)Draisma, Dress, Ferreira, and
  De Haan]draisma2004bivariate
Draisma G, Dress H, Ferreira A, De Haan L (2004) Bivariate tail estimation:
  dependence in asymptotic independence. Bernoulli pp 251–280

[Drees(1998a)]drees1998evindex
Drees H (1998a) A general class of estimators of the extreme value
  index. Journal of Statistical Planning and Inference 66(1):95–112

[Drees(1998b)]drees1998smooth
Drees H (1998b) On smooth statistical tail functionals.
  Scandinavian Journal of Statistics 25(1):187–210

[Eastoe and Tawn(2012)]eastoe2012subasymptotic
Eastoe EF, Tawn JA (2012) Modelling the distribution of the cluster maxima of
  exceedances of subasymptotic thresholds. Biometrika 99(1)

[Einmahl(1997)]einmahl1997poisson
Einmahl JH (1997) Poisson and Gaussian approximation of weighted local
  empirical processes. Stochastic Processes and Their Applications 70(1):31–58

[Einmahl et al(2012)Einmahl, Krajina, Segers et al]einmahl2012m
Einmahl JH, Krajina A, Segers J, et al (2012) An M-estimator for tail
  dependence in arbitrary dimensions. The Annals of Statistics 40(3):1764–1793

[Goix et al(2016)Goix, Sabourin, and Clémençon]goixsparse
Goix N, Sabourin A, Clémençon S (2016) Sparse representation of
  multivariate extremes with applications to anomaly ranking. In: Proceedings
  of the 19th AISTAT conference, pp 287–295

[Goix et al(2017)Goix, Sabourin, and
  Clémençon]goix2017sparse
Goix N, Sabourin A, Clémençon S (2017) Sparse representation of
  multivariate extremes with applications to anomaly detection. Journal of
  Multivariate Analysis 161:12–31

[Ledford and Tawn(1996)]ledford1996statistics
Ledford AW, Tawn JA (1996) Statistics for near independence in multivariate
  extreme values. Biometrika 83(1):169–187

[Peng(1999)]peng1999estimation
Peng L (1999) Estimation of the coefficient of tail dependence in bivariate
  extremes. Statistics & Probability Letters 43(4):399–409

[Pickands III(1975)]pickands1975statistical
Pickands III J (1975) Statistical inference using extreme order statistics. The
  Annals of Statistics pp 119–131

[Qi(1997)]qi1997almost
Qi Y (1997) Almost sure convergence of the stable tail empirical dependence
  function in multivariate extreme statistics. Acta Mathematicae Applicatae
  Sinica (English series) 13(2):167–175

[Ramos and Ledford(2009)]ramos2009new
Ramos A, Ledford A (2009) A new class of models for bivariate joint tails.
  Journal of the Royal Statistical Society: Series B 71(1):219–241

[Resnick(2007)]resnick:2007
Resnick SI (2007) Heavy-Tail Phenomena. Springer Series in Operations Research
  and Financial Engineering, Springer, New York

[Resnick(2008)]resnick2013extreme
Resnick SI (2008) Extreme Values, Regular Variation and Point Processes.
  Springer Series in Operations Research and Financial Engineering, Springer,
  New York, reprint of the 1987 original

[Rockafellar(1970)]rockafellar:1970
Rockafellar RT (1970) Convex Analysis. Princeton Mathematical Series, No. 28,
  Princeton University Press, Princeton, N.J.

[Schlather and Tawn(2003)]schlather2003dependence
Schlather M, Tawn JA (2003) A dependence measure for multivariate and spatial
  extreme values: Properties and inference. Biometrika 90(1):139–156

[Shorack and Wellner(2009)]shorack2009empirical
Shorack GR, Wellner JA (2009) Empirical processes with applications to
  statistics. SIAM

[Stephenson(2003)]stephenson2003simulating
Stephenson A (2003) Simulating multivariate extreme value distributions of
  logistic type. Extremes 6(1):49–59

[Tawn(1990)]tawn1990modelling
Tawn JA (1990) Modelling multivariate extreme value distributions. Biometrika
  77(2):245–253

[van der Vaart(1998)]van2000asymptotic
van der Vaart AW (1998) Asymptotic Statistics, Cambridge Series in Statistical
  and Probabilistic Mathematics, vol 3. Cambridge University Press, Cambridge

[van der Vaart and Wellner(1996)]van1996weak
van der Vaart AW, Wellner JA (1996) Weak Convergence and Empirical Processes.
  Springer, New York



