Extremal attractors of Liouville copulas
Léo R. Belzile∗ and Johanna G. Nešlehová†

arXiv:1704.03377v2 [math.ST] 6 Jul 2017

This is a copyedited, author-produced version of an article accepted for publication following peer-review in the
Journal of Multivariate Analysis, an Elsevier publication, c 2017. This manuscript version is made available under
the CC-BY-NC-ND 4.0 license. Please cite as
Belzile, L. R and J. G. Nešlehová. Extremal attractors of Liouville copulas, Journal of Multivariate
Analysis (2017), 160C, pp. 68–92. doi:10.1016/j.jmva.2017.05.008
Abstract
Liouville copulas introduced in [31] are asymmetric generalizations of the ubiquitous Archimedean copula class.
They are the dependence structures of scale mixtures of Dirichlet distributions, also called Liouville distributions. In
this paper, the limiting extreme-value attractors of Liouville copulas and of their survival counterparts are derived.
The limiting max-stable models, termed here the scaled extremal Dirichlet, are new and encompass several existing
classes of multivariate max-stable distributions, including the logistic, negative logistic and extremal Dirichlet. As
shown herein, the stable tail dependence function and angular density of the scaled extremal Dirichlet model have a
tractable form, which in turn leads to a simple de Haan representation. The latter is used to design efficient algorithms
for unconditional simulation based on the work of [10] and to derive tractable formulas for maximum-likelihood
inference. The scaled extremal Dirichlet model is illustrated on river flow data of the river Isar in southern Germany.

1. Introduction
Copula models play an important role in the analysis of multivariate data and find applications in many areas, including
biostatistics, environmental sciences, finance, insurance, and risk management. The popularity of copulas is rooted
in the decomposition of Sklar [39], which is at the heart of flexible statistical models and various measures, concepts
and orderings of dependence between random variables. According to Sklar’s result, the distribution function of any
random vector X = (X1 , . . . , Xd ) with continuous univariate margins F1 , . . . , Fd satisfies, for any x1 , . . . , xd ∈ R,
Pr(X1 ≤ x1 , . . . , Xd ≤ xd ) = C{F1 (x1 ), . . . , Fd (xd )},
for a unique copula C, i.e., a distribution function on [0, 1]d whose univariate margins are standard uniform. Alternatively, Sklar’s decomposition also holds for survival functions, i.e., for any x1 , . . . , xd ∈ R,
Pr(X1 > x1 , . . . , Xd > xd ) = Ĉ{F̄1 (x1 ), . . . , F̄d (xd )},
where F̄1 , . . . , F̄d are the marginal survival functions and Ĉ is the survival copula of X, related to the copula of X as
follows. If U is a random vector distributed as the copula C of X, Ĉ is the distribution function of 1 − U.
In risk management applications, the extremal behavior of copulas is of particular interest, as it describes the dependence between extreme events and consequently the value of risk measures at high levels. Our purpose is to study
the extremal behavior of Liouville copulas. The latter are defined as the survival copulas of Liouville distributions
[14, 17, 38], i.e., distributions of random vectors of the form R Dα , where R is a strictly positive random variable
independent of the Dirichlet random vector Dα = (D1 , . . . , Dd ) with parameter vector α = (α1 , . . . , αd ). Liouville
copulas were proposed by McNeil and Nešlehová [31] in order to extend the widely used class of Archimedean
copulas and create dependence structures that are not necessarily exchangeable. The latter property means that for
∗ École

Polytechnique Fédérale de Lausanne, EPFL-SB-MATH-STAT, Station 8, CH-1015 Lausanne, Switzerland. leo.belzile@epfl.ch
of Mathematics and Statistics, McGill University, 805 rue Sherbrooke Ouest, Montréal, Québec, H3A 0B9, Canada.
johanna.neslehova@mcgill.ca
† Department

2

L. R. Belzile and J. G. Nešlehová

any u1 , . . . , ud ∈ [0, 1] and any permutation π of the integers 1, . . . , d, C(u1 , . . . , ud ) = C(uπ(1) , . . . , uπ(d) ). When
α = 1d ≡ (1, . . . , 1), Dα = D1d is uniformly distributed on the unit simplex
Sd = {x ∈ [0, 1]d : x1 + · · · + xd = 1}.

(1)

In this special case, one recovers Archimedean copulas. Indeed, according to [30], the latter are the survival copulas
of random vectors R D1d , where R is a strictly positive random variable independent of D1d . When α , 1d , the survival
copula of RDα is not Archimedean anymore. It is also no longer exchangeable, unless α1 = · · · = αd .
In this article, we determine the extremal attractor of a Liouville copula and of its survival counterpart. As a byproduct, we also obtain the lower and upper tail dependence coefficients of Liouville copulas that quantify the strength
of dependence at extreme levels [25]. These results are complementary to [21], where the upper tail order functions
of a Liouville copula and its density are derived when α1 = · · · = αd , and to [19], where the extremal attractor of
RDα is derived when R is light-tailed. The extremal attractors of Liouville copulas are interesting in their own right.
Because non-exchangeability of Liouville copulas carries over to their extremal limits, the latter can be used to model
the dependence between extreme risks in the presence of causality relationships [15]. The limiting extreme-value
models can be embedded in a single family, termed here the scaled extremal Dirichlet, whose members are new, nonexchangeable generalizations of the logistic, negative logistic, and Coles–Tawn extremal Dirichlet models given in
[7]. We examine the scaled extremal Dirichlet model in detail and derive its de Haan spectral representation. The latter
is simple and leads to feasible stochastic simulation algorithms and tractable formulas for likelihood-based inference.
The article is organized as follows. The extremal behavior of the univariate margins of Liouville distributions is
first studied in Section 2. The extremal attractors of Liouville copulas and their survival counterparts are then derived
in Section 3. When α is integer-valued, the results of [27, 31] lead to closed-form expressions for the limiting stable
tail dependence functions, as shown in Section 4. Section 5 is devoted to a detailed study of the scaled extremal
Dirichlet model. In Section 6, the de Haan representation is derived and used for stochastic simulation. Estimation
is investigated in Section 7, where expressions for the censored likelihood and the gradient score are also given. An
illustrative data analysis of river flow of the river Isar is presented in Section 8, and the paper is concluded by a
discussion in Section 9. Lengthy proofs are relegated to the Appendices.
In what follows, vectors in Rd are denoted by boldface letters, x = (x1 , . . . , xd ); 0d and 1d refer to the vectors
(0, . . . 0) and (1, . . . , 1) in Rd , respectively. Binary operations such as x + y or a · x, xa are understood as componentwise operations. k · k stands for the `1 -norm, viz. kxk = |x1 | + · · · + |xd |, ⊥
⊥ for statistical independence. For any x, y ∈ R,
let x ∧ y = min(x, y) and x ∨ y = max(x, y). The Dirac delta function Ii j is 1 if i = j and zero otherwise. Finally, Rd+ is
the positive orthant [0, ∞)d and for any x ∈ R, x+ denotes the positive part of x, max(0, x).

2. Marginal extremal behavior
A Liouville random vector X = RDα is a scale mixture of a Dirichlet random vector Dα = (D1 , . . . , Dd ) with
parameters α = (α1 , . . . , αd ) > 0d . In what follows, R is referred to as the radial variable of X and ᾱ denotes the sum
of the Dirichlet parameters, viz. ᾱ = kαk = α1 + · · · + αd . Recall that Dα has the same distribution as Z/kZk, where
Zi ∼ Ga(αi , 1), i = 1, . . . , d are independent Gamma variables with scaling parameter 1. The margins of X are thus
scale mixtures of Beta distributions, i.e., for i = 1, . . . , d, Xi = RDi with Di ∼ Beta(αi , ᾱ − αi ).
As a first step towards the extremal behavior of Liouville copulas, this section is devoted to the extreme-value
properties of the univariate margins of the vectors X and 1/X, where X is a Liouville random vector with parameters
α and a strictly positive radial part R, i.e., such that Pr(R ≤ 0) = 0. To this end, recall that a univariate random variable
X with distribution function F is in the maximum domain of attraction of a non-degenerate distribution F0 , denoted
F ∈ M(F0 ) or X ∈ M(F0 ), if and only if there exist sequences of reals (an ) and (bn ) with an > 0, such that, for any
x ∈ R,
lim F n (an x + bn ) = F0 (x).

n→∞

By the Fisher–Tippett Theorem, F0 must be, up to location and scale, either the Fréchet (Φρ ), the Gumbel (Λ) or the
Weibull distribution (Ψρ ) with parameter ρ > 0. Further recall that a measurable function f : R+ → R+ is called
regularly varying with index ρ ∈ (−∞, ∞), denoted f ∈ Rρ , if for any x > 0, f (tx)/ f (t) → xρ as t → ∞. If ρ = 0, f is
called slowly varying. For more details and conditions for F ∈ M(F0 ), see, e.g., [12, 35].

Extremal attractors of Liouville copulas

3

Because the univariate margins of X are scale mixtures of Beta distributions, their extremal behavior, detailed in
Proposition 1, follows directly from Theorems 4.1, 4.4. and 4.5 in [20].
Proposition 1 Let X = RDα be a Liouville random vector with parameters α = (α1 , . . . , αd ) and a strictly positive
radial variable R, i.e., Pr(R ≤ 0) = 0. Then the following statements hold for any ρ > 0:
(a) R ∈ M(Φρ ) if and only if Xi ∈ M(Φρ ) for all i = 1, . . . , d.
(b) R ∈ M(Λ) if and only if Xi ∈ M(Λ) for all i = 1, . . . , d.
(c) R ∈ M(Ψρ ) if and only if Xi ∈ M(Ψρ+ᾱ−αi ) for all i = 1, . . . , d.
Proposition 1 implies that the univariate margins of X are all in the domain of attraction of the same distribution if the
latter is Gumbel or Fréchet. This is not the case when R is in the Weibull domain of attraction. Note also that there are
cases not covered by Proposition 1, in which the univariate margins Xi are in the Weibull domain while R is not in the
domain of attraction of any extreme-value distribution. For example, when d = 2, α = (1, 1) and R = 1 almost surely,
the margins of X are standard uniform and hence in the maximum domain of attraction of Ψ1 ; see Example 3.3.15 in
[12]. At the same time, R is clearly neither in the Weibull, nor the Gumbel, nor the Fréchet domain of attraction.
In subsequent sections, we shall also need the extremal behavior of the univariate margins of 1/X. The proposition
below shows that the latter is determined by the properties of 1/R. In contrast to Proposition 1, however, the univariate
margins of 1/X are always in the Fréchet domain. The proof may be found in A.
Proposition 2 Let X = RDα be a Liouville random vector with parameters α = (α1 , . . . , αd ) and a strictly positive
radial variable R with Pr(R ≤ 0) = 0. The following statements hold for any i = 1, . . . , d.
(a) If 1/R ∈ M(Φρ ) for ρ ∈ (0, αi ], then 1/Xi ∈ M(Φρ ).
(b) If E(1/Rαi +ε ) < ∞ for some ε > 0, then 1/Xi ∈ M(Φαi ).

3. Extremal behavior of Liouville copulas
In this section, we will identify the extremal behavior of a Liouville random vector X = RDα and of the random
vector 1/X, assuming that Pr(R ≤ 0) = 0. As a by-product, we will obtain the extremal attractors of Liouville copulas
and their survival counterparts. To this end, recall that a random vector Y with joint distribution function H is in the
maximum domain of attraction of a non-degenerate distribution function H0 , in notation H ∈ M(H0 ) or Y ∈ M(H0 ),
iff there exist sequences of vectors (an ) in (0, ∞)d and (bn ) in Rd such that for all x ∈ Rd ,
lim H n (an x + bn ) = H0 (x).

n→∞

When the univariate margins F1 , . . . , Fd of H are continuous, H ∈ M(H0 ) holds if and only if Fi ∈ M(F0i ) for all
i = 1, . . . , d, where F01 , . . . , F0d are the univariate margins of H0 , and further if the unique copula C of H is in the
domain of attraction of the unique copula C0 of H0 , denoted C ∈ M(C0 ), i.e., iff for all u ∈ [0, 1]d ,
lim C n (u1/n ) = C0 (u).

n→∞

In particular, the univariate margins of the max-stable distribution H0 must each follow a generalized extreme-value
distribution, and C0 must be an extreme-value copula. This means that for all u ∈ [0, 1]d ,
C0 (u) = exp[−`{− log(u1 ), . . . , − log(ud )}],

(2)

where ` : Rd+ → [0, ∞) is a stable tail dependence function, linked to the so-called exponent measure ν viz.
ν{[0d , x)c } = `(1/x), see, e.g., [35]. The latter Rcan be characterized through an angular (or spectral) probability
measure σd on Sd given in Eq. (1) which satisfies S wi dσd (w) = 1/d for all i = 1, . . . , d. For all x ∈ Rd+ , one has
d

`(x) = d

Z
Sd

max(w1 x1 , . . . , wd xd ) dσd (w).

(3)

4

L. R. Belzile and J. G. Nešlehová

Because ` is homogeneous of order 1, i.e., for any c > 0 and x ∈ Rd+ , `(cx) = c`(x), C0 can also be expressed via the
Pickands dependence function A : Sd → [0, ∞) related to ` through `(x) = kxkA(x/kxk). Then at any u ∈ [0, 1]d ,
)#
"
(
log(ud )
log(u1 )
,...,
.
C0 (u) = exp log(u1 · · · ud )A
log(u1 · · · ud )
log(u1 · · · ud )
When d = 2, it is more common to define the Pickands dependence function A : [0, 1] → [0, 1] through `(x1 , x2 ) =
(x1 + x2 )A{x2 /(x1 + x2 )} so that, for all u1 , u2 ∈ [0, 1],
"
(
)#
log(u2 )
C0 (u1 , u2 ) = exp log(u1 u2 )A
.
(4)
log(u1 u2 )
Now consider a Liouville vector X = RDα with a strictly positive radial variable. Theorem 1 specifies when
X ∈ M(H0 ) and identifies H0 . While part (a) follows from regular variation of X, parts (b) and (c) are special cases
of the results discussed in Section 2.2 in [19]. Details of the proof may be found in B.
Theorem 1 Let X = RDα , Dα = (D1 , . . . , Dd ), α > 0d , and Pr(R ≤ 0) = 0. Then the following statements hold.
(a) If R ∈ M(Φρ ) for some ρ > 0, then X ∈ M(H0 ), where H0 is a multivariate extreme-value distribution with
univariate margins F0i = Φρ , i = 1, . . . , d, and a stable tail dependence function given, for all x ∈ Rd+ , by



ρ

Γ(αd )xd Dρd 
Γ(ᾱ + ρ) 
 Γ(α1 )x1 D1

E max 
,...,
`(x) =
 .

Γ(ᾱ)
Γ(α1 + ρ)
Γ(αd + ρ) 
(b) If R ∈ M(Λ), then X ∈ M(H0 ), where for all x ∈ Rd , H0 (x) =

Qd
i=1

Λ(xi ).

(c) If R ∈ M(Ψρ ) for some ρ > 0, then X ∈ M(H0 ), where for all x ∈ Rd , H0 (x) =

Qd
i=1

Ψρ+ᾱ−αi (xi ).

The next result, also proved in B, specifies the conditions under which 1/X ∈ M(H0 ) and gives the form of the limiting
extreme-value distribution H0 .
Theorem 2 Let X = R Dα , Dα = (D1 , . . . , Dd ), α > 0d , and assume that Pr(R ≤ 0) = 0. Let αM = max(α1 , . . . , αd ).
The following cases can be distinguished:
P
(a) If 1/R ∈ M(Φρ ) for ρ ∈ (0, αM ], set I1 = {i : αi ≤ ρ}, I2 = {i : αi > ρ} and ᾱ2 = i∈I2 αi . Then
1/X ∈ M(H0 ), where the univariate margins of H0 are F0i = Φρ∧αi , i = 1, . . . , d, and the stable tail dependence
function is given, for all x ∈ Rd+ , by






−ρ 


X
e−ρ 




Γ(ᾱ2 − ρ) 
Γ(ᾱ − ρ) 
 Γ(αi )xi Di 
 X
 Γ(αi )xi D
i 
 ,
E max 
xi +
E max 
xi +
`(x) =
 =






i∈I2  Γ(αi − ρ) 
i∈I2  Γ(αi − ρ) 
Γ(ᾱ)
Γ(ᾱ2 )
i∈I1

i∈I1

ei , i ∈ I2 ) is a Dirichlet random vector with parameters (αi , i ∈ I2 ) if |I2 | > 1 and D
ei ≡ 1 if I2 = {i}.
where (D


Q
(b) If E 1/Rβ < ∞ for β > αM , then 1/X ∈ M(H0 ), where for all x ∈ Rd , H0 (x) = di=1 Φαi (xi ).
Remark 1 Note that in the case of asymptotic independence between the components of X (Theorem 1 (b–c)) or
1/X (Theorem 2 (b)), dependence between component-wise maxima of finitely many vectors may still be present.
Refinements of asymptotic independence are then needed, but these considerations surpass the scope of this paper.
One option would be to consider triangular arrays as in [23]; extremes of arrays of Liouville vectors can be obtained
as a special case of extremes of arrays of weighted Dirichlet distributions developed in [18]. Another avenue worth
exploring might be the limits of scaled sample clouds, as in [2] and [32].
The stable tail dependence functions appearing in Theorems 1 and 2 will be investigated in greater detail in the
subsequent sections. Before proceeding, we introduce the following terminology, emphasizing that they can in fact be
embedded in one and the same parametric class.

Extremal attractors of Liouville copulas

5

Definition 1 For any α > 0 and ρ ∈ (−α, ∞), let c(α, ρ) = Γ(α + ρ)/Γ(α) denote the rising factorial. For d ≥ 2
and α1 , . . . , αd > 0 and let (D1 , . . . , Dd ) denote a Dirichlet random vector with parameters α = (α1 , . . . , αd ) and set
ᾱ = α1 + · · · + αd . For any − min(α1 , . . . , αd ) < ρ < ∞, the scaled extremal Dirichlet stable tail dependence function
with parameters ρ and α is given, for all x ∈ Rd+ , by



ρ

xd Dρd 


 x 1 D1
,
.
.
.
,
`D (x; ρ, α) = c(ᾱ, ρ)E max 
 ,
 c(α1 , ρ)
c(αd , ρ) 

(5)

when ρ , 0 and by max(x1 , . . . , xd ) when ρ = 0. For any ρ > 0, the positive scaled extremal Dirichlet stable tail
dependence function `pD with parameters ρ and α is given, for all x ∈ Rd+ , by `pD (x; ρ, α) = `D (x; ρ, α), while for any
0 < ρ < min(α1 , . . . , αd ), the negative scaled extremal Dirichlet stable tail dependence function `nD is given, for all
x ∈ Rd+ , by `nD (x; ρ, α) = `D (x; −ρ, α).
Remark 2 As will be seen in Section 5, distinguishing between the positive and negative scaled extremal Dirichlet
models makes the discussion of their properties slightly easier because the sign of ρ impacts the shape of the corresponding angular measure. When ρ → 0, `D (x; ρ, α) becomes max(x1 , . . . , xd ), the stable tail dependence function
corresponding to comonotonicity, while when ρ → ∞, `D (x; ρ, α) becomes x1 + · · · + xd , the stable tail dependence
function corresponding to independence. Note also that ρ ∈ (−∞, ∞) can be allowed, with the convention that all
variables whose indices i are such that ρ ≤ −αi are independent, i.e., `nD is then of the form given in Theorem 2 (a).
From Theorems 1 and 2, we can now easily deduce the extremal behavior of Liouville copulas and their survival
counterparts. To this end, recall that a Liouville copula C is defined as the survival copula of a Liouville random
vector X = RDα with Pr(R ≤ 0) = 0. The following corollary follows directly from Theorem 2 upon noting that C is
also the unique copula of 1/X.
Corollary 1 Let C be the unique survival copula of a Liouville random vector X = RDα with Pr(R ≤ 0) = 0. Let
αM = max(α1 , . . . , αd ) and set I1 = {i : αi ≤ ρ}, I2 = {i : αi > ρ}. Then the following statements hold.
(a) If 1/R ∈ M(Φρ ) for ρ ∈ (0, αM ] and |I2 | > 1, then C ∈ M(C0 ), where C0 is an extreme-value copula of the form
(2) whose stable tail dependence function is given, for all x ∈ Rd+ , by
X
xi + `nD (x{2} ; ρ, α{2} ),
`(x) =
i∈I1

where and x{2} = (xi , i ∈ I2 ), α{2} = (αi , i ∈ I2 ). If 1/R ∈ M(Φρ ) for ρ ∈ (0, αM ] and |I2 | ≤ 1, then C ∈ M(Π),
where Π is the independence copula given, for all u ∈ [0, 1]d , by Π(u) = u1 · · · ud .


(b) If E 1/Rβ < ∞ for β > αM , then C ∈ M(Π).
Remark 3 Observe that Corollary 1 (a) in particular implies that when d = 2, α1 < α2 and 1/R ∈ M(Φρ ) for
α1 ≤ ρ < α2 , C ∈ M(Π). Also note that when α1 = · · · = αd ≡ α and 1/R ∈ M(Φρ ) for ρ ∈ (0, α), the result in
Corollary 1 (a) can be derived from formula (5) in Proposition 3 in [21] by relating the tail order function to the stable
tail dependence function when the tail order equals 1.
The survival counterpart Ĉ of a Liouville copula C is given as the distribution function of 1 − U, where U is a
random vector distributed as C. As C is the unique survival copula of X, Ĉ is the unique copula of X. The following
result thus follows directly from Theorem 1.
Corollary 2 Let Ĉ be the unique copula of a Liouville random vector X = R Dα with Pr(R ≤ 0) = 0. Then the
following statements hold.
(a) If R ∈ M(Φρ ) for ρ > 0, then Ĉ ∈ M(C0 ), where C0 is an extreme-value copula of the form (2) with the positive
scaled extremal Dirichlet stable tail dependence function given, for all x ∈ Rd+ , by `pD (x; ρ, α).
(b) If R ∈ M(Λ) or R ∈ M(Ψρ ) with ρ > 0, then Ĉ ∈ M(Π), where Π is the independence copula.

6

L. R. Belzile and J. G. Nešlehová

4. The case of integer-valued Dirichlet parameters
When α is integer-valued, Liouville distributions are particularly tractable because their survival function is explicit.
In this section, we will use this fact to derive closed-form expressions for the positive and negative scaled extremal
Dirichlet stable tail dependence functions. To this end, first recall the notion of the Williamson transform. The latter
is related to Weyl’s fractional integral transform and was used to characterize d-monotone functions in [43]; it was
adapted to non-negative random variables in [30].
Definition 2 Let X be a non-negative random variable with distribution function F, and let k ≥ 1 be an arbitrary
integer. The Williamson k-transform of X is given, for all x > 0, by
Wk F(x) =

∞

Z
x



x k−1
x k−1
dF(r) = E 1 −
.
1−
r
X +

For any k ≥ 1, the distribution of a positive random variable X is uniquely determined by its Williamson k-transform,
the formula for the inverse transform being explicit [30, 43]. If ψ = Wk F, then, for all x > 0,
F(x) = W−1
k ψ(x) = 1 −

k−2
X
(−1) j x j ψ( j) (x)

j!

j=0

−

(−1)k−1 xk−1 ψ+(k−1) (x)
,
(k − 1)!

where for j = 1, . . . , k − 2, ψ( j) is the jth derivative of ψ and ψ(k−1)
is the right-hand derivative of ψ(k−2) . These deriva+
tives exist because a Williamson k-transform ψ is necessarily k-monotone [43]. This means that ψ is differentiable up
to order k − 2 on (0, ∞) with derivatives satisfying (−1) j ψ( j) ≥ 0 for j = 0, . . . , k − 2 and such that (−1)k−2 ψ(k−2) is
non-increasing and convex on (0, ∞). Moreover, ψ(x) → 0 as x → ∞ and if F(0) = 0, ψ(x) → 1 and x → 0.
Now let C be a Liouville copula corresponding to a Liouville random vector X = RDα with integer-valued
parameters α = (α1 , . . . , αd ) and a strictly positive radial part R, i.e., Pr(R ≤ 0) = 0. Let ψ be the Williamson
ᾱ-transform of R and set Iα = {0, . . . , α1 − 1} × · · · × {0, . . . , αd − 1}. By Theorem 2 in [31], one then has, for all x ∈ Rd+ ,
Pr(X > x) = H̄(x) =

X

ψ( j1 +···+ jd ) (x1 + · · · + xd ) Y ji
xi .
j1 ! · · · jd !
i=1
d

(−1) j1 +···+ jd

( j1 ,..., jd )∈Iα

(6)

In particular, the margins of X have survival functions satisfying, for all x > 0 and i = 1, . . . , d,
Pr(Xi > x) = H̄i (x) =

α
i −1
X
j=0

(−1) j x j ψ( j) (x)
= 1 − W−1
αi ψ(x).
j!

(7)

By Sklar’s Theorem for survival functions, the Liouville copula C is given, for all u ∈ [0, 1]d , by
C(u) = H̄{H̄1−1 (u1 ), . . . , H̄d−1 (ud )}.
Although this formula is not explicit, it is clear from Equations (6) and (7) that C depends on the distribution of X
only through the Williamson ᾱ-transform ψ of R and the Dirichlet parameters α. For this reason, we shall denote
the Liouville copula in this section by Cψ,α and refer to ψ as its generator, reiterating that ψ must be an ᾱ-monotone
function satisfying ψ(1) = 0 and ψ(x) → 0 as x → ∞. When α = 1d , Cψ,1 is the Archimedean copula with generator
ψ, given, for all u ∈ [0, 1]d by Cψ,1 (u) = ψ{ψ−1 (u1 ) + · · · + ψ−1 (ud )}. Because the relationship between ψ and R is
one-to-one [30, Proposition 3.1], we will refer to R as the radial distribution corresponding to ψ.
Now suppose that 1/R ∈ M(Φρ ) with ρ ∈ (0, 1). By Theorem 2 in [27], this condition is equivalent to 1 − ψ(1/·) ∈
R−ρ . It further follows from Corollary 1 (a) that Cψ,α ∈ M(C0 ) where C0 is an extreme-value copula with the negative
scaled extremal Dirichlet stable tail dependence function `nD (·; ρ, α). This is because ρ < 1 ≤ min(α1 , . . . , αd ) so that
I1 = ∅ in Corollary 1 (a). Eq. (6) and the results of [27] can now be used to derive the following explicit expression
for `nD , as detailed in C.

7

Extremal attractors of Liouville copulas

Proposition 3 Let Cψ,α be a Liouville copula with integer-valued parameters α = (α1 , . . . , αd ) and generator ψ. If
1 − ψ(1/·) ∈ R−ρ for some ρ ∈ (0, 1), then Cψ,α ∈ M(C0 ), where C0 is an extreme-value copula with scaled negative
extremal Dirichlet stable tail dependence function `nD as given in Definition 1. Furthermore, for all x ∈ Rd+ ,


 n
o1/ρ  ji 
 ρ 
 d (
)
xi


d
1/ρ
X Γ( j1 + · · · + jd − ρ) Y

 
 
X
xj
1
c(αi ,−ρ)

  .
 1 − ρ
`nD (x; ρ, α) = Γ(1 − ρ) 

c(α j , −ρ)
Γ(1 − ρ)
Γ( ji + 1)  Pd n xk o1/ρ  
( j1 ,..., jd )∈Iα
i=1
j=1
k=1 c(α ,−ρ)
( j1 ,..., jd ),(0,...,0)

k

When α = 1d , the index set Iα reduces to the singleton {0}, and the expression for `nD given in Proposition 3 simplifies,
for all x ∈ Rd+ , to the stable tail dependence function of the Gumbel–Hougaard copula, viz.

`nD (x; ρ, 1d ) = x11/ρ + · · · + xd1/ρ ρ .
The Liouville copula Cψ,1d , which is the Archimedean copula with generator ψ, is thus indeed in the domain of
attraction of the Gumbel–Hougaard copula with parameter 1/ρ, as shown, e.g., in [6, 27].
Remark 4 When α = 1d and 1 − ψ(1/·) ∈ R−1 , it is shown in Proposition 2 of [27] that Cψ,1 is in the domain of
attraction of the independence copula. However, when α is integer-valued but such that max(α1 , . . . , αd ) > 1, regular
variation of 1 − ψ(1/·) does not suffice to characterize those cases in Corollary 1 that are not covered by Proposition 3.
This is because by Theorem 2 of [27], 1/R ∈ M(Φρ ) for ρ ≥ 1, 1/R ∈ M(Λ) and 1/R ∈ M(Ψρ ) for ρ > 0 all imply
that 1 − ψ(1/·) ∈ R−1 . At the same time, by Corollary 1, Cψ,α ∈ M(Π) clearly does not hold in all these cases.
Next, let Ĉψ,α be the survival copula of a Liouville copula Cψ,α , i.e., the distribution function of 1 − U, where U is a
random vector with distribution function Cψ,α . The results of [27] can again be used to restate the conditions under
which Ĉψ,α ∈ M(C0 ) in terms of ψ and to give an explicit expression for the stable tail dependence function of C0 .
Proposition 4 Let Ĉψ,α be the survival copula of a Liouville copula Cψ,α with integer-valued parameters α and a
generator ψ. Then the following statements hold.
(a) If ψ ∈ R−ρ for some ρ > 0, then Ĉψ,α ∈ M(C0 ), where C0 has a positive scaled extremal Dirichlet stable tail
dependence function `pD as given in Definition 1. The latter can be expressed, for all x ∈ Rd+ , as

−1/ρ −ρ


d
k


X
X
X



x
Γ(1
+
ρ)
i


j

(−1)k+1 
`pD (x; ρ, α) =

 ×





Γ(ρ) k=1 1≤i <···<i ≤d
 j=1 c(αi j , ρ) 
1

k

 n
o−1/ρ  jm

xim

k 
Y


Γ( j1 + · · · + jk + ρ)
c(αim , ρ)




 .
−1/ρ 
 Pk
xi j
j
!
·
·
·
j
!
1
k

( j1 ,..., jk )∈I(αi ,...,αi )
m=1
1
k
X

j=1

c(αi j , ρ)

(b) If ψ ∈ M(Λ) or ψ ∈ M(Ψρ ) for some ρ > 0, C̄ψ,α ∈ M(Π), where Π is the independence copula.
When α = 1d , the expression for `pD in part (a) of Proposition 4 simplifies, for all x ∈ Rd+ , to


X
X −1/ρ −ρ
pD
|A|+1 
` (x; ρ, 1d ) =
(−1)
x

 ,
i

A⊆{1,...,d}, A,∅

i∈A

which is the stable tail dependence function of the Galambos copula [24]. When ψ ∈ R−ρ for some ρ > 0, Ĉψ,1d is thus
indeed in the domain of attraction of the Galambos copula, as shown, e.g., in [27].

5. Properties of the scaled extremal Dirichlet models
In this section, the scaled extremal Dirichlet model with stable tail dependence function given in Definition 1 is
investigated in greater detail. In Section 5.1 we derive formulas for the so-called angular density and relate the
positive and negative scaled extremal Dirichlet models to classical classes of stable tail dependence functions. In
Section 5.2 we focus on the bivariate case and derive explicit expressions for the stable tail dependence functions and,
as a by-product, obtain formulas for the tail dependence coefficients of Liouville copulas.

8

L. R. Belzile and J. G. Nešlehová
scaled Dirichlet, ρ < 0

2.0
1.5
0.0

0.0

0.5

1.0

Angular density

2.0
1.0

Angular density

3.0

2.5

3.0

scaled Dirichlet, ρ > 0

0.0

0.2

0.4

0.6
w

0.8

1.0

0.0

0.2

0.4

0.6

0.8

1.0

w

Figure 1: Angular density of the scaled extremal Dirichlet model. Left panel: ρ = 4/5 and α = (2, 1/2) (black full),
ρ = 1/4 and α = (1/10, 1/10) (red dashed), ρ = 1/4 and α = (1/2, 1/2) (blue dotted). Right panel: ρ = −1/4 and
α = (2, 1/2) (black full), α = (2/5, 2/5) (red dashed) and α = (1/2, 1/2) (blue dotted).
5.1. Angular density
The first property worth noting is that the positive and negative scaled extremal Dirichlet models are closed under
marginalization. Indeed, letting xi → 0 for some arbitrary 1 ≤ i ≤ d, we can easily derive from Lemma 2 that for
any α > 0d , ρ > 0, and any x ∈ Rd+ , `pD (x; ρ, α) → `pD (x−i ; ρ, α−i ) as xi → 0, where for any y ∈ Rd , y−i denotes the
vector (y1 , . . . , yi−1 , yi+1 , . . . , yd ). Similarly, for any 1 ≤ i ≤ d, α > 0d , 0 < ρ < min(α1 , . . . , αd ), and any x ∈ Rd+ ,
`nD (x; ρ, α) → `nD (x−i ; ρ, α−i ) as xi → 0.
Since none of the scaled extremal Dirichlet models places mass on the vertices or facets of the simplex Sd when
ρ , 0, the density of the angular measure σd completely characterizes the stable tail dependence function and hence
also the associated extreme-value copula. This so-called angular density of the scaled extremal Dirichlet models is
given below and derived in D.
Proposition 5 Let d ≥ 2 and set α1 , . . . , αd > 0 and ᾱ = α1 + · · · + αd . For any ρ > − min(α1 , . . . , αd ), let also
c(α, ρ) = (c(α1 , ρ), . . . , c(αd , ρ)), where c(α, ρ) is as in Definition 1. Then for any − min(α1 , . . . , αd ) < ρ < ∞, ρ , 0,
the angular density of the scaled extremal Dirichlet model with parameters ρ > 0 and α is given, for all w ∈ Sd , by
 d
−ρ−ᾱ d
Y
X n
o1/ρ 
Γ(
ᾱ
+
ρ)


hD (w; ρ, α) =
c(α
,
ρ)w
{c(αi , ρ)}αi /ρ wiαi /ρ−1 .
Q
j
j


d
d−1
d|ρ|
i=1 Γ(αi ) j=1
i=1
The angular density of the positive scaled extremal Dirichlet model with parameters ρ > 0 and α is given, for all
w ∈ Sd , by hpD (w; ρ, α) = hD (w; ρ, α), while the angular density of the negative scaled extremal Dirichlet model with
parameters 0 < ρ < min(α1 , . . . , αd ) and α is given, for all w ∈ Rd+ , by hnD (w; ρ, α) = hD (w; −ρ, α).
From Proposition 5, it is easily seen that when α = 1d , the angular density hpD reduces, for any ρ > 0 and w ∈ Sd ,
to the angular density of the symmetric negative logistic model; see, e.g., Section 4.2 in [7]. In general, the angular
density hpD is not symmetric unless α = α1d .
The positive scaled Dirichlet model can thus be viewed as a new asymmetric generalization of the negative logistic
model which does not place any mass on the vertices or facets of Sd , unless at independence or comonotonicity, i.e.,
when ρ → ∞ and ρ → 0, respectively. Furthermore, hpD can also be interpreted as a generalization of the Coles–
Tawn extremal Dirichlet model. Indeed, hpD (x; 1, α) is precisely the angular density of the latter model given, e.g., in
Equation (3.6) in [7]. Similarly, the negative scaled extremal Dirichlet model is a new asymmetric generalization of
Gumbel’s logistic model [16]. Indeed, when α = 1d , hnD simplifies to the logistic angular density, given, e.g., on p.
381 in [7].

9

Extremal attractors of Liouville copulas
scaled Dirichlet,
α = (1, 1/2, 1/5), ρ = 1/5

scaled Dirichlet,
α = (5/4, 2, 1), ρ = −2/5

w3 = 1

w3 = 1

w1 = 1

w1 = 1

w2 = 1

w1 = 1

w2 = 1

scaled Dirichlet,
α = (1/5, 1/5, 1/5), ρ = 1/5

scaled Dirichlet,
α = (5/4, 5/4, 5/4), ρ = −2/5

w3 = 1

w3 = 1

w2 = 1

w1 = 1

w2 = 1

Figure 2: Angular density of the scaled extremal Dirichlet model with α = (1, 1/2, 1/5), ρ = 1/5 (top left) and
α = (1/5, 1/5, 1/5), ρ = 1/5 (bottom left), α = (5/4, 2, 1), ρ = −2/5 (top right) and α = (5/4, 5/4, 5/4), ρ = −2/5
(bottom right). The colors correspond to log density values and range from red (high density) to blue (low density).

Figures 1 and 2 illustrate the various shapes of hpD and hnD that obtain through various choices of α and ρ. The
asymmetry when α , α1d is clearly apparent. For the same value of ρ, the shapes of the angular density can be quite
different depending on α. In view of the aforementioned closure of both the positive and negative scaled extremal
Dirichlet models under marginalization, this means that these models are able to capture strong dependence in some
pairs of variables (represented by a mode close to 1/2 of the angular density) and at the same time weak dependence
in others pairs (represented by a bathtub shape).
5.2. The bivariate case
When d = 2, the stable tail dependence functions of the positive and negative scaled extremal Dirichlet models have
a closed-form expression in terms of the incomplete beta function given, for any t ∈ (0, 1) and α1 , α2 > 0, by
B(t; α1 , α2 ) =

Z
0

t

xα1 −1 (1 − x)α2 −1 dx.

When t = 1, this integral is the beta function, viz. B(α1 , α2 ) = Γ(α1 )Γ(α2 )/Γ(α1 + α2 ). A direct calculation yields the
corresponding Pickands dependence function, for any t ∈ [0, 1], ApD (t; ρ, α1 , α2 ) = `pD (1 − t, t; ρ, α1 , α2 ), i.e.,

10

L. R. Belzile and J. G. Nešlehová

scaled Dirichlet, ρ < 0

1.0

1.0

0.9

0.9

0.8

0.8

A(t)

A(t)

scaled Dirichlet, ρ > 0

0.7

0.7

0.6

0.6

0.5

0.5
0.0

0.2

0.4

0.6

0.8

1.0

0.0

t

0.2

0.4

0.6

0.8

1.0

t

Figure 3: Pickands dependence function of the scaled extremal Dirichlet model. Left panel: ρ = 4/5 and α = (2, 1/2)
(black full), ρ = 1/4 and α = (1/10, 1/10) (red dashed), ρ = 1/4 and α = (1/2, 1/2) (blue dotted). Right panel:
ρ = −1/4 and α = (2, 1/2) (black full), α = (2/5, 2/5) (red dashed) and α = (1/2, 1/2) (blue dotted).

ApD (t; ρ, α1 , α2 ) =

"
#
(1 − t)
{c(α2 , ρ)(1 − t)}1/ρ
B
;
α
,
α
+
ρ
2
1
B(α2 , α1 + ρ) {c(α2 , ρ)(1 − t)}1/ρ + {c(α1 , ρ)t}1/ρ
"
#
t
{c(α1 , ρ)t}1/ρ
B
+
; α1 , α2 + ρ .
B(α1 , α2 + ρ) {c(α2 , ρ)(1 − t)}1/ρ + {c(α1 , ρ)t}1/ρ

When α1 = α2 = 1, ApD becomes the Pickands dependence function of the Galambos copula, viz. ApD (t; ρ, 1, 1) =

1 − t−1/ρ + (1 − t)−1/ρ −ρ , as expected given that the positive scaled extremal Dirichlet model becomes the symmetric
negative logistic model in this case.
Similarly, for any t ∈ [0, 1], the Pickands dependence function AnD (t; ρ, α1 , α2 ) = `nD (1 − t, t; ρ, α1 , α2 ) equals
AnD (t; ρ, α1 , α2 ) =

"
#
(1 − t)
{(1 − t)c(α2 , −ρ)}1/ρ
B
;
α
−
ρ,
α
1
2
B(α1 − ρ, α2 ) {tc(α1 , −ρ)}1/ρ + {(1 − t)c(α2 , −ρ)}1/ρ
#
"
t
{tc(α1 , −ρ)}1/ρ
; α2 − ρ, α1 .
+
B
B(α2 − ρ, α1 ) {tc(α1 , −ρ)}1/ρ + {(1 − t)c(α2 , −ρ)}1/ρ

When α1 = α2 = 1, AnD simplifies to the stable tail dependence function of the Gumbel extreme-value copula, viz.

AnD (t; ρ, 1, 1) = t1/ρ + (1 − t)1/ρ ρ . This again confirms that the negative scaled extremal Dirichlet model becomes
the symmetric logistic model when α1 = α2 = 1. The Pickands dependence functions ApD and AnD are illustrated in
Figure 3, for the same choices of parameters and the corresponding angular density shown in Figure 1.
The above formulas for ApD and AnD now easily lead to expressions for their upper tail dependence coefficients.
Recall that for an arbitrary bivariate copula C, the lower and upper tail dependence coefficients of [25] are given by
C(u, u)
,
u→0
u

λ` (C) = lim

C(u, u) − 1
Ĉ(u, u)
= lim
,
u→1
u→0
u−1
u

λu (C) = 2 − lim

where Ĉ is the survival copula of C, provided these limits exist. When C is bivariate extreme-value with Pickands
dependence function A, it follows easily from (4) that λ` (C) = 0 and λu (C) = 2 − 2A(1/2).

Extremal attractors of Liouville copulas

11

pD

Now suppose that Cρ,α is a bivariate extreme-value copula with positive scaled extremal Dirichlet Pickands dependence function ApD and parameters ρ > 0 and α1 , α2 > 0. Then
pD
λu (Cρ,α )

(
)
c(α2 , ρ)1/ρ
1
B
; α2 , α1 + ρ
=2−
B(α2 , α1 + ρ)
c(α2 , ρ)1/ρ + c(α1 , ρ)1/ρ
(
)
1
c(α1 , ρ)1/ρ
−
B
;
α
,
α
+
ρ
. (8)
1
2
B(α1 , α2 + ρ)
c(α2 , ρ)1/ρ + c(α1 , ρ)1/ρ

nD
Similarly, if Cρ,α
is a bivariate extreme-value copula with negative scaled extremal Dirichlet Pickands dependence
nD
function A and parameters α1 , α2 > 0 and 0 < ρ < min(α1 , α2 ), then

nD
λu (Cρ,α
)=2−

(
)
1
c(α2 , −ρ)1/ρ
B
;
α
−
ρ,
α
1
2
B(α1 − ρ, α2 )
c(α1 , −ρ)1/ρ + c(α2 , −ρ)1/ρ
)
(
1
c(α1 , −ρ)1/ρ
;
α
−
ρ,
α
+
B
2
1 . (9)
B(α2 − ρ, α1 )
c(α1 , −ρ)1/ρ + c(α2 , −ρ)1/ρ

In the symmetric case α1 = α2 ≡ α, Expressions (8) and (9) simplify to
pD
λu (Cρ,α )

!
1
2
=2−
B ; α, α + ρ ,
B(α, α + ρ) 2

nD
λu (Cρ,α
)

!
2
1
=2−
B ; α − ρ, α .
B(α − ρ, α) 2

Formulas (8) and (9) lead directly to expressions for the tail dependence coefficients of Liouville copulas. This is
because if C ∈ M(C0 ), where C0 is an extreme-value copula with Pickands tail dependence function A0 , λu (C) =
2 − 2A0 (1/2) [29, Proposition 7.51]. Similarly, if Ĉ ∈ M(C0∗ ), where C0∗ is an extreme-value copula with Pickands
tail dependence function A∗0 , λ` (C) = 2 − 2A∗0 (1/2). The following corollary is thus an immediate consequence of
Corollaries 1 and 2.
Corollary 3 Suppose that C is the survival copula of a Liouville random vector RDα with parameters α > 0 and a
radial part R such that Pr(R ≤ 0) = 0. Then the following statements hold.
pD

(a) If R ∈ M(Φρ ) for some ρ > 0, λ` (C) = λu (Cρ,α ) is given by Eq. (8).
(b) If R ∈ M(Λ) or R ∈ M(Ψρ ) for some ρ > 0, λ` (C) = 0.
nD
(c) If 1/R ∈ M(Φρ ) for some 0 < ρ < α1 ∧ α2 , λu (C) = λu (Cρ,α
) is given by Eq. (9).

(d) If 1/R ∈ M(Φρ ) for ρ > α1 ∧ α2 or if E(1/Rβ ) < ∞ for β > α1 ∨ α2 , λu (C) = 0.
The role of the parameters α and ρ is best explained if we consider the reparametrization ∆α = |α1 − α2 | and
Σα = α1 + α2 . As is the case for the Dirichlet distribution, the level of dependence is higher for large values of Σα .
Furthermore, λu is monotonically decreasing in ρ. Higher levels of extremal asymmetry, as measured by departures
from the diagonal on the copula scale, are governed by both Σα and ∆α . The larger Σα , the lower the asymmetry.
Likewise, the larger ∆α , the larger the asymmetry. Contrary to the case of extremal dependence, the behavior in ρ is
not monotone. For the negative scaled extremal Dirichlet model, asymmetry is maximal when ρ ≈ α1 ∧ α2 . When
Σα is small, smaller values of ρ induce larger asymmetry, but this is not the case for larger values of Σα where the
asymmetry profile is convex with a global maximum attained for larger values of ρ.

6. de Haan representation and simulation algorithms
Random samples from the scaled extremal Dirichlet model can be drawn efficiently using the algorithms recently
developed in [10]. We first derive the so-called de Haan representation in Section 6.1 and adapt the algorithms from
[10] to the present setting in Section 6.2.

12

L. R. Belzile and J. G. Nešlehová

6.1. de Haan representation
First, introduce the following family of univariate distributions, which we term the scaled Gamma family and denote
by sGa(a, b, c). It has three parameters a, c > 0 and b , 0 and a density given, for all x > 0, by
(  b )
|b| −bc bc−1
x
f (x; a, b, c) =
a x
exp −
.
(10)
Γ(c)
a
d

Observe that when Z ∼ Ga(c, 1) is a Gamma variable with shape parameter c > 0 and scaling parameter 1, Y = aZ 1/b is
scaled Gamma sGa(a, b, c). Consequently, E (Y) = aΓ(c + 1/b)/Γ(c) < ∞ provided that b < −1/c. The scaled Gamma
family includes several well-known distributions as special cases, notably the Gamma when b = 1, the Weibull when
c = 1 and b > 0, the inverse Gamma when b = −1, and the Fréchet when c = 1 and b < 0. When b > 0, the scaled
Gamma is the generalized Gamma distribution of [41], albeit in a different parametrization.
Now consider the parameters α = (α1 , . . . , αd ) with α > 0d and ρ > − min(α1 , . . . , αd ), ρ , 0. Let V be a random
vector with independent scaled Gamma margins Vi ∼ sGa{1/c(αi , ρ), 1/ρ, αi }, where for α > 0, c(α, ρ) = Γ(α+ρ)/Γ(α)
as in Definition 1. If Z is a random vector with independent Gamma margins Zi ∼ Ga(αi , 1) then for all i = 1, . . . , d,
d
Vi = Ziρ /c(αi , ρ). Furthermore, recall that kZk ∼ Ga(ᾱ, 1) is independent of Z/kZk, which has the same distribution
as the Dirichlet vector Dα = (D1 , . . . , Dd ). One thus has, for all x ∈ Rd+ ,
)#
)#
"
(
"
(


xi Dρi
xi Ziρ
ρ
= E (kZk ) E max
= `D (x; ρ, α),
(11)
E max (xi Vi ) = E max
1≤i≤d c(αi , ρ)
1≤i≤d
1≤i≤d c(αi , ρ)
where `D is as in Definition 1, given that E (kZkρ ) = c(ᾱ, ρ).
When ρ = 1, the positive scaled Dirichlet extremal model becomes the Coles–Tawn Dirichlet extremal model,
Vi ∼ Ga(αi , 1) and Eq. (11) reduces to the representation derived in [37]. When α = 1d , `D becomes the stable tail
dependence function of the negative logistic model, Vi is Weibull and Eq. (11) is the representation in Appendix A.2.4
of [10]. Similarly, when ρ < 0 and α = 1d , the negative scaled Dirichlet extremal model becomes the logistic model,
Vi is Fréchet and Eq. (11) is the representation in Appendix A.2.4 of [10]. The requirement that ρ > − min(α1 , . . . , αd )
ensures that the expectation of Vi is finite for all i ∈ {1, . . . , d}.
Eq. (11) implies that the max-stable random vector Y with unit Fréchet margins and extreme-value copula with
stable tail dependence function `D (·; ρ, α) admits the de Haan [9] spectral representation
d

Y = max ζk V k ,
k∈N

(12)

−2
where Z = {ζk }∞
dζ and V k is an i.i.d. sequence of rank=1 is a Poisson point process on (0, ∞) with intensity ζ
dom vectors independent of Z. Furthermore, the univariate margins of V k are independent and such that Vk j ∼
sGa{1/c(α j , ρ), 1/ρ, α j } for j = 1, . . . , d with E (V k ) = 1d for all k ∈ N.

6.2. Unconditional simulation
The de Haan representation (12) offers, among other things, an easy route to unconditional simulation of max-stable
random vectors that follow the scaled Dirichlet extremal model, as laid out in [10] in the more general context of maxstable processes. To see how this work applies in the present setting, fix an arbitrary j0 ∈ {1, . . . , d} and recall that the
j0 th extremal function φ+j0 is given, almost surely, as ζk V k such that Y j0 = ζk Vk j0 . From Eq. (12) and Proposition 1 in
d

[10] it then directly follows that φ+j0 /Y j0 = (W j0 1 /W j0 j0 , . . . , W j0 d /W j0 j0 ), where W j0 = (W j0 1 , . . . , W j0 d ) is a random
vector with density given, for all x ∈ Rd+ , by
d
Y
h
i
h
i
|1/ρ|
|1/ρ|
α j /ρ
α /ρ−1
c(α j0 , ρ)α j0 /ρ x j0 0 exp −{c(α j0 , ρ)x j0 }1/ρ ×
c(α j , ρ)α j /ρ x j j
exp −{c(α j , ρ)x j }1/ρ .
Γ(α j0 )
Γ(α j )
j=1, j, j
0

This means that the components of W j0 are independent and such that W j0 j ∼ sGa{1/c(α j , ρ), 1/ρ, α j } when j , j0
and W j0 j0 ∼ sGa{1/c(α j0 , ρ), 1/ρ, α j0 + ρ}. In other words, W j0 j0 ∼ Z ρj0 /c(α j0 , ρ) where Z j0 ∼ Ga(α j0 + ρ, 1), while for
all j , j0 , W j0 j = Z ρj /c(α j , ρ) where Z j ∼ Ga(α j , 1).
d

Extremal attractors of Liouville copulas

13

The exact distribution of φ+j0 /Y j0 given above now allows for an easy adaptation of the algorithms in [10]. To
draw an observation from the extreme-value copula with the scaled Dirichlet stable tail dependence function `D with
parameters α > 0d and ρ > − min(α1 , . . . , αd ), ρ , 0, one can follow Algorithms 1 and 2 below. The first procedure
corresponds to Algorithm 1 in [10] and relies on [36]; the second is an adaptation of Algorithm 2 in [10].
Algorithm 1 Exact simulations from the extreme-value copula based on spectral densities.
1: Simulate E ∼ Exp(1) .
2: Set Y = 0.
3: while 1/E > min(Y1 , . . . , Yd ) do
4:
Simulate J from the uniform distribution on {1, . . . , d}.
5:
Simulate independent Z j ∼ Ga(α j , 1) for j ∈ {1, . . . , d} \ J and Z J ∼ Ga(α J + ρ, 1).
6:
Set W j ← Z ρj /c(α j , ρ), j = 1, . . . , d.
7:
Set S ← W/kWk.
8:
Update Y ← max{Y, dS/E}.
9:
Simulate E ∗ ∼ Exp(1) and update E ← E + E ∗ .
10: return U = exp(−1/Y).

Algorithm 2 Exact simulations based on sequential sampling of the extremal functions.
Simulate Z1 ∼ Ga(α1 + ρ, 1) and Z j ∼ Ga(α j , 1), j = 2, . . . , d.
Compute W where W j ← Z ρj /c(α j , ρ), j = 1, . . . , d.
Simulate E1 ∼ Exp(1).
Set Y ← W/(W1 E1 ).
for k = 2, . . . , d do
Simulate Ek ∼ Exp(1).
while 1/Ek > Yk do
Simulate independent Zk ∼ Ga(αk + ρ, 1) and Z j ∼ Ga(α j , 1), j = 1, . . . , d, j , k .
Set W = (W1 , . . . , Wd ) where W j ← Z ρj /c(α j , ρ), j = 1, . . . , d.
if Wi /(Wk Ek ) < Yi for all i = 1, . . . , k − 1 then
Update Y ← max{Y, W/(Wk Ek )}.
12:
Simulate E ∗ ∼ Exp(1) and update Ek ← Ek + E ∗ .
13: return U = exp(−1/Y).
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:

Note that S obtained in Step 7 of Algorithm 1 has the angular distribution σd of `D ; see Theorem 1 in [10].
Similar algorithms for drawing samples from the angular distribution of the extremal logistic and Dirichlet models
were obtained in [3]. Algorithm 2 requires a lower number of simulations and is more efficient on average, cf. [10].
Both algorithms are easily implemented using the function rmev in the mev package within the R Project for Statistical
Computing [34], which returns samples of max-stable scaled extremal Dirichlet vectors with unit Fréchet margins,
i.e., Y in Algorithms 1 and 2.

7. Estimation
The scaled extremal Dirichlet model can be used to model dependence between extreme events. To this end, several
schemes can be envisaged. For example, one can consider the block-maxima approach, given that max-stable distributions are the most natural for such data. Another option is peaks-over-threshold models. Yet another alternative,
used in [13] for the Brown–Resnick model, is to approximate the conditional distribution of a random vector with unit
Fréchet margins given that the j0 th component exceeds a large threshold by the distribution of φ+j0 /Y j0 discussed in
Section 6.2.
Here, we focus on the multivariate tail model of [28]; see also Section 16.4 in [29]. To this end, let X1 , . . . , Xn be a
random sample from some unknown multivariate distribution H with continuous univariate margins which is assumed
to be in the maximum domain of attraction of a multivariate extreme-value distribution H0 . To model the tail of H, its
margins F j , j = 1, . . . , d can first be approximated using the univariate peaks-over-threshold method. For all x above

14

L. R. Belzile and J. G. Nešlehová

some high threshold u j , one then has
F j (x) ≈ F̃ j (x; η j , ξ j ) = 1 − ν j 1 + ξ j

(x − u j )
ηj

!−1/ξ j
,

(13)

+

where ν j = 1 − F j (u j ), and η j > 0 and ξ j are the parameters of the generalized Pareto distribution. Furthermore, for
w sufficiently close to 1d , the copula of H can be approximated by the extreme-value copula C0 of H0 , so that, for
x ≥ u, H(x) ≈ H̃(x) = C0 {F̃1 (x1 ), . . . , F̃d (xd )}. The parameters of this multivariate tail model, i.e., the parameters
θ of the stable tail dependence function `0 of C0 as well as the marginal parameters ν, η and ξ can be estimated
using likelihood methods; this allows, e.g., for Bayesian inference, generalized additive modeling of the parameters
and model selection based on likelihood-ratio tests. For a comprehensive review of likelihood inference methods for
extremes, see, e.g., [22].
The multivariate tail model can be fitted in low-dimensions using the censored likelihood L(X; ν, η, ξ, θ) =
Qn
i=1 Li (X i ; ν, η, ξ, θ), where for i = 1, . . . , n,
Li (Xi ; u, ν, η, ξ, θ) =

∂mi H̃(y1 , . . . , yd )
∂y j1 · · · ∂y jmi

=
y=max(Xi ,u)

∂mi exp{−`0 (1/y)}
∂y j1 · · · ∂y jmi

mi
Y

J jk (Xi jk )

(14)

y=t{max(Xi ,u)} k=1

In this expression, the indices j1 , . . . , jmi are those of the components of Xi exceeding the thresholds u and for x ≥ u,
t(x) = (t1 (x1 ), . . . , td (xd )), where for j = 1, . . . , d,
(x j − u j )
νj
1 + ξj
J j (x j ) =
ηj
ηj

1
,
t j (x j ) = −
log{F̃ j (x j ; η j , ξ j )}

!−1/ξ j −1

1
.
[log{F̃ j (x j ; η j , ξ j )}]2 F̃ j (x j ; η j , ξ j )

(15)

The censored likelihood L(X; ν, η, ξ, θ) can be maximized either over all parameters at once, or the marginal parameters ν, η and ξ can be estimated from each univariate margin separately, so that only the estimate of θ is obtained
through maximizing L. When d is large, one can also maximize the likelihood in [40] that uses the tail approximation
H̄(x) ≈ 1 − `(1/x). In either case, `0 and the higher-order partial derivatives of `0 (1/x) need to be computed.
When `0 is the scaled extremal Dirichlet stable tail dependence function `D (·; ρ, α) given in Definition 1 with
parameters α > 0 and ρ > − min(α1 , . . . , αd ), ρ , 0, its expression is not explicit. However, `D can be calculated numerically using adaptive numerical cubature algorithms for integrals of functions defined on the simplex, as
implemented in, e.g., the R package SimplicialCubature. Given the representation in Eq. (5), `D is also easily
approximated using Monte Carlo methods. Instead of employing Eq. (5) directly and sampling from the Dirichlet
vector Dα , one can use the more efficient importance sampling estimator
h
i
−1 ρ
B max
X
1≤ j≤d {c(α j , ρ)u j } Di j
1
`bD (1/u, ρ, α) =
,
1 Pd
−1 ρ
B i=1
j=1 c(α j , ρ) Di j
d
P
where Di ∼ d−1 dj=1 Dir(α + I j ρ1d ) is sampled from a Dirichlet mixture.
The partial derivatives of `D can be calculated using the following result, shown in E.
Proposition 6 Let `D be the scaled extremal Dirichlet stable tail dependence function with parameters α > 0d and
− min(α1 , . . . , αd ) < ρ < ∞, ρ , 0. Then, for any x ∈ Rd+ ,
 d
−ρ−ᾱ d
Y
X n
o1/ρ 
∂d `D (1/x)
Γ(ᾱ + ρ)
D


= −dh (x; ρ, α) = −
c(α j , ρ)x j 
{c(αi , ρ)}αi /ρ xiαi /ρ−1 ,
Qd
d−1
∂x1 · · · ∂xd
|ρ|
Γ(α
)
i
i=1
j=1
i=1
where hD is as given in Proposition 5. Furthermore, for all k = 1, . . . , d − 1 and x ∈ Rd+ ,
∂k `D (1/x)
=−
∂x1 · · · ∂xk

Z

∞

tk
0

k
Y
i=1

f xi t;

!Y
!
d
1
1
1
1
, , αi
F xi t;
, , αi dt,
c(αi , ρ) ρ
c(αi , ρ) ρ
i=k+1

(16)

15

Extremal attractors of Liouville copulas

where f (; a, b, c) and F(; a, b, c) denote, respectively the density and distribution function of the Rscaled Gamma
x
distribution with parameters a, c > 0 and b , 0 given in Eq. (10). Furthermore, if γ(c, x) = 0 tc−1 e−t dt denotes the lower incomplete gamma function, then for x > 0, F(x; a, b, c) = γ{c, (x/a)b }/Γ(c) when b > 0 while
F(x; a, b, c) = 1 − γ{c, (x/a)b }/Γ(c) when b < 0.
Other estimating equations could be used to circumvent the calculation of `D (1/x) and its partial derivatives. An
interesting alternative to likelihoods in the context of proper scoring functions is proposed in [8]. Specifically, the
authors advocate the use of the gradient score, adapted by them for the peaks-over-threshold framework,


(
) 
d 
X
 ∂2 log h(x) 1 ∂ log h(x) 2 
∂wi (x) ∂ log h(x)


2

+
+ wi (x) 
δw (x) =
2wi (x)
∂xi
∂xi
2
∂xi
∂xi2
i=1
for a differentiable weighting function w(x), unit Fréchet observations x and density h(x) that would correspond in
the setting of the scaled extremal Dirichlet to dhD (x; ρ, α). Explicit expressions for the derivatives of log dhD may be
P
found in E. The parameter estimates are obtained as the solution to argmaxθ∈Θ ni=1 δw (xi )IR(xi /u)>1 , where θ = (ρ, α)
is the vector of parameters of the model and R is a differentiable risk functional, usually the ` p norm for some p ∈ N.
Although the gradient score is not asymptotically most efficient, weighting functions can be designed to reproduce
approximate censoring, lending the method robustness and tractability.

8. Data illustration
In this section, we illustrate the use of the scaled extremal Dirichlet model on a trivariate sample of daily river
flow data of the river Isar in southern Germany; this dataset is a subset of the one analyzed in [1]. All the code
can be downloaded from https://github.com/lbelzile/ealc. For this analysis, we selected data measured at
Lenggries (upstream), Pupplinger Au (in the middle) and Munich (downstream). To ensure stationarity of the series
and given that the most extreme events occur during the summer, we restricted our attention to the measurements
for the months of June, July and August. Since the sites are measuring the flow of the same river, dependence at
extreme levels is likely to be present, as is indeed apparent from Figure 4. Directionality of the river may further lead
to asymmetry in the asymptotic dependence structure, suggesting that the scaled extremal Dirichlet model may be
well suited for these data. Furthermore, given that other well-known models like the extremal Dirichlet, logistic and
negative logistic are nested within this family, their adequacy can be assessed through likelihood ratio tests.
To remove dependence between extremes over time, we decluster each series and retain only the cluster maxima
based on three-day runs. Rounding of the measurements has no impact on parameter estimates and is henceforth
neglected. The multivariate tail model outlined in Section 7 is next fitted to the cluster maxima. The thresholds
u = (u1 , u2 , u3 ) were selected to be the 92% quantiles using the parameter stability plot of [42] (not shown here).
Next, set θ = (η, ξ, α, ρ), where η and ξ are the marginal parameters of the generalized Pareto distribution in Eq. (13)
and ρ and α are the parameters of the scaled Dirichlet model. To estimate θ, the trivariate censored likelihood (14)
could be used. To avoid numerical integration and because of the relative robustness to misspecification, we employed
the pairwise composite log-likelihood lC of [28] instead; the loss of efficiency in this trivariate example is likely small.
Specifically, we maximized
lC (θ) =

n X
d−1 X
d h
X

i
log g{t j (xi j ), tk (xik ); θ, t j (u j ), tk (uk )} + I xi j >u j log J j (xi j ) + I xik >uk log Jk (xik ) ,

i=1 j=1 k= j+1

where



exp{−`(1/u j , 1/uk )},






−∂`(1/y j , 1/uk )/∂y j exp{−`(1/y j , 1/uk )},
g(y j , yk ; θ, u j , uk ) = 


−∂`(1/u j , 1/yk )/∂yk exp{−`(1/u j , 1/yk )},


on
o
hn
i



 ∂`(1/y j , 1/yk )/∂y j ∂`(1/y j , 1/yk )/∂yk − dhD (y j , yk ) exp{−`(1/y j , 1/yk )},
where ` = `D and for all j = 1, . . . , d, t j and J j are as in Eq. (15).

yj
yj
yj
yj

≤ u j , yk
> u j , yk
≤ u j , yk
> u j , yk

≤ uk
≤ uk
> uk
> uk

16

L. R. Belzile and J. G. Nešlehová
Isar daily river flow
800

1958-11-01 / 2014-01-18

Muenchen

800

600

600

400

400

200

200

400 Pupplinger Au

400

300

300

200

200

100

100

400

400

Lenggries
300

300

200

200

100

100

1958

1962

1964

1966

1968

1970

1972

1974

1976

1978

1980

1982

1984

1986

1988

1990

1992

1994

1996

1998

2000

2002

2004

2006

2008

2010

2012

2014

Figure 4: Daily river flow of the Isar river at three sites
Uncertainty assessment can be done in the same way as for general estimating equations. Specifically, let g(θ)
denote an unbiased estimating function and define the variability matrix J, the sensitivity matrix H and the Godambe
information matrix G as
!
!
∂g(θ) ∂g(θ) >
∂2 g(θ)
J=E
,
H = −E
,
G = HJ−1 H.
(17)
∂θ ∂θ
∂θ∂θ>
The maximum composite likelihood estimator is strongly consistent and asymptotically normal, centered at the true
parameter θ with covariance matrix given by the inverse Godambe matrix G−1 .
Using the pairwise composite log-likelihood lC , we fitted the scaled extremal Dirichlet model as well as the
logistic and negative logistic models that correspond to the negative and positive scaled extremal Dirichlet models,
respectively, and the parameter restriction α = 1d . The estimates of the marginal generalized Pareto parameters η and
ξ are given in Table 1. As the estimates were obtained by maximizing `C , their values depend on the fitted model; the
line labeled “Marginal" corresponds to fitting the generalized Pareto distribution to threshold exceedances of each one
of the three series separately. The marginal Q-Q plots displayed in Figure 5 indicate a good fit of the model as well.

Scaled Dirichlet
Neg. logistic
Logistic
ext. Dirichlet
Marginal

η1

η2

η3

ξ1

ξ2

ξ3

123.2 (7.5)
117.1 (6.8)
117.3 (6.8)
114.4 (6.8)
129.1 (14.5)

84.4 (5)
86.2 (5.1)
86.6 (5.1)
84.3 (4.9)
95.1 (10.6)

68.1 (4.2)
70 (4.3)
70.4 (4.3)
68.2 (4.1)
76 (8.7)

0.05 (0.04)
0.08 (0.04)
0.08 (0.04)
0.12 (0.04)
-0.01 (0.08)

-0.03 (0.04)
-0.05 (0.04)
-0.05 (0.04)
-0.02 (0.04)
-0.15 (0.08)

0.02 (0.04)
0 (0.04)
0 (0.04)
0.04 (0.04)
-0.08 (0.08)

Table 1: Generalized Pareto parameter estimates and standard errors (in parenthesis) for the trivariate river example
for four different models.

17

Extremal attractors of Liouville copulas

600

700

500
400
300

500
100

200

100

200

400

200

300

400

Sample quantiles

Sample quantiles

600

1200
1000
800
600

Sample quantiles

Q-Q plot, Lenggries
700

Q-Q plot, Pupplinger Au

1400

Q-Q plot, Muenchen

200 300 400 500 600 700 800 900

200

Theoretical quantiles

300

400

500

100

200

Theoretical quantiles

300

400

Theoretical quantiles

Figure 5: Marginal Q-Q plots for the three sites based on pairwise composite likelihood estimates for the scale and
shape parameters obtained from the scaled Dirichlet model, retaining marginal exceedances of the 92% quantiles. The
pointwise confidence intervals were obtained from the transformed Beta quantiles of the order statistics.

Scaled Dirichlet
Neg. logistic
Logistic
ext. Dirichlet
Gradient score

α1

α2

α3

ρ

0.76 (0.3)
1
1
3.34 (0.52)
1

1.65 (0.82)
1
1
10.2 (2.84)
2.72

2.03 (1.15)
1
1
12.78 (3.93)
2.66

−0.32 (0.1)
0.36 (0.02)
0.28 (0.01)
1
−0.39

Table 2: Dependence parameters estimates and standard errors (in parenthesis) for the trivariate river example.

w1 = 1

Negative logistic

Logistic

scaled
Dirichlet

w3 = 1

w3 = 1

w3 = 1

w2 = 1

w1 = 1

w2 = 1

w1 = 1

w2 = 1

Figure 6: Angular density plots for the three models, the negative logistic (left), logistic (middle) and scaled Dirichlet
(right). The colours correspond to log density values and range from red (high density) to blue (low density).
The estimates of the dependence parameters α and ρ are given in Table 2. The last line displays the maximum gradient score estimates were obtained from the raw data, i.e., ignoring the clustering, after transforming the observations
to the standard Fréchet scale using the probability integral transform. We retained only the 10% largest values based
on the ` p norm with p = 20; this risk functional is essentially a differentiable approximation of `∞ . We selected the
weight function w(x, u) = x[1 − exp{−(kxk p /u − 1)}] based on [8] to reproduce approximate censoring. The estimates
are similar to the composite maximum likelihood estimators, though not efficient.

18

L. R. Belzile and J. G. Nešlehová

The angular densities of the fitted logistic, negative logistic and scaled extremal Dirichlet models are displayed
in Figure 6. The right panel of this figure shows asymmetry caused by a few extreme events that only happened
downstream. Whether this asymmetry is significant can be assessed through composite likelihood ratio tests; recall
that the logistic model, the negative logistic model and the extremal Dirichlet model of [7] are all nested within
the scaled extremal Dirichlet model. To this end, consider a partition of θ = (ψ, λ) into a q dimensional parameter of
interest ψ and a 3d +1−q dimensional nuisance parameter λ, and the corresponding partitions of the matrices H, J and
G. Let b
θC = (b
ψC , b
λC ) denote the maximum composite likelihood parameter estimates and b
θ0 = (ψ0 , b
λ0 ) the restricted
parameter estimates under the null hypothesis that the simpler model is adequate. The asymptotic distribution of
P
the composite likelihood ratio test statistic 2{log lC (θbC ) − log lC (θb0 )} is equal to qi=1 ci Zi where Zi are independent
−1
χ21 variables and ci are the eigenvalues of the q × q matrix (Hψψ − Hψλ H−1
λλ Hλψ )Gψψ ; see [26]. We estimated the
−1
inverse Godambe information matrix, G , by the empirical covariance of B nonparametric bootstrap replicates. The
sensitivity matrix H was obtained from the Hessian matrix at the maximum composite likelihood estimate and the
variability matrix J from Eq. (17). Since the Coles–Tawn extremal Dirichlet, negative logistic and logistic models are
nested within the scaled Dirichlet family, we test for a restriction to these simpler models; the respective approximate
P-values were 0.003, 0.74 and 0.78. These values suggest that while the Coles–Tawn extremal Dirichlet model is
clearly not suitable, there is not sufficient evidence to discard the logistic and negative logistic models. The effects of
possible model misspecification are also visible for the Coles–Tawn extremal Dirichlet model, as the parameter values
of α1 , α2 and α3 are very large (viz. Table 2) and this induces negative bias in the shape parameter estimates, as can
be seen from Table 1.

9. Discussion
In this article, we have identified extremal attractors of copulas and survival copulas of Liouville random vectors
RDα , where Dα has a Dirichlet distribution on the unit simplex with parameters α, and R is a strictly positive random
variable independent of Dα . The limiting stable tail dependence functions can be embedded in a single family,
which can capture asymmetry and provides a valid model in dimension d. The latter is novel and termed here the
scaled extremal Dirichlet; it includes the well-known logistic, negative logistic as well as the Coles–Tawn extremal
Dirichlet models as special cases. In particular, therefore, this paper is first to provide an example of a random vector
attracted to the Coles–Tawn extremal Dirichlet model, which was derived by enforcing moment constraints on a
simplex distribution rather than as the limiting distribution of a random vector.
A scaled extremal Dirichlet stable tail dependence function `D has d + 1 parameters, ρ and α. The parameter
vector α is inherited from Dα and induces asymmetry in `D . The parameter ρ comes from the regular variation of
R at zero and infinity, respectively; this is reminiscent of the extremal attractors of elliptical distributions [33]. The
magnitude of ρ has impact on the strength of dependence while its sign changes the overall shape of `D . Having
d + 1 parameters, the scaled extremal Dirichlet model may not be sufficiently rich to account for spatial dependence,
unlike the Hüsler–Reiss or the extremal Student-t models, which have one parameter for each pair of variables and
are thus easily combined with distances. Also, it is less flexible than Dirichlet mixtures [4], which are however hard
to estimate in high dimensions and require sophisticated machinery. To achieve greater flexibility, the scaled extremal
Dirichlet model could perhaps be extended by working with more general scale mixtures, such as of the weighted
Dirichlet distributions considered, e.g., in [18].
Nonetheless, the scaled extremal Dirichlet model may naturally find applications whenever asymmetric extremal
dependence is suspected; the latter may be caused, e.g., by causal relationships between the variables [15]. The
stochastic structure of the scaled extremal Dirichlet model has several major advantages, that make the model easy
to interpret, estimate and simulate from. Its angular density has a simple form; in contrast to the asymmetric generalizations of the logistic and negative logistic models, this model does not place any mass on the vertices and
lower-dimensional facets of the unit simplex. Another plus is the tractability of the de Haan representation and of
the extremal functions, both expressible in terms of independent scaled Gamma variables; this allows for feasible
inference and stochastic simulation. While the scaled extremal Dirichlet stable tail dependence function `D does not
have a closed form in general, closed-form algebraic expressions exist when α is integer-valued and in the bivariate
case. Model selection for well-known families of extreme-value distributions can be performed through likelihood
ratio tests. Another potentially useful feature is that ρ ∈ (−∞, ∞) can be allowed, with the convention that all variables
whose indices i are such that −ρ ≤ −αi are independent.

Extremal attractors of Liouville copulas

19

Acknowledgment
Funding in partial support of this work was provided by the Natural Sciences and Engineering Research Council
(RGPIN-2015-06801, CGSD3-459751-2014), the Canadian Statistical Sciences Institute, and the Fonds de recherche
du Québec – Nature et technologies (2015–PR–183236). We thank the acting Editor-in-Chief, Richard Lockhart, the
Associate Editor and two anonymous referees for their valuable suggestions.

References
[1] P. Asadi, A. C. Davison, and S. Engelke. Extremes on river networks. Ann. Appl. Stat., 9(4):2023–2050, 2015.
[2] G. Balkema and N. Nolde. Asymptotic independence for unimodal densities. Adv. in Appl. Probab., 42(2):411–
432, 2010.
[3] M.-O. Boldi. A note on the representation of parametric models for multivariate extremes. Extremes, 12(3):211–
218, 2009.
[4] M.-O. Boldi and A. C. Davison. A mixture model for multivariate extremes. J. Roy. Stat. Soc. B Met., 69(2):217–
229, 2007.
[5] L. Breiman. On some limit theorems similar to the arc-sin law. Theory Probab. Appl., 10(2):323–331, 1965.
[6] A. Charpentier and J. Segers. Tails of multivariate Archimedean copulas. J. Multivariate Anal., 100:1521–1537,
2009.
[7] S. G. Coles and J. A. Tawn. Modelling extreme multivariate events. J. Roy. Stat. Soc. B Met., 53(2):377–392,
1991.
[8] R. de Fondeville and A. C. Davison. High-dimensional peaks-over-threshold inference for the Brown-Resnick
process. ArXiv e-prints, 2016.
[9] L. de Haan. A spectral representation for max-stable processes. Ann. Probab., 12(4):1194–1204, 1984.
[10] C. Dombry, S. Engelke, and M. Oesting. Exact simulation of max-stable processes. Biometrika, 103(2):303–317,
2016.
[11] P. Embrechts and C. M. Goldie. On closure and factorization properties of subexponential and related distributions. J. Austral. Math. Soc. Ser. A, 29(2):243–256, 1980.
[12] P. Embrechts, C. Klüppelberg, and T. Mikosch. Modelling Extremal Events for Insurance and Finance. Springer,
New York, 1997.
[13] S. Engelke, A. Malinowski, Z. Kabluchko, and M. Schlather. Estimation of Hüsler–Reiss distributions and
Brown–Resnick processes. J. Roy. Stat. Soc. B, 77(1):239–265, 2015.
[14] K.-T. Fang, S. Kotz, and K. W. Ng. Symmetric Multivariate and Related Distributions. Chapman & Hall,
London, 1990.
[15] C. Genest and J. G. Nešlehová. Assessing and modeling asymmetry in bivariate continuous data. In P. Jaworski,
F. Durante, and W. K. Härdle, editors, Copulae in Mathematical and Quantitative Finance, pages 91–114.
Springer, Berlin, 2013.
[16] É. J. Gumbel. Distributions des valeurs extrêmes en plusieurs dimensions. Publ. Inst. Statist. Univ. Paris,
9:171–173, 1960.
[17] R. D. Gupta and D. S. P. Richards. Multivariate Liouville distributions. J. Multivariate Anal., 23(2):233–256,
1987.
[18] E. Hashorva. Extremes of weighted dirichlet arrays. Extremes, 11(4):393–420, 2008.
[19] E. Hashorva. Extremes of aggregated dirichlet risks. J. Multivariate Anal., 133:334 – 345, 2015.
[20] E. Hashorva and A. G. Pakes. Tail asymptotics under beta random scaling. J. Math. Anal. Appl., 372(2):496–514,
2010.
[21] L. Hua. A note on upper tail behaviour of Liouville copulas. Risks, 4(4):40, 2016.
[22] R. Huser, A. C. Davison, and M. G. Genton. Likelihood estimators for multivariate extremes. Extremes,
19(1):79–103, 2016.
[23] J. Hüsler and R.-D. Reiss. Maxima of normal random vectors: between independence and complete dependence.
Statist. Probab. Lett., 7(4):283–286, 1989.
[24] H. Joe. Families of min-stable multivariate exponential and multivariate extreme value distributions. Statist.
Probab. Lett., 9:75–81, 1990.
[25] H. Joe. Multivariate dependence measures and data analysis. Comput. Statist. Data Anal., 16:279–297, 1993.

20

L. R. Belzile and J. G. Nešlehová

[26] J. T. Kent. Robust properties of likelihood ratio test. Biometrika, 69(1):19–27, 1982.
[27] M. Larsson and J. Nešlehová. Extremal behavior of Archimedean copulas. Adv. Appl. Probab., 43:195–216,
2011.
[28] A. W. Ledford and J. A. Tawn. Statistics for near independence in multivariate extreme values. Biometrika,
83(1):169–187, 1996.
[29] A. J. McNeil, R. Frey, and P. Embrechts. Quantitative risk management. Princeton Series in Finance. Princeton
University Press, Princeton, NJ, revised edition, 2015.
[30] A. J. McNeil and J. Nešlehová. Multivariate Archimedean copulas, d-monotone functions and `1 -norm symmetric distributions. Ann. Statist., 37:3059–3097, 2009.
[31] A. J. McNeil and J. Nešlehová. From Archimedean to Liouville copulas. J. Multivariate Anal., 101:1772–1790,
2010.
[32] N. Nolde. The effect of aggregation on extremes from asymptotically independent light-tailed risks. Extremes,
17(4):615–631, 2014.
[33] T. Opitz. Extremal processes: Elliptical domain of attraction and a spectral representation. J. Multivariate Anal.,
122:409 – 413, 2013.
[34] R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2016.
[35] S. I. Resnick. Extreme values, regular variation, and point processes. Springer-Verlag, Berlin; New York, 1987.
[36] M. Schlather. Models for stationary max-stable random fields. Extremes, 5(1):33–44, 2002.
[37] J. Segers. Max-stable models for multivariate extremes. REVSTAT, 10(1):61–82, 2012.
[38] B. D. Sivazlian. On a multivariate extension of the gamma and beta distributions. SIAM J. Appl. Math.,
41(2):205–209, 1981.
[39] A. Sklar. Fonctions de répartition à n dimensions et leurs marges. Publ. Inst. Statist. Univ. Paris, 8:229–231,
1959.
[40] R. L. Smith, J. A. Tawn, and S. G. Coles. Markov chain models for threshold exceedances. Biometrika,
84(2):249–268, 1997.
[41] E. W. Stacy. A generalization of the gamma distribution. Ann. Math. Stat., 33(3):1187–1192, 1962.
[42] J. L. Wadsworth. Exploiting structure of maximum likelihood estimators for extreme value threshold selection.
Technometrics, 58(1):116–126, 2016.
[43] R. E. Williamson. Multiply monotone functions and their Laplace transforms. Duke Math. J., 23:189–207, 1956.

21

Extremal attractors of Liouville copulas

Appendix A. : Proofs from Section 2
Proof of Proposition 2. To prove parts (a) and (b), recall that 1/Xi is distributed as 1/(RDi ), where Di ∼ Beta(αi , ᾱ−αi )
is independent of R. Furthermore, it is easy to show that 1/Di ∈ M(Φαi ), which implies that E(1/Dβi ) < ∞ for any
β < αi . The extremal behavior of 1/Xi will thus be determined by the extremal behavior of either 1/R or 1/Di ,
depending on which one has a heavier tail. Indeed, Breiman’s Lemma [5] implies that 1/Xi ∈ M(Φρ ) if 1/R ∈ M(Φρ )
for some ρ < αi and that 1/Xi ∈ M(Φαi ) if E(1/Rαi +ε ) < ∞ for some ε > 0. Finally, the fact that 1/Xi ∈ M(Φαi ) when
1/R ∈ M(Φαi ) follows directly from the Corollary to Theorem 3 in [11].
The following lemma is a side result of Proposition 2, which is needed in the subsequent proofs.
Lemma 1 Suppose that X = RDα . If 1/R ∈ M(Φαi ) for some i ∈ {1, . . . , d}, then
lim

x→∞

Pr (1/R > x)
= 0.
Pr (1/Xi > x)

Proof of Lemma 1.
Because 1/R ∈ M(Φαi ), Pr(1/R > x) is regularly varying with index −αi . In particular, for any b ∈ (0, 1),
Pr(1/R > xb)/ Pr(1/R > x) → b−αi as x → ∞. An application of Fatou’s lemma thus gives
Z 1
Z 1 −1
Pr(1/R > xb)
b (1 − b)ᾱ−αi −1
Pr(1/Xi > x)
= lim inf
fDi (b) db ≥
db = ∞
lim inf
x→∞
x→∞ Pr(1/R > x)
B(αi , ᾱ − αi )
0 Pr(1/R > x)
0
and hence the result.

Appendix B. : Proofs from Section 3
First recall the following property of the Dirichlet distribution, which is easily shown using the transformation formula
for Lebesgue densities.
Lemma 2 Let Dα be a Dirichlet random vector with parameters α. Then for any 2 ≤ k ≤ d and any collection of
distinct indices 1 ≤ i1 < · · · < ik ≤ d,
(Di1 , . . . , Dik ) = Bi1 ,...,ik × D(αi1 ,...,αik ) ,
d

where Bi1 ,...,ik ∼ Beta(αi1 + · · · + αik , ᾱ − (αi1 + · · · + αik )) is independent of the k-variate Dirichlet vector D(αi1 ,...,αik )
with parameters (αi1 , . . . , αik ).
d

d

Proof of Theorem 1. In order to prove part (a), recall that kXk = R is independent of X/kXk = Dα . Because
R ∈ M(Φρ ), there exists a sequence (bn ) of constants in (0, ∞) such that, for any Borel set B ⊆ Sd and any r > 0,
!
X
lim n Pr kXk > bn r,
∈ B = lim n Pr(R > bn r) Pr(Dα ∈ B) = r−ρ Pr( Dα ∈ B).
n→∞
n→∞
kXk
By Corollary 5.18 in [35], X ∈ M(H0 ) where for all x ∈ Rd+ ,
 
 ρ

Dρd 
 
 D1




H0 (x) = exp −E 
max  xρ , . . . , xρ 
 .
1
d
Let B(·, ·) denote the Beta function. The univariate margins of H0 are given, for all i = 1, . . . , d and x > 0, by
)
(
)
(
n
o
B(ρ + αi , ᾱ − αi )
Γ(αi + ρ)Γ(ᾱ)
F0i (x) = exp −x−ρ E(Dρi ) = exp −x−ρ
= exp −x−ρ
B(αi , ᾱ − αi )
Γ(ᾱ + ρ)Γ(αi )
for i ∈ {1, . . . , d}. The copula of H0 then satisfies, for all u ∈ [0, 1]d ,
C0 (u) =

−1
−1
H0 {F01
(u1 ), . . . , F0d
(ud )}

"
(
)#!
(− log ui )Γ(αi )Dρi
Γ(ᾱ + ρ)
= exp −
E max
.
1≤i≤d
Γ(ᾱ)
Γ(αi + ρ)

22

L. R. Belzile and J. G. Nešlehová

By Eq. (2), the stable tail dependence function of C0 thus indeed equals, for all x ∈ Rd+ ,



ρ

xd Γ(αd )Dρd 
Γ(ᾱ + ρ) 
 x1 Γ(α1 )D1

E max 
,...,
`(x) =
 .

Γ(ᾱ)
Γ(α1 + ρ)
Γ(αd + ρ) 
The part (b) follows directly from Proposition 2.2 in [19] upon setting p = 1 and taking, for i = 1, . . . , d and
j = 1, . . . , d, λi j = 1 whenever i = j and λi j = 0 otherwise.
To prove part (c), recall first that from Proposition 1, for i = 1, . . . , d, Xi ∈ M(Ψρ+ᾱ−αi ) and hence there exist
sequences (ani ) ∈ (0, ∞), (bni ) ∈ R, such that for all x ∈ R,
lim n Pr(Xi > ani x + bni ) = − log{Ψρ+ᾱ−αi (x)}.

n→∞

Next, observe that as in the proof of Proposition 5.27 in [35], X ∈ M(H0 ) follows if for all 1 ≤ i < j ≤ d and xk such
that Ψρ+ᾱ−αk (xk ) > 0 for k = i, j,
lim n Pr(Xi > ani xi + bni , X j > an j x j + bn j ) = 0.

n→∞

(B.1)

To prove that (B.1) indeed holds, it suffices to assume that d = 2. This is because for arbitrary indices 1 ≤ i <
d
d
j ≤ d, Lemma 2 implies that (Xi , X j ) = R∗ (B, 1 − B), where B ∼ Beta(αi , α j ), R∗ = RY is independent of B and
∗
Y ∼ Beta(αi + α j , ᾱ − αi − α j ) is independent of B and R. Because Pr(R ≤ 0) = 0, Theorem 4.5 in [20] implies that
R∗ ∈ M(Ψρ+ᾱ−αi −α j ) when R ∈ M(Ψρ ) for some ρ > 0. Thus suppose that d = 2 and write (D1 , D2 ) ≡ (B, 1 − B),
where B ∼ Beta(α1 , α2 ). Fix arbitrary x1 , x2 ∈ R are such that Ψρ+ᾱ−αi (xi ) > 0 for i = 1, 2. Then because for any
a, c > 0 and b ∈ (0, 1), max{a/b, c/(1 − b)} ≥ a + c, one has
(
!)
an1 x1 + bn1 an2 x2 + bn2
0 ≤ Pr(X1 > an1 x1 + bn1 , X2 > an2 x2 + bn2 ) = Pr R > max
,
B
1−B
≤ Pr(R > an1 x1 + bn1 + an2 x2 + bn2 ).
In order to prove Eq. (B.1), it thus suffices to show that
lim n Pr(R > an1 x1 + bn1 + an2 x2 + bn2 ) = 0.

n→∞

(B.2)

This however follows immediately from the fact that if R ∈ M(Ψρ ) for some ρ > 0, the upper end-point r of R,
viz. r = sup{x : Pr(R ≤ x) < 1}, is finite. Because for i = 1, 2, r is also the upper endpoint of Xi , ani xi + bni → r
as n → ∞. This means that there exists n0 ∈ N so that for all n ≥ n0 , an1 x1 + bn1 + an2 x2 + bn2 > r and Pr(R >
an1 x1 + bn1 + an2 x2 + bn2 ) = 0. This proves Eq. (B.1) and hence also Theorem 1 (c). Note that alternatively, part (c)
could be proved using Theorem 2.1 in [19] similarly to the proof of Proposition 2.2 therein.
The proof of Theorem 2 requires the following technical lemma.
Lemma 3 Suppose that Dα = (D1 , . . . , Dd ) is a Dirichlet random vector with parameters α. Further let R be a
positive random variable independent of Dα such that Pr(R ≤ 0) = 0, and let X = R Dα . Then for any 1 ≤ i < j ≤ d
and any xi , x j ∈ (0, ∞),
!
1
1
lim n Pr
> ani xi ,
> an j x j = 0
n→∞
Xi
Xj
if either:
(i) 1/R ∈ M(Φρ ) with ρ ∈ [αi ∧ α j , α1 ∨ α2 ], and for k = i, j, (ank ) is a sequence of positive constants such that
n Pr(1/Xk > ank xk ) → xk−(αk ∧ρ) as n → ∞;
(ii) E(1/Rβ ) < ∞ for some β > αi ∨ α j and for k = i, j, (ank ) is a sequence of positive constants such that
n Pr(1/Xk > ank xk ) → xk−αk as n → ∞.

23

Extremal attractors of Liouville copulas
d

Proof of Lemma 3. Observe first that when d > 2, Lemma 2 implies that (Xi , X j ) = R∗ (B, 1 − B), where R∗ ⊥
⊥ B,
B ∼ Beta(αi , α j ) and R∗ = RY, with Y ⊥
⊥ R and Y ∼ Beta(αi + α j , ᾱ − αi − α j ). Now note that 1/Y ∈ M(Φαi +α j ). Thus
if 1/R ∈ M(Φρ ) for ρ ∈ [αi ∧ α j , α1 ∨ α2 ], ρ < αi + α j and Breiman’s Lemma implies that 1/R∗ ∈ M(Φρ ). Further, if
E(1/Rβ ) < ∞ for some β ∈ (αi ∨ α j , αi + α j ), E{1/(R∗ )β } < ∞ given that E(1/Y β ) < ∞. We can thus assume without
loss of generality that d = 2 and α1 ≤ α2 ; we shall also write (D1 , D2 ) ≡ (B, 1 − B), where B ∼ Beta(α1 , α2 ).
To prove part (i), note first that the existence of the sequences (ank ), k = i, j, follows from Proposition 2, by which
1/Xk ∈ M(Φρ∧αk ) for k = i, j, and the Poisson approximation [12, Proposition 3.1.1]. Next, observe that for any
constants a, c > 0 and any b ∈ (0, 1),
ac
≤ ab ∨ c(1 − b) < a ∨ c.
(B.3)
a+c
Indeed, when b < c/(a + c), ab ∨ c(1 − b) = c(1 − b) and c(1 − b) ∈ (ac/(a + c), c), while when b ≥ c/(a + c),
ab ∨ c(1 − b) = ab and ab ∈ [ac/(a + c), a). To show the claim in part (i), distinguish the cases below:
d

Case I. α1 = α2 . Here, ρ = α1 = α2 and X1 = X2 , so that an1 /an2 → 1 by the Convergence to Types Theorem [35].
By Eq. (B.3),
#
!
"
!
1
an1 an2 x1 x2
1
1
1
> max{an1 x1 B, an2 x2 (1 − B)} ≤ n Pr
>
> an1 x1 ,
> an2 x2 = n Pr
. (B.4)
0 ≤ n Pr
X1
X2
R
R an1 x1 + an2 x2
Because (x1 x2 )/{(an1 /an2 )x1 + x2 } → (x1 x2 )/(x1 + x2 ) as n → ∞,
!
!−α1
1
an1 an2 x1 x2
x1 x2
lim n Pr
>
=
.
n→∞
X1 an1 x1 + an2 x2
x1 + x2
Furthermore, by Lemma 1, given that (an1 an2 x1 x2 )/(an1 x1 + an2 x2 ) → ∞ as n → ∞,


Pr 1/R > aan1n1xa1n2+axn21 xx22
 = 0,
lim 
n→∞ Pr 1/X > an1 an2 x1 x2
1
an1 x1 +an2 x2
so that the right-hand side in Eq. (B.4) tends to 0 as n → ∞, and this implies the claim.
Case II. α1 < α2 and ρ = α2 . Then for i = 1, 2, there exists a slowly varying function Li such that ani = n1/αi Li (n).
Hence an2 /an1 → 0 and (x1 x2 )/{x1 + x2 (an2 /an1 )} → x2 as n → ∞. Consequently,
!
1
an1 an2 x1 x2
lim n Pr
>
= x2−α2 .
n→∞
X2 an1 x1 + an2 x2
Moreover, by Lemma 1, given that (an1 an2 x1 x2 )/(an1 x1 + an2 x2 ) → ∞ as n → ∞,


Pr 1/R > aan1n1xa1n2+axn21 xx22
 = 0,
lim 
n→∞ Pr 1/X > an1 an2 x1 x2
2
an1 x1 +an2 x2
so that again the right-hand side in Eq. (B.4) tends to 0 as n → ∞.
Case III. α1 < α2 and ρ ∈ [α1 , α2 ). In this case, 1/X1 ∈ M(Φα1 ) and 1/X2 ∈ M(Φρ ). Therefore, either directly
when ρ > α1 or by Lemma 1, one can easily deduce that
lim

x→∞

Pr(1/R > x)
= 0.
Pr(1/X1 > x)

(B.5)

At the same time, Breiman’s Lemma [5] implies that
lim

x→∞

Pr(1/R > x)
1
B(α1 , α2 )
=
=
.
ρ
Pr(1/X2 > x) E{1/(1 − B) } B(α1 , α2 − ρ)

(B.6)

24

L. R. Belzile and J. G. Nešlehová

Hence, for any b ∈ (0, 1), the limit of n Pr {1/R > an2 x2 (1 − b)} as n → ∞ equals
lim n Pr {1/X2 > an2 x2 (1 − b)}

n→∞

B(α1 , α2 )
Pr {1/R > an2 x2 (1 − b)}
= {x2 (1 − b)}−ρ
Pr {1/X2 > an2 x2 (1 − b)}
B(α1 , α2 − ρ)

so that
lim

n→∞

Z
0

1

bα1 −1 (1 − b)α2 −1
db = n Pr(1/X2 > an2 x2 ) = x2−ρ
n Pr{1/R > an2 x2 (1 − b)}
B(α1 , α2 )
Z 1 α1 −1
Z 1
b
(1 − b)α2 −ρ−1
bα1 −1 (1 − b)α2 −1
= x2−ρ
db =
lim n Pr{1/R > an2 x2 (1 − b)}
db. (B.7)
B(α1 , α2 − ρ)
B(α1 , α2 )
0
0 n→∞

Given that for any b ∈ (0, 1), Pr{1/R > an1 x1 b, 1/R > an2 x2 (1 − b)} ≤ Pr{1/R > an2 x2 (1 − b)},
Z 1

 bα1 −1 (1 − b)α2 −1
lim inf n Pr{1/R > an2 x2 (1 − b)} − Pr{1/R > an1 x1 b, 1/R > an2 x2 (1 − b)}
db
n→∞
B(α1 , α2 )
0
Z 1

 bα1 −1 (1 − b)α2 −1
n Pr{1/R > an2 x2 (1 − b)} − Pr{1/R > an1 x1 b, 1/R > an2 x2 (1 − b)}
≤ lim inf
db
n→∞
B(α1 , α2 )
0
by Fatou’s Lemma. Because of Eq. (B.7), this inequality simplifies to
Z 1

 bα1 −1 (1 − b)α2 −1
db
x2−ρ −
lim sup n Pr{1/R > an1 x1 b, 1/R > an2 x2 (1 − b)}
B(α1 , α2 )
n→∞
0
Z 1
bα1 −1 (1 − b)α2 −1
−ρ
≤ x2 − lim sup
n Pr{1/R > an1 x1 b, 1/R > an2 x2 (1 − b)}
db
B(α1 , α2 )
n→∞
0
and hence

0 ≤ lim sup n Pr(1/X1 > an1 x1 , 1/X2 > an2 x2 )
n→∞

Z
≤
0

1


 bα1 −1 (1 − b)α2 −1
db.
lim sup n Pr{1/R > an1 x1 b, 1/R > an2 x2 (1 − b)}
B(α1 , α2 )
n→∞

To show the desired claim, it thus suffices to show that for arbitrary b ∈ (0, 1),
lim n Pr{1/R > an1 x1 b, 1/R > an2 x2 (1 − b)} = 0.

(B.8)

n→∞

To this end, fix b ∈ (0, 1) and observe that an1 /an2 → ∞. Indeed, if ρ > α1 , this follows directly from the fact
that an1 = n1/α1 L1 (n) and an2 = n1/ρ L2 (n) for some slowly varying functions L1 , L2 . When ρ = α1 , suppose that
lim inf n→∞ an1 /an2 were finite. Then there exists a subsequence ank 1 /ank 2 such that ank 1 /ank 2 → a as k → ∞ for some
a ∈ [0, ∞). Hence, for a fixed ε > 0 and all k ≥ k0 , ank 1 /ank 2 ≤ a + ε. Using the latter observation and Eq. (B.6),
lim nk Pr(1/R > ank 1 ) ≥ lim nk Pr{1/R > ank 2 (a + ε)}

k→∞

k→∞

= lim nk Pr{1/X2 > ank 2 (a + ε)}
k→∞

Pr{1/R > ank 2 (a + ε)}
B(α1 , α2 )
= (a + ε)−ρ
> 0.
Pr{1/X2 > ank 2 (a + ε)}
B(α1 , α2 − ρ)

At the same time, by Eq. (B.5),
lim nk Pr(1/R > ank 1 ) = lim nk Pr(1/X1 > ank 1 )

k→∞

k→∞

Pr(1/R > ank 1 )
=0
Pr(1/X1 > ank 1 )

and hence a contradiction. Therefore, lim inf n→∞ an1 /an2 = ∞ and hence an1 /an2 → ∞ as n → ∞. Because an1 b >
an2 (1 − b) if and only if b > an2 /(an1 + an2 ) and an2 /(an1 + an2 ) → 0 as n → ∞, there exists n0 such that for all n ≥ n0 ,
n Pr{1/R > an1 x1 b, 1/R > an2 x2 (1 − b)} = n Pr(1/R > an1 x1 b) = n Pr(1/X1 > an1 x1 b)

Pr(1/R > an1 x1 b)
.
Pr(1/X1 > an1 x1 b)

25

Extremal attractors of Liouville copulas
The last expression tends to 0 as n → ∞ by Eq. (B.5) and hence Eq. (B.8) indeed holds.

To prove part (ii), first recall that by Proposition 2 (b), 1/Xi ∈ M(Φαi ), i = 1, 2, and hence the scaling sequences
(an1 ) and (an2 ) indeed exist. Recall that for i = 1, 2, ani = n1/αi Li (n) for some slowly varying function Li . As in
the proof of part (i), n Pr(1/X1 > an1 x1 , 1/X2 > an2 x2 ) can be bounded above by the right-hand side in Eq. (B.4).
Markov’s inequality further implies that for β ∈ (α2 , α1 + α2 ) such that E(1/Rβ ) < ∞,


)β
!

 (an1 x1 + an2 x2 )β E 1/Rβ (
x1
an1 an2 x1 x2
x2
1
β
=
>
≤ nE 1/R
+
n Pr
R an1 x1 + an2 x2
(an1 an2 x1 x2 )β
(x1 x2 )β n1/α2 −1/β L2 (n) n1/α1 −1/β L1 (n)
The right-most expression tends to 0 as n → ∞ because for any i = 1, 2 and ρ > 0, nρ Li (n) → ∞.
Proof of Theorem 2. First note that a positive random vector Y is in the maximum domain of attraction of a multivariate
extreme-value distribution H0 with Fréchet margins if and only if there exist sequences of positive constants (ani ) ∈
(0, ∞), i = 1, . . . , d, so that, for all y ∈ Rd+ ,
lim n{1 − Pr(Y1 ≤ an1 y1 , . . . , Yd ≤ and yd )} =
 d

X

X
lim n 

n→∞ 

n→∞

k=1 1≤i1 <···<ik ≤d





= − log H0 (y).
(−1)k+1 Pr(Yi1 > ani1 yi1 , . . . , Yik > anik yik )



(B.9)

This multivariate version of the Poisson approximation holds by the same argument as in the univariate case [12,
Proposition 3.1.1].
To prove part (a), suppose that 1/R ∈ M(Φρ ) for some ρ ∈ (0, α M ]. By Proposition 2, one then has that for
any i ∈ I1 , 1/(RDi ) ∈ M(Φαi ). For any i ∈ I1 , let (ani ) be a sequence of positive constants such that, for all
x > 0, n Pr{1/(RDi ) > ani x} → x−αi as n → ∞; such a sequence exists by the univariate Poisson approximation [12,
Proposition 3.1.1]. The same result also guarantees the existence of a sequence (an ) of positive constants such that,
for all x > 0, n Pr(1/R > an x) → x−ρ as n → ∞. Now set, for any i ∈ I2 ,

 Γ(αi − ρ)Γ(ᾱ − αi )
Γ(ᾱ)
Γ(ᾱ)/Γ(ᾱ − ρ)
×
=
,
bi = E D−ρ
=
i
Γ(ᾱ − ρ)
Γ(αi )Γ(ᾱ − αi ) Γ(αi )/Γ(αi − ρ)

(B.10)

and define, for any i ∈ I2 and n ∈ N, ani = b1/ρ
i an . As detailed in the proof of Proposition 2 (a), Breiman’s Lemma
then implies that, for all i ∈ I2 and x > 0,
(
)
(
) Pr n 1 > a (b1/ρ x)o
n i
1
1
RD
−ρ
n i
o = x−ρ b−1
lim n Pr
> ani x = lim n Pr
> an (b1/ρ
i bi = x ,
i x)
1
n→∞
n→∞
RDi
R
Pr R > an (b1/ρ
x)
i
given that for all i ∈ I2 , Di ∼ Beta(αi , ᾱ − αi ).
Next, fix an arbitrary x ∈ (0, ∞)d , k ∈ {2, . . . , d} and indices 1 ≤ i1 < · · · < ik ≤ d. To calculate the limit of
n Pr(1/(RDi1 ) > ani1 xi1 , . . . , 1/(RDik ) > anik xik ), two cases must be distinguished:
Case I. {i1 , . . . , ik } ∩ I1 , ∅. In this case, suppose, without loss of generality, that i1 ∈ I1 . Then
!
!
1
1
1
1
0 ≤ n Pr
> ani1 xi1 , . . . ,
> anik xik ≤ n Pr
> ani1 xi1 ,
> ani2 xi2 .
RDi1
RDik
RDi1
RDi2
Now either i2 ∈ I1 , in which case ρ ≥ αi1 ∨ αi2 , or i2 ∈ I2 , so that αi1 ≤ ρ < αi2 . Either way, Lemma 3 implies that
!
1
1
lim n Pr
> ani1 xi1 ,
> ani2 xi2 = 0
n→∞
RDi1
RDi2
and consequently n Pr{1/(RDi1 ) > ani1 xi1 , . . . , 1/(RDik ) > anik xik } → 0 as n → ∞.

26

L. R. Belzile and J. G. Nešlehová

Case II. {i1 , . . . , ik } ∩ I1 = ∅. In this case, let Zi1 ,...,ik = max(xi1 (bi1 )1/ρ Di1 , . . . , xik (bik )1/ρ Dik ) and observe that for
any ε > 0 such that ρ + ε < min(α1 , . . . , αd ),




 1 
 1 
−ρ−ε −(ρ+ε)/ρ 



E  ρ+ε  < ∞.
E  ρ+ε  ≤ xi1 bi1
Zi1 ,...,ik
Di1
Therefore, by Breiman’s Lemma,
lim n Pr

n→∞

!
!


1
1
1
> ani1 xi1 , . . . ,
> anik xik = lim n Pr
> an = E Zi−ρ
,...,i
1
k
n→∞
RDi1
RDik
RZi1 ,...,ik

"(
)−ρ #


1/ρ
= E max xi j bi j Di j
= E  min


−ρ 






 xi j Di j 

 .





1≤ j≤k 
bi j



1≤ j≤k

Putting the above calculations together, one then has, for any x ∈ Rd+ ,
1
1
lim n 1 − Pr
≤ an1 x1 , . . . ,
≤ and xd
n→∞
RD1
RDd
(

!)
=

X

xi−αi +

|I2 |
X

X

k=1 {i1 ,...,ik }⊆I2
i1 <···<ik

i∈I1


−ρ 







x
D


i
i
j
j


(−1)k+1 E  min 





1≤ j≤k 
bi j



Furthermore, one can readily establish by induction that for any t ∈ Rd ,
|I2 |
X

X

(−1)k+1 min(ti1 , . . . , tik ) = max(ti ).
i∈I2

k=1 {i1 ,...,ik }⊆I2
i1 <···<ik

Hence, for any x ∈ Rd+ ,
(
lim n 1 − Pr
n→∞

1
1
≤ an1 x1 , . . . ,
≤ and xd
RD1
RDd

!)
=

X
i∈I1

)#
(
"
(xi Di )−ρ
.
xi−αi + E max
i∈I2
bi

By the multivariate Poisson approximation (B.9), 1/X ∈ M(H0 ), where for all x ∈ Rd+ ,

"
)#
−ρ
 X

(x
D
)
i
i
−α
 .
H0 (x) = exp −
xi i − E max
i∈I2
bi
i∈I1

The univariate margins of H0 are given, for all i ∈ I1 , by F0i (x) = x−αi and for all i ∈ I2 , F0i (x) = exp(−x−ρ ). By
Sklar’s Theorem, the unique copula of H0 is given, for all u ∈ [0, 1]d , by (2), where for all x ∈ Rd+ ,



−ρ 

X

 xi Di 


 .
`(x) =
xi + E 
max 

 i∈I2  bi 

i∈I1

The first expression for ` follows immediately from Eq. (B.10). The second expression is readily verified using
Lemma 2, given the fact that if B ∼ Beta(ᾱ2 , ᾱ − ᾱ2 ), E(B−ρ ) = Γ(ᾱ2 − ρ)Γ(ᾱ)/Γ(ᾱ − ρ)Γ(ᾱ2 ).
To prove part (b), recall that by Proposition 2 (b), 1/Xi ∈ M(Φαi ). Hence, there exist sequences of positive
constants (ani ), i = 1, . . . , d, such that for all i = 1, . . . , d and all x > 0, n Pr(1/(RDi ) > ani x) → x−αi as n → ∞. By
Lemma 2 (ii), it also follows that for arbitrary x ∈ (0, ∞)d , k ∈ {2, . . . , d} and indices 1 ≤ i1 < · · · < ik ≤ d,
!
!
1
1
1
1
0 ≤ lim n Pr
> ani1 xi1 , . . . ,
> anik xik ≤ lim n Pr
> ani1 xi1 ,
> ani2 xi2 = 0.
n→∞
n→∞
RDi1
RDik
RDi1
RDi2
Thus, by Eq. (B.9), 1/X is in the domain of attraction of the multivariate extreme-value distribution given, for all
x ∈ Rd+ , by H0 (x) = exp(−x1−α1 − · · · − xd−αd ), as was to be showed.

27

Extremal attractors of Liouville copulas

Appendix C. : Proofs from Section 4
Proof of Proposition 3. In view of Corollary 1 and Theorem 2 in [27], it only remains to derive the explicit expression
for `nD . Because 1−ψ(1/·) ∈ R−ρ , there exists a slowly varying function L such that for all x > 0, 1−ψ(1/x) = x−ρ L(x).
Given that the distribution function ψ(1/·) is in the domain of attraction of Φρ , the Poisson approximation implies that
there exists a sequence (an ) of positive constants such that, for all x > 0,


lim n 1 − ψ{1/(an x)} = lim n(an x)−ρ L(an x) = x−ρ .
(C.1)
n→∞

n→∞

Furthermore, by Equation (A6) in the proof of Theorem 2 (a) in [27], one has, for any j = 1, . . . , ᾱ − 2,
(−1) j x− j ψ( j) (1/x)
= 1,
x→∞
κ j x−ρ L(x)

(C.2)

lim

where κ j = ρΓ( j − ρ)/Γ(1 − ρ). Now for all i = 1, . . . , d, Eq. (7) yields, for any x > 0,


!)
!
(


α
i −1


j
− j ( j)
X


1
(−1)
(a
x)
ψ
(1/a
x)
1


n
n
1
−
.
= n(an x)−ρ L(an x) 
n Pr
> an x = n 1 − H̄i



−ρ L(a x)


Xi
an x
j!(a
x)


n
n
j=1
Given that an → ∞ as n → ∞, the last expression converges by Equations (C.1) and (C.2) as n → ∞ to








α
α
i −1
i −1




X
X




κ
Γ( j − ρ)
c(αi , −ρ)
j


−ρ 
1
−
x−ρ 
=
x
1
−
ρ
= x−ρ
.











j!
Γ(
j
+
1)Γ(1
−
ρ)
Γ(1 − ρ)




j=1
j=1
The Poisson approximation thus implies that, as n → ∞, for all i = 1, . . . , d and x > 0,
!
(
)
1
n
−ρ c(αi , −ρ)
H̄i
→ exp −x
.
an x
Γ(1 − ρ)

(C.3)

For any x ∈ (0, ∞)d , let 1/(an x) = {1/(an x1 ), . . . , 1/(an xd )} and denote by x̄H the harmonic mean of x, viz. x̄H =
d/(1/x1 + · · · + 1/xd ). From Eq. (6) one then has





− j1 −···− jd






(
!)
!
d
j1 +···+ jd an x̄H
( j1 +···+ jd )


d
j




i


X
Y
−ρ
(−1)
ψ


1
an x̄H
an x̄H 
x̄

d
an x̄H
H
1
−
n 1 − H̄
=n
L






−ρ

a
x̄
a
x̄


n H
n H
an x
d
d 
dxi 

j
!
·
·
·
j
!
L


(
j
,...,
j
)∈I
i=1
1
d
1
d
α


d
d


( j ,..., j ),0
1

d

d

By Eq. (C.2), the right most expression in the curly brackets converges, as n → ∞, to
1−ρ

X
( j1 ,..., jd )∈Iα
( j1 ,..., jd ),0d

d

Γ( j1 + · · · + jd − ρ) Y x̄H
Γ(1 − ρ) j1 ! · · · jd ! i=1 dxi

! ji

=1−ρ

X
( j1 ,..., jd )∈Iα
( j1 ,..., jd ),0d

d

Γ( j1 + · · · + jd − ρ) Y
1
1/xi
Γ(1 − ρ)
Γ(
j
+
1)
1/x
+
· · · + 1/xd
i
1
i=1

Furthermore, Eq. (C.1) implies that, as n → ∞,
!ρ
 a x̄ −ρ  a x̄ 
1
1
n H
n H
n
L
→
+ ··· +
.
d
d
x1
xd
Consequently, as n → ∞, n{1 − H̄(1/an x)} → − log H0 (x), where









j
!ρ 
i

d

X Γ( j1 + · · · + jd − ρ) Y
 1/xi  


1
1 
1






− log H0 (x) =
+ ··· +
1
−
ρ
.


P






d
1


x1
xd 
Γ(1 − ρ)
Γ( ji + 1)



j=1 x j
i=1
(
j
,...,
j
)∈I
1
d
α




( j ,..., j ),0
1

d

d

! ji
.

28

L. R. Belzile and J. G. Nešlehová

By Eq. (B.9), 1/X ∈ M(H0 ). From Eq. (C.3), the univariate margins of H0 are scaled Fréchet, and Sklar’s theorem
implies that the unique copula of H0 is of the form (2) with stable tail dependence function as in Proposition 3.
Proof of Proposition 4. In view of Corollary 2 and Theorem 1 in [27], it only remains to compute the expression for
`pD given in part (a). Suppose that ψ ∈ R−ρ for some ρ > 0. This means that there exists a slowly varying function
such that for all x > 0, ψ(x) = x−ρ L(x). Because ψ is itself a survival function, ψ ∈ M(Φρ ) and by the univariate
Poisson approximation, there exists a sequence (an ) of strictly positive constants such that, for all x > 0,
lim nψ(an x) = x−ρ .

(C.4)

n→∞

Furthermore, by Equation (A1) in the proof of Theorem 1 (a) in [27], one has, for any j = 1, . . . , ᾱ − 1,
(−1) j x j ψ( j) (x)
= c( j, ρ).
x→∞
ψ(x)

(C.5)

lim

Now let X be the Dirichlet random vector with parameters α and radial part R whose Williamson ᾱ-transform is ψ.
Denote the distribution function of X by H and its univariate margins by Fi , i = 1, . . . , d. Then for all i = 1, . . . , d,
Equations (C.4) and (C.5) imply that
lim nF̄i (an x) = lim n

n→∞

n→∞

α
i −1
X
j=0

α −1

i
X
Γ( j + ρ)
x−ρ c(αi , ρ)
(−1) j (an x) j ψ( j) (an x)
= x−ρ
=
j!
Γ(ρ)Γ( j + 1)
Γ(ρ + 1)
j=0

(C.6)

and hence, by the Poisson approximation, Fin (x) → exp{−x−ρ c(αi , ρ)/Γ(ρ + 1)} as n → ∞.
Next, for arbitrary k = 1, . . . , d and 1 ≤ i1 < · · · < ik ≤ d, let I(αi1 ,...,αik ) = {0, . . . , αi1 − 1} × · · · × {0, . . . , αik − 1}.
For any x ∈ (0, ∞)d , Equations (6), (C.4) and (C.5) imply that

n→∞

n→∞

(−1) j1 +···+ jk

( j1 ,..., jk )∈I(αi

1

,...,αi )
k

= (xi1 + · · · + xik )−ρ

X

(−1)

k+1

k=1 1≤i1 <···<ik ≤d

k

X
( j1 ,..., jk )∈I(αi

Therefore, for any x ∈ (0, ∞)d ,
 d


X
lim n 

n→∞ 

ψ( j1 +···+ jk ) {an (xi1 + · · · + xik )} Y
(an xim ) jm
j1 ! · · · jk !
m=1
k

X

lim n Pr(Xi1 > xi1 , . . . , Xik > xik ) = lim n

,...,αi )
1
k

x im
Γ( j1 + · · · + jk + ρ) Y
Γ(ρ) j1 ! · · · jk ! m=1 xi1 + · · · + xik

! jm
.





Pr(Xi1 > an xi1 , . . . , Xik > an xik )
= − log H0 (x),



where
− log H0 (x) =

d
X

X

(−1)k+1 (xi1 + · · · + xik )−ρ

k=1 1≤i1 <···<ik ≤d

k

X
( j1 ,..., jk )∈I(αi

,...,αi )
1
k

xim
Γ( j1 + · · · + jk + ρ) Y
Γ(ρ) j1 ! · · · jk ! m=1 xi1 + · · · + xik

! jm
.

By Eq. (B.9), X ∈ M(H0 ). As argued above, the univariate margins of H0 are given, for all i = 1, . . . , d and x > 0, by
exp{−x−ρ c(αi , ρ)/Γ(ρ + 1)}. Sklar’s theorem thus implies that the unique copula of H0 is of the form (2) with stable
tail dependence function indeed as given by the expression in part (a).

Appendix D. : Proofs from Section 5
Proof of Proposition 5. First, we show that for any ρ > − min(α1 , . . . , αd ), ρ , 0,
 d
−ρ−ᾱ d
)#
"
(
Z
Y
X

xi Dρi
Γ(ᾱ)
1/ρ

E max
=
max(xi ti )  {c(αi , ρ)ti } 
{c(αi , ρ)}αi /ρ (ti )αi /ρ−1 dt.
Qd
d−1
1≤i≤d c(αi , ρ)
|ρ|
Γ(α
)
S
i
d
i=1
i=1
i=1

(D.1)

29

Extremal attractors of Liouville copulas
d

Indeed, using the fact that (D1 , . . . , Dd ) = Z/kZk, where Zi ∼ Ga(αi , 1), i = 1, . . . , d are independent,
)# Z
)
(
"
(
d
Y
e−zi zαi i −1
xi zρi
xi Dρi
=
(z1 + · · · + zd )−ρ
max
dz.
E max
1≤i≤d c(αi , ρ)
Γ(αi )
Rd+ 1≤i≤d c(αi , ρ)
i=1
P
P
Make a change of variable ti = {zρi /c(αi , ρ)}/ dj=1 zρj /c(α j , ρ) for i = 1, . . . , d − 1 and w = dj=1 zρj /c(α j , ρ). For ease of
P
1/ρ
notation, set also td = 1 − d−1
and the absolute value of the Jacobian is
i=1 ti . Then, for i = 1, . . . , d, zi = {c(αi , ρ)ti w}
d

1 d/ρ−1 Y
c(αi , ρ)1/ρ ti1/ρ−1 .
w
|ρ|d
i=1

|J| =
Therefore,
xi Dρi
E max
1≤i≤d c(αi , ρ)
"

(

)#

−ρ d
 d
 Y
X
 {c(αi , ρ)ti }1/ρ 

c(αi , ρ)αi /ρ tiαi /ρ−1
=
max
(x
t
)
Q
i
i


d
d
1≤i≤d
|ρ|
i=1 Γ(αi ) Sd
i=1
i=1
Z ∞
1/ρ Pd
1/ρ
×
wᾱ/ρ−1 e−w i=1 {c(αi ,ρ)ti } dw dt.
1

Z

0

Eq. (D.1) now follows from the fact that
Z

∞

ᾱ/ρ−1 −w1/ρ

w

e

1/ρ
i=1 {c(αi ,ρ)ti }

Pd

0

 d
−ᾱ
X

1/ρ

dw = |ρ| Γ(ᾱ)  {c(αi , ρ)ti }  .
i=1

The expression for hD now follows directly from Eqs. (3) and (D.1), while the formulas for hpD and hnD obtain upon
setting ρ = ρ and ρ = −ρ, respectively.

Appendix E. : Proofs from Section 7
Proof of Proposition 6. For k = 1, . . . , d, the formula for the kth order mixed partial derivatives of `D (1/x) can
be established from Eq. (12). Indeed, if V denotes a random vector with independent scaled Gamma components
Vi ∼ sGa{1/c(αi , ρ), 1/ρ, αi }, then the point process representation Eq. (12) implies that, for all x ∈ Rd+ ,
(
)
Z ∞ 
Z ∞ 
d

Y


Vi
1
1
D
1 −
` (1/x) =
Pr
F xi t;
(18)
> xi for at least one i ∈ {1, . . . , d} dt =
, , αi  dt.
t
c(αi , ρ) ρ
0
0
i=1
For any k = 1, . . . , d, the expression on the right-hand side of Eq. (18) can be differentiated with respect to x1 , . . . , xk
under the integral sign. This gives the formulas for ∂`D (1/x)/∂x1 . . . ∂xk . When k = d, Eq. (18) implies that
∂d `D (1/x)
=−
∂x1 · · · ∂xd

(
d
Y
t
f xi t;

)
1
1
, , αi dt
c(αi , ρ) ρ
0
i=1


Z
d
d
X


1 Y c(α j , ρ){c(α j , ρ)x j }α j /ρ−1 ∞ α j /ρ
1/ρ
1/ρ
= d
t
exp −t
{c(α j , ρ)x j }  dt
Γ(α j )
ρ j=1
0
j=1

=

Z

∞

d

Γ(ᾱ + ρ)
ρd−1

d
Y
c(α j , ρ){c(α j , ρ)x j }α j /ρ−1

iᾱ+ρ
d
1/ρ
j=1
j=1 {c(α j , ρ)x j }

hP

Γ(α j )

,

P
where the last equality follows upon making the change of variable u = dj=1 {c(α j , ρ)x j }1/ρ t1/ρ . Alternatively, Theorem 1 in [7] implies that that the dth order mixed partial derivative of ` D (1/x) equals −dkxk−d−1 hD (x/kxk; ρ, α), which
indeed simplifies to −dhD (x; ρ, α) given that hD (x/kxk; ρ, α) = kxkd+1 hD (x; ρ, α).
Finally, the formulas for F(x; a, b, c) follow immediately from the fact that the scaled Gamma distribution is also
the distribution of the random variable aZ 1/b , where Z is Gamma with shape c and unit scaling.

30

L. R. Belzile and J. G. Nešlehová
Derivation of the gradient score. Straightforward calculations show that
!
(ᾱ + ρ)c(αi , ρ)1/ρ xi1/ρ−1
∂ log dhD (x)
1
αi
= − Pd
−1
+
∂xi
ρ
xi
ρ j=1 {c(α j , ρ)x j }1/ρ


!
!
c(αk , ρ)1/ρ xk1/ρ−1 
(ᾱ + ρ)c(αi , ρ)1/ρ xi1/ρ−1  1
αi
Iik
Iik
∂2 log dhD (x)


−1 2.
= − Pd
− Pd
 − 1
 −
1/ρ
1/ρ
∂xi ∂xk
ρ
xi ρ j=1 {c(α j , ρ)x j }
ρ
xi
ρ j=1 {c(α j , ρ)x j }

