Simulation algorithms for hierarchical Archimedean copulas beyond the completely monotone case

Abstract Two simulation algorithms for hierarchical Archimedean copulas in the case when intra-group generators are not necessarily completely monotone are presented. Both generalize existing algorithms for the completely monotone case. The underlying stochastic models for both algorithms arise as a particular instance of a more general probability space studied recently in Ressel, P. (2018): A multivariate version of Williamson’s theorem, ℓ1-symmetric survival functions, and generalized Archimedean copulas. Depend. Model. 6, 356–368. On this probability space the inter-group dependence need not be Archimedean, however, we highlight two particular circumstances that guarantee that a hierarchical Archimedean copula is obtained.


Introduction
Copulas are a popular tool to model dependence outside the Gaussian world. Archimedean copulas form one of the best-studied (semi-) parametric family of copulas, due to their convenient analytical expression and the availability of e cient and exact simulation algorithms. The concept of hierarchical Archimedean copulas as a generalization of Archimedean copulas was initiated in [14, p. 88]. The idea is to model exchangeable subgroups via well-understood Archimedean copulas and introduce dependence between them again by Archimedean copulas. A full understanding of the (non-exchangeable) dependence structure can then be decomposed into accessible (exchangeable) Archimedean building blocks. As a result, hierarchical Archimedean copulas have been found useful in mathematical nance, e.g. [7,13,23,27], and insurance risk, e.g. [1]. Furthermore, there is extensive literature on the design and estimation of hierarchical Archimedean copulas, see [2,4,5,19,24,25].
The construction of hierarchical Archimedean copulas requires the parameterizing functions of the involved Archimedean copulas to satisfy non-trivial compatibility conditions. However, [17] presented such su cient and useful compatibility conditions in arbitrary dimension, which were investigated further in [6,8,10,11], who in particular constructed many examples of compatible parameterizing functions. More recently, [1] present a second, di erent compatibility condition. Both of these compatibility conditions are based on parameterizing functions that are all completely monotone. It is well-known that this assumption in the su cient condition of [17] is unnecessarily strict. For instance, [9, Remark 7.3.1, p. 153] and [12] indicate that the completely monotone conditions on the parameterizing functions can be weakened. Indeed, [22,Corollary 2.3] derives su cient compatibility conditions which are signi cantly weaker than those of [17]. Un-fortunately, the weaker conditions of [22], which dispense with the unnecessarily strict completely monotone condition, are derived on purely analytical grounds. Simulation algorithms based on these weaker conditions are unknown to date. In contrast, the compatibility conditions in [1,17] are derived from a probabilistic construction (based on complete monotonicity), which can be leveraged to construct simulation algorithms for the respective hierarchical Archimedean copulas.
This article illustrates that two di erent special cases of a probabilistic model in [21] give rise to hierarchical Archimedean copulas. In the rst case we obtain copulas based on parameterizing functions satisfying conditions weaker than those of [17], and stronger than those of [22]. In the second case we achieve a signi cant weakening of the conditions in [1]. This insight implies the following contributions to the existing literature: • Simulation algorithms beyond the well-known completely monotone case are derived based on the presented probabilistic model. This paves the way to making use of these copulas in Monte Carlo studies. Until now, it was not possible to simulate random vectors from these copulas. • By weakening the conditions of [1], we construct hierarchical Archimedean copulas which are outside the realm of the conditions in [22], see Example 3.4 below. • Within the exchangeable groups, the presented probabilistic model agrees with the well-known -norm symmetric model of [18] for exchangeable Archimedean copulas, in which dependence is caused by random radial variables. Archimedean dependence between groups relies on a tricky, non-arbitrary, dependence model of these radial variables. This sheds light on how hierarchical Archimedean copulas arise as dependence structures of a natural and easy-to-grasp probabilistic model.
The remainder of the article is organized as follows. Section 2 repeats known simulation algorithms and required mathematical background. Section 3 presents simulation algorithms based on weaker assumptions than previous literature. Section 4 nally concludes.

Mathematical preliminaries and existing algorithms
Throughout, d = denotes equality in distribution, the symbol ∼ means "is distributed according to", and f (k) denotes the k-th derivative of a real-valued function f in one variable, provided it exists. We further denote by U[ , ] the uniform distribution on the unit interval [ , ] and by Bernoulli(p) the Bernoulli distribution with success probability p ∈ [ , ]. Finally, all random variables to be introduced in the sequel are formally de ned on some probability space (Ω, F, P), and we denote the expected value with respect to the probability measure P by E[.].

. Mathematical background
We recall that a function φ : being non-increasing and convex on (a, b). For further reference we introduce the notation for φ ∈ Φn, where φ − denotes the generalized inverse of φ. We recall from [18] that for each φ ∈ Φn there is a random variable R taking values in ( , ∞), uniquely determined in law, such that and φ is called the Williamson-n-transform of R, named after [26]. The Williamson inversion formula states the distribution function of R in terms of φ as where φ (n− ) denotes the right-continuous version of the almost everywhere existing derivative of φ (n− ) . If S = (S , . . . , Sn) is independent of R and uniformly distributed on the n-dimensional unit simplex, then the random vector has distribution function Cφ,n, which implies a simulation algorithm for Archimedean copulas.  . . , u J ) accordingly. We will later denote other partitioned random vectors using the same notation, which will be made clear from the context. As a convenience, we will sometimes denote the group of all indices corresponding to the sub-vector is called a hierarchical Archimedean copula, provided it is a proper distribution function, which depends on properties of the functions φ ,J , φ ,d . . . , φ J,d J . The function φ ,J is called an outer generator, since it speci es the inter-group dependence, and the functions φ j,d j are called inner generators, since they specify the intra-group dependence. Known simulation algorithms for C rely on the assumption that all involved generators are in Φ∞, hence completely monotone. The present article relaxes this assumption by presenting simulation algorithms for the case when the inner generators are allowed to be proper d j -monotone.
For later reference, we recall that a function g : [21]. Furthermore, if g is n-alternating for all n ≥ and g( ) = , then g is called a Bernstein function. Recall that a function g is a Bernstein function if and only if exp(−m g) is a Laplace transform for each m > , i.e. the probability law on [ , ∞] speci ed by the Laplace transform exp(−g) is in nitely divisible.

. Overview of known and new results
The following list contains analytical conditions on inner and outer generators which are relevant with respect to compatibility for nesting Archimedean copulas: • Conditions on the outer generators: ,J ∈ Φ∞, and a random variable M with Laplace transform φ ,J is assumed to take values in the set of natural numbers N.
• Conditions on the inner generators (given some outer generator φ ,J ): ) an unbounded Bernstein function that is continuous¹, for each j = , . . . , J. Figure 1 illustrates the di erent levels of generality of these conditions. In particular, we have O ⊂ O ⊂ O and I ⊂ I ⊂ I , as well as I ⊂ I . The following list provides an overview of known results ( ) and contributions of the present article (*) with regards to the nesting and simulation of hierarchical Archimedean copulas: [17] shows that O ∩ I is su cient, and provides a simulation algorithm under these conditions. [22] shows that O ∩ I is su cient, but no simulation algorithm under these conditions is known to date. * Lemma 3.1 below provides a simulation algorithm under the su cient condition O ∩ I , which lies between O ∩ I and O ∩ I . [1] show that O ∩ I is su cient, and provide a simulation algorithm under these conditions. * Lemma 3.3 below shows that O ∩ I is su cient, and provides a simulation algorithm under these conditions. re-interpreted and explored further in [6], leading to the simulation of a random vector U ∼ C as follows: (i) Simulate M with Laplace transform φ ,J , and independently, a list of d iid unit exponentials ϵ = ( ϵ , . . . , ϵ J ) = (ϵ , . . . , ϵ d ).
(ii) For j = , . . . , J, simulate a positive random variable M j with Laplace transform exp(−M g j ). (Notice that M j is positive and nite, since g j is assumed to be unbounded and continuous.) On the other hand, [1] consider the conditions O ∩ I . Under these, a random vector U ∼ C can be simulated as follows: For both conditions O ∩ I and O ∩ I , all involved generators are completely monotone. Consequently, the respective probabilistic models must be seen as hierarchical generalizations of the conditionally iid model (3) rather than the more general model (2). In the simulation algorithms derived below, we stick with the completely monotone assumptions O and O on the outer generators, but weaken the conditions I and I on the inner generators. Both of our resulting simulation algorithms arise as two distinct special cases of a more general stochastic model considered in [21]. In general this model only leads to something "similar" to a hierarchical Archimedean copula. More precisely, Let R = (R , . . . , R J ) be an arbitrary random vector on ( , ∞) J and de ne Then the survival copula of with S j , j = , . . . , J, independent and each S j uniformly distributed on the d j -dimensional unit simplex (independent of R) is given by In general, copulas of the form C f do not coincide with hierarchical Archimedean copulas, but there is an intersection. Whereas the survival copulas of the sub-vectors X j = R j S j are always Archimedean, the dependence between components from distinct groups needs not be, but can be, Archimedean, relying on a clever modeling of the dependence between the components of R. On a high level, the simulation algorithm relying on O ∩ I relies on the choice with M ∼ φ ,J and, independently, {B (k) } k∈N , . . . , {B (k) J } k∈N independent iid sequences of random variates with Laplace transforms given by h , . . . , h J , respectively, and ϵ = ( ϵ , . . . , ϵ J ) a partitioned vector of d = d + . . . + d J iid unit exponentials independent of all other variables. As can be seen from (6) and (7), both models coincide when the probability laws of B ( ) , . . . , B ( ) J are in nitely divisible and M takes values in N, see Remark 2.1. In the upcoming section, we embed both cases in a more general model, which intuitively is a hierarchical generalization of (2) rather than the conditionally iid model (3).

Simulation algorithms beyond the completely monotone case
In this section, we present two di erent simulation algorithms. In section 3.1 we generalize the algorithm based on conditions O ∩ I by relaxing the condition I to I , and include the possibility for inner generators φ j,d j that are not completely monotone. In section 3.
has as distribution function the hierarchical Archimedean copula C in (4).
Proof. We rst observe that the function lies in Φ d j for each j = , . . . , J and each m > . To see this, observe that the functions exp(m .) and −g j (−.) are d j -absolutely monotone, so is their composition by [20,Theorem 12]. This implies that c j,m (−.) is d j -absolutely monotone, which means that c j,m is d j -monotone. Since further g j is assumed to be a bijection, in particular is unbounded, we realize c j,m ∈ Φ d j . Consequently, there exists a unique random variable R j,m whose Williamson-d j -transform equals c j,m . We consider X = ( X , . . . , X J ), where We compute the survival function of X as follows, where (*) follows from the general result on (exchangeable) Archimedean copulas in [18]: where φ j,d j = φ ,J • g j is used in the last equality. This implies for the one-dimensional margin i ∈ G j that Observing that U i = φ ,J • g j (X i ) for each i = , . . . , d proves the claim.

Example 3.2.
To illustrate we consider the case J = d = d = . We specify g j (x) := min{x, a j + b j x} for constants a j > and > b j > , then g j is clearly not a Bernstein function, but is concave, and hence -alternating. Therefore, it satis es I but not I . Further, we observe from (1) that R j,m with Williamson-transform exp(−m g j ) satis es Apparently, and so R j,m has an atom at a j /( − b j ). Since the position of this atom a j /( − b j ) does not depend on m, the resulting inner Archimedean copula with generator φ j,d j will exhibit a singular component, see Figure 2. If F n,λ denotes the Erlang distribution function with parameters λ > and n ∈ N, that is the law of the sum of n iid exponential variates with rate λ, then we observe that Thus, an exact simulation algorithm for R j,m is given as follows: and distinguish three cases: Now assume that φ , (x) = λ/(λ + x), so that the outer generator is of Clayton type with parameter . Notice that the choice of λ has no in uence on the outer Archimedean copula Cφ , , but it does have an e ect on the inner Archimedean copulas Cφ , and Cφ , . Figure 2 visualizes bivariate scatter plots of some two-dimensional margins from the constructed hierarchical Archimedean copula. Inter-group: (Clayton family) Intra-group: has as distribution function the hierarchical Archimedean copula C in (4).
Proof. We observe that the function lies in Φ d j for each j = , . . . , J and each m ∈ N. To see this, observe that the functions x → x m for m ∈ N and h j (−.) are d j -absolutely monotone, so is their composition by [20,Theorem 12]. This implies that c j,m (−.) is d j -absolutely monotone, which means that c j,m is d j -monotone, and therefore c j,m ∈ Φ d j . Consequently, there exists a unique random variable R j,m whose Williamson-d j -transform equals c j,m . Considering X = ( X , . . . , X J ), where we compute the survival function of X as follows, where (*) follows from the general result on (exchangeable) Archimedean copulas in [18]: • h j is used in the last equality. This implies for the one-dimensional margin i ∈ G j that Observing  For m = , we obtain which are both distribution functions that can be simulated exactly and quickly via the inversion method, in particular . Now assume that φ , (x) = − − exp(−x) /θ (Joe's family), which means that M ∼ φ , can be simulated according to the algorithm presented in [10], see also [15,Algorithm 2.2,p. 71]. Figure 4 visualizes bivariate scatter plots of some two-dimensional margins from the constructed hierarchical Archimedean copula. , a 2 =3.00 Figure 4: Bivariate scatter plots ( samples) to illustrate the hierarchical Archimedean copula of Example 3.4. Inter-group dependence is visualized in the scatter plot for (U , U ), speci ed by φ , (x) = − − exp(−x) /θ , corresponding to Joe's family. Intra-group dependence is based on generators that are not completely monotone, and is visualized by the scatter plots for U = (U , U ) and U = (U , U ).

Conclusion
Simulation algorithms have been presented for two-level hierarchical Archimedean copulas, for which the intra-group dependence can be parameterized by generators that are not completely monotone, thus generalizing existing simulation algorithms. The algorithms have been illustrated along with two explicit examples for hierarchical Archimedean copulas. The outer generators, however, were still assumed to be completely monotone. The simulation from hierarchical Archimedean copulas with outer generators that are not completely monotone remains an interesting open problem, relevant also for the introduction of further levels of nesting.