Dependence properties and Bayesian inference for asymmetric multivariate copulas

TL;DR

This paper explores asymmetric copulas, providing new theoretical insights, dependence properties, and an iterative Bayesian inference method, with applications demonstrated through simulations.

Contribution

It introduces a new iterative representation for Liebscher copulas that ensures uniform margins and facilitates Bayesian inference.

Findings

Derived exact tail dependence expressions.

Characterized copulas with arbitrary singular components.

Developed an ABC sampling scheme for inference.

Abstract

We study a broad class of asymmetric copulas introduced by Liebscher (2008) as a combination of multiple - usually symmetric - copulas. The main thrust of the paper is to provide new theoretical properties including exact tail dependence expressions and stability properties. A subclass of Liebscher copulas obtained by combining Fr\'echet copulas is studied in more details. We establish further dependence properties for copulas of this class and show that they are characterized by an arbitrary number of singular components. Furthermore, we introduce a novel iterative representation for general Liebscher copulas which de facto insures uniform margins, thus relaxing a constraint of Liebscher's original construction. Besides, we show that this iterative construction proves useful for inference by developing an Approximate Bayesian computation sampling scheme. This inferential procedure is demonstrated on simulated data.

Tables

Table 1. Table 1: First two rows: average relative errors ( 22 ) for Kendall’s distribution function and Spearman’s ρ 𝜌 \rho between the observed sample and the samples retained by the ABC procedure, for varying K 𝐾 K (columns). Third row: fraction of times that the same decision is taken (‘reject’ vs ‘do not reject’, at the 5% level) based on 𝑿 m subscript 𝑿 𝑚 \bm{X}_{m} and 𝑿 obs subscript 𝑿 obs \bm{X}_{\text{obs}} . The results are averaged over 20 independent repetitions. Standard deviations in parentheses. All values are in % percent \% .

$K$	2	3	4	5
${\overline{\eta }}_{𝒦}$	1.99 (0.16)	2.32 (0.49)	2.38 (0.40)	2.37 (0.56)
${\overline{\eta }}_{\rho }$	5.44 (3.98)	8.16 (5.62)	6.34 (3.74)	9.66 (4.12)
${\overline{\text{f}}}_{\text{test}}$	16.8 (18.1)	19.4 (22.0)	9.2 (13.2)	13.3 (18.1)

Table 2. Table 2: First two rows: average relative errors ( 22 ) for Kendall’s distribution function and Spearman’s ρ 𝜌 \rho between the observed sample and the samples retained by the ABC procedure, for growing noise (columns). Third row: fraction of times that the same decision is taken at the 5% level based on 𝑿 m subscript 𝑿 𝑚 \bm{X}_{m} and 𝑿 obs subscript 𝑿 obs \bm{X}_{\text{obs}} . All values are in % percent \% .

${\sigma }_a^2$	0	${10}^{-5}$	${10}^{-4}$	${10}^{-3}$	${10}^{-2}$	${10}^{-1}$
${\eta }_{𝒦}$	1.7	1.68	1.77	1.68	2.03	2.81
${\eta }_{\rho }$	4.07	3.77	4.66	4.30	4.88	13.00
${\text{f}}_{\text{test}}$	9.00	1.33	15.7	3.33	2.33	7.00

Equations

$$F(\mathbf{x})=C(F_1(x_1),\dots ,F_d(x_d)),\mathbf{x}=(x_1,\dots ,x_d)\in {\mathbb{R}}^d.$$

$$\mathbf{u}\in [0,1{]}^d\mapsto \tilde{C}(\mathbf{u})=\prod _{k=1}^K{C}_k(g_1^{(k)}(u_1),\dots ,g_d^{(k)}(u_d))$$

$$\prod _{k=1}^K{g}_j^{(k)}(u)=u\text{ for all }u\in [0,1],j\in \{1,\dots ,d\}.$$

$${\Lambda }_L(C;x,y)=\underset{\epsilon \to 0}{\lim }\frac{C(\epsilon x,\epsilon y)}{\epsilon },\text{and}{\Lambda }_U(C;x,y)=x+y+\underset{\epsilon \to 0}{\lim }\frac{C(1-\epsilon x,1-\epsilon y)-1}{\epsilon },$$

$${\lambda }_L(C)=\underset{u\to 0}{\lim }\frac{C(u,u)}{u},{\lambda }_U(C)=2-\underset{u\to 1}{\lim }\frac{C(u,u)-1}{u-1}.$$

$${\Lambda }_L(\tilde{C};x,y)=\prod _{k=1}^K{\Lambda }_L(C_k;x^{{\gamma }^{(k)}},y^{{\gamma }^{(k)}})$$

$${\Lambda }_L(\tilde{C};x,y)=0.$$

$${\Lambda }_U(\tilde{C};x,y)=\sum _{k=1}^K{\Lambda }_U(C_k;{\text{d}}_1^{(k)}x,{\text{d}}_2^{(k)}y)$$

$${\Lambda }_U(\tilde{C};x,y)=\sum _{k=1}^K{\text{d}}^{(k)}{\Lambda }_U(C_k;x,y)$$

$$\mathrm{Pr}(X\le x_1,Y\le x_2)\mathrm{Pr}(X\le y_1,Y\le y_2)\ge \mathrm{Pr}(X\le x_1,Y\le y_2)\mathrm{Pr}(X\le y_1,Y\le x_2).$$

$$\mathrm{Pr}(X\le x,Y\le y)\ge \mathrm{Pr}(X\le x)\mathrm{Pr}(Y\le y)\text{mbox}forall(x,y).$$

$$\mathrm{Pr}(X\le x|Y\le y)\text{ is a decreasing function of }y\text{ for all }x,$$

$$\mathrm{Pr}(Y\le y|X\le x)\text{ is a decreasing function of }x\text{ for all }y.$$

$$C(u,v)/u\text{mbox}isdecreasinginu\text{mbox}forallv\in [0,1],$$

$$C(u,v)/v\text{mbox}isdecreasinginv\text{mbox}forallu\in [0,1],$$

$$\mathrm{Pr}(X>x|Y=y)\text{ is an increasing function of }y\text{ for all }x,$$

$$\mathrm{Pr}(Y>y|X=x)\text{ is an increasing function of }x\text{ for all }y.$$

$$C(u,v)\text{mbox}isaconcavefunctionofu\text{mbox}forallv\in [0,1],$$

$$C(u,v)\text{mbox}isaconcavefunctionofv\text{mbox}forallu\in [0,1],$$

$$g_j^{(k)}(t)=t^{p_j^{(k)}},p_j^{(k)}\in (0,1),\sum _{\ell =1}^K{p}_j^{(\ell )}=1.$$

$$C_{\#}^n(u_1^{1/n},\dots ,u_d^{1/n})=C_{\#}(u_1,\dots ,u_d).$$

$$C(\cdot |{\theta }_k):=\prod _{i=1}^m{\phi }_i^{{\theta }_{ik}}(\cdot ),$$

$${\tilde{C}}^{(K)}(\cdot )=C(\cdot |{\tilde{\theta }}_K),$$

$${\tilde{\theta }}_{iK}=\sum _{k=1}^K{\theta }_{ik}(p^{(k)}{)}^{{\lambda }_i}.$$

$$C_k(\mathbf{u})=C(\mathbf{u}|{\theta }_k)=\prod _{j=1}^d{u}_j\exp \left({\theta }_k\prod _{j=1}^d\ln (1/u_j)\right)={\phi }_1^{{\theta }_{1k}}(\mathbf{u}){\phi }_2^{{\theta }_{2k}}(\mathbf{u})$$

$${\phi }_1(\mathbf{u})=\prod _{j=1}^d{u}_j\text{mbox}and{\phi }_2(\mathbf{u})=\exp \left(\prod _{j=1}^d\ln (1/u_j)\right).$$

$$C_{\#}(u_1,\dots ,u_d)=\exp (-L(\ln (1/u_1),\dots ,\ln (1/u_d))),$$

$${\tilde{\theta }}_{1K}=\sum _{k=1}^K{p}^{(k)}=1$$

$${\tilde{C}}^{(K)}(\mathbf{u}):=\tilde{C}(\mathbf{u})=C^K(u_1^{1/K},\dots ,u_d^{1/K}),$$

$${\tilde{C}}^{(1)}(\mathbf{u})$$

$${\tilde{C}}^{(k)}(\mathbf{u})$$

$$\underset{k=i}{\overset{j}{⨀}}{f}^{(k)}:=f^{(i)}\circ \dots \circ f^{(j)}\text{mbox}ifi\le j\text{mbox}and\underset{k=i}{\overset{j}{⨀}}{f}^{(k)}:=\mathrm{Id}\text{mbox}otherwise.$$

$${\tilde{C}}^{(K)}(\mathbf{u})=\prod _{k=1}^K{C}_k\left(g_1^{(K-k+1,K)}(u_1),\dots ,g_d^{(K-k+1,K)}(u_d)\right)$$

Full text

††journal: Journal of Multivariate Analysis

1 Introduction

Let $\bm{X}=(X_{1},\ldots,X_{d})$ be a continuous random vector with $d$-variate cumulative distribution function (cdf) $F$, and let $F_{j}$, $j\in\{1,\ldots,d\}$, be the marginal cdf of $X_{j}$. According to Sklar’s theorem [36], there exists a unique $d$-variate function $C:[0,1]^{d}\rightarrow[0,1]$ such that

[TABLE]

The function $C$ is referred to as the copula associated with $F$. It is the $d$-dimensional cdf of the random vector $(F_{1}(X_{1}),\ldots,F_{d}(X_{d}))$ with uniform margins on $[0,1]$.

A copula is said to be symmetric (or exchangeable) if for any $\bm{u}\in[0,1]^{d}$, and for any permutation $(\sigma_{1},\ldots,\sigma_{d})$ of the first $d$ integers $\{1,\ldots,d\}$, it holds that $C(u_{1},\ldots,u_{d})=C(u_{\sigma_{1}},\ldots,u_{\sigma_{d}}).$ The assumption of exchangeability may be unrealistic in many domains, including quantitative risk management [7], reliability modeling [39], and oceanography [41]. The urge for asymmetric copula models in order to better account for complex dependence structures has recently stimulated research in several directions, including Rodrıguez-Lallena and Úbeda-Flores [32], Alfonsi and Brigo [1], Durante [8], Wu [39], Durante et al. [9]. We focus here on a simple yet general method for building asymmetric copulas introduced by Liebscher [22, Theorem 2.1, Property (i)]:

Theorem 1**.**

[Liebscher, 22]*

Let $C_{1},\ldots,C_{K}:[0,1]^{d}\to[0,1]$ be copulas, $g_{j}^{(k)}:[0,1]\to[0,1]$ be increasing functions such that $g_{j}^{(k)}(0)=0$ and $g_{j}^{(k)}(1)=1$ for all $k\in\{1,\ldots,K\}$ and $j\in\{1,\ldots,d\}$. Then,*

[TABLE]

is also a copula under the constraint that

[TABLE]

Theorem 1 provides a generic way to construct an asymmetric copula $\tilde{C}$, henceforth referred to as Liebscher copula, starting from a sequence of symmetric copulas $C_{1},\ldots,C_{K}$. This mechanism was first introduced by [19] in the particular case where $K=2$ and with the functions $g_{j}^{(k)}$ assumed to be power functions, for each $j\in\{1,\ldots,d\}$ and $k\in\{1,\ldots,K\}$, that satisfy condition (2). The class of Liebscher copulas covers a broad range of dependencies and benefits from tractable bounds on dependence coefficients of the bivariate marginals [22, 23, 27]. However, there are two main reasons why the practical implementation of this approach is not straightforward: (i) it is not immediate to construct functions that satisfy condition (2); and (ii) the product form complicates the density computation even numerically, which makes it difficult to perform likelihood inference on the model parameters [27].

The aim of this paper is to deepen the understanding of Liebscher’s construction in order to overcome drawbacks (i) and (ii). Our contributions in this regard are three-fold. First, we provide theoretical properties of the asymmetric copulas in (1), including exact expressions of tail dependence indices, thus complementing the partial results of Liebscher [22, 23]. Second, we give an iterative representation of (1) which has the advantage to relax assumption (2) by automatically satisfying it. Third, we develop an inferential procedure and a sampling scheme that rely on the newly developed iterative representation.

The Bayesian paradigm proves very useful for inference in our context as it overcomes the problematic computation of the maximum likelihood estimate, which requires the maximization of a very complicated likelihood function (see recent contributions [38, 29]). General Bayesian sampling solutions in the form of Markov chain Monte Carlo are not particularly well-suited neither since they require the evaluation of that complex likelihood. Instead, we resort to Approximate Bayesian computation (ABC), a technique dedicated to models with complicated, or intractable, likelihoods (see [25, 31, 18] for recent reviews). ABC requires the ability to sample from the model, which is straightforward with our iterative representation of Liebscher copula. The adequacy of ABC for inference in copula models was leveraged by Grazian and Liseo [13], although in the different context of empirical likelihood estimation. A reversed approach to ours is followed by Li et al. [21], who make use of copulas in order to adapt ABC to high-dimensional settings.

Since its introduction, the construction by Liebscher has received much attention in the copula literature (e.g., [34, 10, 20]). However, most studies have been limited to simple cases where the product in (1) has only two terms. We hope that our paper will contribute to the further spreading of Liebscher’s copulas, because it allows to exploit their full potential by: (i) better understanding their properties; (ii) providing a novel construction, which facilitates their use with an arbitrary number $K$ of terms in (1); and (iii) giving a strategy to make inference on them.

On top of what has been presented above, an additional contribution of this paper is to derive specific results for the subclass of Liebscher’s copula when two or more comonotonic copulas are combined, which we call comonotonic-based Liebscher copula. This subclass is characterized by an arbitrary number of singular components. To the best of our knowledge, this is the first paper to investigate this copula’s properties and to provide an inference procedure.

The paper is organized as follows. Section 2 provides some theoretical results concerning the properties of asymmetric Liebscher copulas, also presenting the novel iterative construction. In Section 3, we introduce and analyze the comonotonic-based Liebscher copula. Section 4 is dedicated to the inference strategy. It demonstrates our approach on simulated data and provides a comparison with a likelihood-based approach for a class of Liebscher copulas where maximum likelihood estimation is feasible. We conclude with a short discussion in Section 5. Proofs are postponed to the Appendix.

2 Properties of the copula

In this section, some new properties of the copula $\tilde{C}$ are established, complementing the ones in Liebscher [22, 23]. Sections 2.1 and 2.2 are dedicated to (tail) dependence properties. For the sake of simplicity, we focus on the case $d=2$ of bivariate copulas. Some stability properties of Liebscher’s construction are highlighted in Section 2.3. Finally, an alternative construction to Liebscher copula (1) is introduced in Section 2.4.

2.1 Tail dependence

The lower and upper tail dependence functions, denoted by $\Lambda_{L}(C;\cdot)$ and $\Lambda_{U}(C,\cdot)$ respectively, are defined for all $(x,y)\in\mathbb{R}_{+}^{2}$ by

[TABLE]

where $C$ is a given bivariate copula, see for instance [17]. Note that these limits exist under a bivariate regular variation assumption, see [30], Section 5.4.2 for details. When they exist, these functions are homogeneous ([17], Proposition 2.2), i.e., for all $t\in(0,1]$ and $(x,y)\in\mathbb{R}_{+}^{2}$, $\Lambda(C;tx,ty)=t\Lambda(C;x,y)$, where $\Lambda$ is equal to $\Lambda_{L}$ or $\Lambda_{U}$. The lower and upper tail dependence coefficients, denoted by $\lambda_{L}(C)$ and $\lambda_{U}(C)$ respectively, are defined as the conditional probabilities that a random vector associated with a copula $C$ belongs to lower or upper tail orthants given that a univariate margin takes extreme values:

[TABLE]

These coefficients can also be interpreted in terms of the tail dependence functions: $\lambda_{L}(C)=\Lambda_{L}(C;1,1)$ and $\lambda_{U}(C)=\Lambda_{U}(C;1,1)$. Conversely, in view of the homogeneity property, the behavior of the tail dependence functions on the diagonal is determined by the tail dependence coefficients: $\Lambda_{L}(C;t,t)=\lambda_{L}(C)t$ and $\Lambda_{U}(C;t,t)=\lambda_{U}(C)t$ for all $t\in(0,1]$. The tail dependence functions for Liebscher copula are provided by Proposition 1 which, in view of the previous remarks, allows us to derive the tail dependence coefficients in Corollary 1. Some of these results rely on the notion of (univariate) regular variation. Recall that a positive function $g$ is said to be regularly varying with index $\gamma$ if $g(xt)/g(x)\to t^{\gamma}$ as $x\to\infty$ for all $t>0$, see [4].

Proposition 1**.**

Let $(x,y)\in\mathbb{R}_{+}^{2}$ and consider $\tilde{C}$ the bivariate copula defined by (1) with $d=2$.

(i)* Lower tail, symmetric case. Assume that $g^{(k)}_{1}=g^{(k)}_{2}$ is a regularly varying function with index $\gamma^{(k)}>0$ for all $k\in\{1,\dots,K\}$. Then,*

[TABLE]

and, necessarily, $\sum_{k=1}^{K}\gamma^{(k)}=1$.

(ii)* Lower tail, asymmetric case. Suppose there exists $k_{0}\in\{1,\dots,K\}$ such that $g_{1}^{(k_{0})}(\varepsilon)/g_{2}^{(k_{0})}(\varepsilon)\to 0$ as $\varepsilon\to 0$. Then,*

[TABLE]

(iii)* Upper tail, general case. Assume that, for all $k\in\{1,\dots,K\}$, $g^{(k)}_{1}$ and $g^{(k)}_{2}$ are differentiable at 1, with derivative at 1 denoted by $\text{\emph{d}}_{1}^{(k)}$ and $\text{\emph{d}}_{2}^{(k)}$ respectively. Then*

[TABLE]

and, necessarily, $\sum_{k=1}^{K}\text{\emph{d}}_{j}^{(k)}=1$, for $j\in\{1,2\}$.

(iv)* Upper tail, particular case. If, in addition to *(iii), $\text{\emph{d}}_{1}^{(k)}=\text{\emph{d}}_{2}^{(k)}\eqqcolon\text{\emph{d}}^{(k)}$ for all $k\in\{1,\dots,K\}$, then

[TABLE]

and, necessarily, $\sum_{k=1}^{K}\text{\emph{d}}^{(k)}=1$.

Let us note that the functions $g_{j}^{(k)}$ considered by Liebscher [22] and indexed by (I-III) in his Section 2.1 all satisfy the assumptions of Proposition 1. The following result complements [22, Proposition 2.3] and [23, Proposition 0.1] which provide bounds on the tail dependence coefficients. Here instead, explicit calculations are provided.

Corollary 1**.**

Let $\tilde{C}$ be the bivariate copula defined by (1) with $d=2$.

(i)* Lower tail, symmetric case. Under the assumptions of Proposition *1(i), $\lambda_{L}(\tilde{C})=\prod_{k=1}^{K}\lambda_{L}(C_{k}).$

(ii) Lower tail, asymmetric case. Under the assumptions of Proposition 1(ii), $\lambda_{L}(\tilde{C})=0.$

(iii) Upper tail, general case. Under the assumptions of Proposition 1(iii), $\lambda_{U}(\tilde{C})=\sum_{k=1}^{K}\Lambda_{U}(C_{k};\text{\emph{d}}_{1}^{(k)},\text{\emph{d}}_{2}^{(k)})$ and, necessarily, $\sum_{k=1}^{K}\text{\emph{d}}_{j}^{(k)}=1$, for $j\in\{1,2\}$.

(iv) Upper tail, particular case. Under the assumptions of Proposition 1(iv), $\lambda_{U}(\tilde{C})=\sum_{k=1}^{K}\text{\emph{d}}^{(k)}\lambda_{U}(C_{k})$ and, necessarily, $\sum_{k=1}^{K}\text{\emph{d}}^{(k)}=1$.

It appears that the lower and upper tail dependence coefficients have very different behaviors. In the case (i) of a symmetric copulas, the lower tail dependence coefficient $\lambda_{L}$ is the product of the lower tail dependence coefficients associated with the components. Besides, $\lambda_{L}=0$ as soon as a component $k_{0}$ has functions $g_{1}^{(k_{0})}$ and $g_{2}^{(k_{0})}$ with different behaviors at the origin (case (ii)). At the opposite, the upper tail dependence coefficient does not vanish even though a component $k_{0}$ has functions $g_{1}^{(k_{0})}$ and $g_{2}^{(k_{0})}$ with different behaviors at 1 (case (iii)). In the particular situation where all components $k\in\{1,\ldots,K\}$ have functions $g_{1}^{(k)}$ and $g_{2}^{(k)}$ with the same behavior at 1 (case (iv)), $\lambda_{U}$ is a convex combination of the upper tail dependence coefficients associated with the components.

2.2 Dependence

Let $(X,Y)$ be a pair of random variables with continuous margins and associated copula $C$.

$X$ and $Y$ are said to be totally positive of order 2, $\mathrm{TP}_{2}$ (see [16]), if for all $x_{1}<y_{1},\,x_{2}<y_{2}$,

[TABLE]

Since this can be equivalently written in terms of $C$, we will write in short that $C$ is $\mathrm{TP}_{2}$.

• 2.

$X$ and $Y$ are said to be positively quadrant dependent (PQD) if

[TABLE]

Since this property can be characterized by the copula property $C\geq\Pi$ where $\Pi$ denotes the independence copula, see for instance [28, Paragraph 5.2.1], we shall write for short that $C$ is PQD. The negatively quadrant dependence (NQD) property is similarly defined by $C\leq\Pi$.

• 3.

$X$ and $Y$ are said to be left-tail decreasing (LTD) if

[TABLE]

From Nelsen [28, Theorem 5.2.5], this property can be characterized by the copula properties

[TABLE]

and we shall thus write that $C$ is LTD. The left-tail increasing property (LTI) is similarly defined by reversing the directions of the monotonicity in (LABEL:eq:LTI-def) and (4).

• 4.

$X$ and $Y$ are said to be stochastically increasing (SI) if

[TABLE]

From Nelsen [28, Corollary 5.2.11], this property can be characterized by the copula properties

[TABLE]

and we shall thus write that $C$ is SI. The stochastically decreasing (SD) property is similarly defined by replacing increasing by decreasing in (LABEL:eq:SI-def) and concave by convex in (6).

In the next proposition, we show that under mild conditions, the above dependence properties are preserved under Liebscher’s construction, thus complementing $\mathrm{LTD}$ and $\mathrm{TP}_{2}$ properties established in Liebscher [22, Proposition 2.2].

Proposition 2**.**

If copulas $C_{1},\ldots,C_{K}$ all satisfy any of the properties defined above, $\mathrm{TP}_{2}$, $\mathrm{PQD}$, $\mathrm{NQD}$, $\mathrm{LTD}$, $\mathrm{LTI}$, $\mathrm{SI}$ or $\mathrm{SD}$, then the same is satisfied for the Liebscher copula $\tilde{C}$ defined in (1)—for $\mathrm{SI}$ (respectively $\mathrm{SD}$), the $g_{j}^{(k)}$ functions in Theorem 1 are additionally required to be concave functions (respectively convex functions) and the copulas $C_{k}$ to be twice differentiable, $k\in\{1,\ldots,K\}$, $j\in\{1,\ldots,d\}$.

2.3 Stability properties

Let us focus on the situation where the functions $g_{j}^{(k)}$ of Theorem 1 are power functions: for all $j\in\{1,\dots,d\}$, $k\in\{1,\dots,K\}$ and $t\in[0,1]$, let

[TABLE]

Recall that a copula $C_{\#}$ is said to be max-stable if for all integer $n\geq 1$ and $(u_{1},\dots,u_{d})\in[0,1]^{d}$:

[TABLE]

From [11, Proposition 3], it is clear that associating max-stable copulas $C_{k}$ with power functions (7) in Liebscher construction (1) still yields a max-stable copula. The goal of this paragraph is to investigate to what extent this result can be generalized. Our first result establishes the stability of a family of Liebscher copulas built from homogeneous functions. More specifically, each copula $C_{k}(\cdot)$ in (1) is rewritten as $C(\cdot\;|\;\theta_{k})$ where

[TABLE]

with $\theta_{k}=(\theta_{1k},\dots,\theta_{mk})^{T}$ and $\varphi_{i}:[0,1]^{d}\to[0,1]$, $i\in\{1,\dots,m\}$.

Proposition 3**.**

For all $j\in\{1,\dots,d\}$ and $t\in[0,1]$, let $g_{j}^{(k)}$ be given by (7) where $p_{j}^{(k)}=p^{(k)}$ for all $k\in\{1,\dots,K\}$. Let $m>0$ and for all $i\in\{1,\dots,m\}$ introduce $\varphi_{i}:[0,1]^{d}\to[0,1]$ such that $\ln\mathop{\circ}\varphi_{i}\mathop{\circ}\exp$ is homogeneous of degree $\lambda_{i}$. For all $k\in\{1,\ldots,K\}$, assume that $C(\cdot\;|\;\theta_{k})$ in (8) is a copula for some $\theta_{k}\in{\mathbb{R}}^{m}$. Then, copula (1) is given for all $K\geq 1$ by

[TABLE]

with $\tilde{\theta}_{K}=(\tilde{\theta}_{1K},\dots,\tilde{\theta}_{mK})^{T}$, and for all $i\in\{1,\dots,m\}$,

[TABLE]

Example 1** (Gumbel-Barnett copula $C_{k}$).**

Let $C_{k}$ be the Gumbel-Barnett copula [28, Table 4.1]. It can be written as

[TABLE]

with $\theta_{1k}=1$, $\theta_{2k}=\theta_{k}\geq 1$,

[TABLE]

It thus fulfills the assumptions of Proposition 3 with $m=2$, $\lambda_{1}=1$ and $\lambda_{2}=d$.

Example 2** (Extreme-value copula $C_{k}$).**

Extreme-value copulas exactly correspond to max-stable copulas and are characterized by their tail-dependence function $L$ as:

[TABLE]

where $L:{\mathbb{R}}_{+}^{d}\to{\mathbb{R}}_{+}$ is homogeneous of degree 1, see for instance [14]. It is thus clear that every max-stable copula $C_{k}$ fulfills the assumptions of Proposition 3 with $m=1$, $\theta_{1k}=1$, $\lambda_{1}=1$ and $\ln\circ\varphi_{1}\circ\exp({\bm{t}})=-L(-{\bm{t}})$, $t\in{\mathbb{R}}_{-}^{d}$. Moreover,

[TABLE]

and thus $\tilde{C}^{(K)}=C$ for all $K\geq 1$.

It appears that max-stable copulas can be considered as fixed-points of Liebscher’s construction (1). The next result shows that, under mild assumptions, they are the only copulas verifying this property.

Proposition 4**.**

For all $j\in\{1,\dots,d\}$, let $g_{j}^{(k)}$ be given by (7) where $p_{j}^{(k)}=p^{(k)}$ for all $k\in\{1,\dots,K\}$. Assume $C_{k}=C$ for all $k\in\{1,\dots,K\}$ and let $\tilde{C}^{(K)}$ be the copula defined by (1). Then, $\tilde{C}^{(K)}=C$ for all $K\geq 1$ and for all sequences $p^{(1)},\ldots,p^{(K)}\in(0,1)$ such that $\sum_{k=1}^{K}p^{(k)}=1$ if and only if $C$ is max-stable.

To complete the links with max-stable copulas, let us consider the situation where $C_{k}=C$ and $p_{j}^{(k)}=1/K$ for all $j\in\{1,\dots,d\}$ and $k\in\{1,\dots,K\}$. Liebscher’s construction thus yields

[TABLE]

which is the normalized cdf associated with the maximum of $K$ independent uniform random vectors distributed according to the cdf $C$. Therefore, as $K\to\infty$, $\tilde{C}^{(K)}$ converges to a max-stable copula under standard extreme-value assumptions on $C$.

2.4 An iterative construction

Let ${\cal F}$ be the class of increasing functions $f:[0,1]\to[0,1]$ such that $f(0)=0$, $f(1)=1$ and $\mathrm{Id}/f$ is increasing, where $\mathrm{Id}$ denotes the identity function. For all $k\geq 1$ let $C_{k}$ be a $d$-variate copula and $f^{(k)}_{j}\in{\cal F}$ for all $j\in\{1,\dots,d\}$, with the assumption $f^{(1)}_{j}(t)=1$ for all $t\in[0,1]$. We propose the following iterative construction of copulas. For all $\bm{u}\in[0,1]^{d}$, consider the sequence defined by

[TABLE]

Lemma 2 in B shows that $\tilde{C}^{(k)}$ is a $d$-variate copula, for all $k\geq 1$.

Let $K\geq 1$. For all functions $f^{(1)},\dots,f^{(K)}$ : $[0,1]\to[0,1]$ and $i,j\in\{1,\dots,K\},$ let us introduce the notation

[TABLE]

The next result shows that there is a one-to-one correspondence between copulas built by the iterative procedure (9), (10) and Liebscher copulas, reported in Theorem 1.

Proposition 5**.**

The copula $\tilde{C}^{(K)}$, $K\geq 1$ defined iteratively by (9), (10) is a Liebscher copula. It can be rewritten as

[TABLE]

for all $\bm{u}\in[0,1]^{d}$ where, for all $j\in\{1,\dots,d\}$ and $K\geq 1$,

[TABLE]

Conversely, each Liebscher copula (defined in Theorem 1) can be built iteratively from (9), (10).

Let us note that the iterative construction (9), (10) thus provides a way to build functions (13), (14) that automatically fulfill Liebscher’s constraints (2) of Theorem 1. As a consequence, the construction (9), (10) also gives an iterative way to sample from a Liebscher copula (1), described in detailed in Algorithm 1 (see the proof of Lemma 2 in B for a theoretical justification).

Example 3** (Power functions $f_{j}^{(k)}$).**

Let functions $f_{j}^{(k)}$ be power functions in the form of

[TABLE]

for all $j\in\{1,\dots,d\}$ and $t\in[0,1]$. From Proposition 5, $g_{j}^{(k)}(t)=g_{j}^{(k,K)}(t)=t^{p_{j}^{(k,K)}}$ with

[TABLE]

for all $K\geq 1$ and $j\in\{1,\dots,d\}$. Note that, by construction,

[TABLE]

for all $j\in\{1,\dots,d\}$ and thus $\left(p_{j}^{(1,K)},\dots,p_{j}^{(K,K)}\right)$ can be interpreted as a discrete probability distribution on $\{1,\dots,K\}$. Besides, let $\tilde{a}_{j}^{(k,K)}\coloneqq a_{j}^{(K+1-k)}$ for all $k\in\{1,\dots,K\}$. Equations (16) can be rewritten as

[TABLE]

which corresponds to the so-called stick-breaking construction [35].

3 The comonotonic-based Liebscher copula

We analyze here in more details the Liebscher copula obtained by combining $K\geq 2$ comonotonic (also called Fréchet) copulas defined by ${C}(u,v)=\min(u,v)$. We here focus on the bivariate case $d=2$, although some of the derivations carry over to the general $d$-dimensional case. We consider the specific case of Example 3, where the functions in Liebscher’s construction are power functions, $g_{j}^{(k)}(t)=t^{p_{j}^{(k,K)}}$ with $p_{j}^{(k,K)}\in[0,1]$, $j\in\{1,2\}$ as in (16). Assuming that $K$ is fixed and limiting ourselves to $d=2$, we denote for notation simplicity $p_{k}\coloneqq p_{1}^{(k,K)}$ and $q_{k}\coloneqq p_{2}^{(k,K)}$ for $k\in\{1,\dots,K\}$. Recall that, in view of (2), $\sum_{k=1}^{K}p_{k}=\sum_{k=1}^{K}q_{k}=1$. Under the above assumptions, the comonotonic-based Liebscher copula denoted by $\tilde{C}_{\mathrm{CL}}$ has the form

[TABLE]

In the particular case where $K=2$, it is referred to as the BC2 copula by [24] and it is proved that any bivariate extreme-value copula with arbitrary discrete dependence measure can be represented as the geometric mean of BC2 copulas, which corresponds to the situation where $K$ is even in (17). We also refer to [37] for further links with extreme-value theory.

3.1 Geometric description of $\tilde{C}_{\mathrm{CL}}$

For all $k\in\{1,\dots,K\}$, introduce $r_{k}=p_{k}/q_{k}\in[0,\infty]$. For notation simplicity, we shall let $r_{0}=0$ and $r_{K+1}=\infty$. Since the above product (17) is commutative, one can assume without loss of generality that the sequence $(r_{k})_{0\leq k\leq K+1}$ is nondecreasing. The copula $\tilde{C}_{\mathrm{CL}}$ can be easily expressed on the partition of the unit square $[0,1]^{2}$ defined by the following moon shaped subsets (see the illustration in Fig. 1 with $K=2$)

[TABLE]

Proposition 6**.**

Let $\tilde{C}_{\mathrm{CL}}$ be the comonotonic-based Liebscher copula defined in (17). Then, for all $(u,v)\in[0,1]^{2}$,

[TABLE]

where $\bar{x}_{k}=x_{1}+\ldots+x_{k}$, with the convention that $\bar{x}_{0}=0$. Moreover, the singular component of $\tilde{C}_{\mathrm{CL}}$ is

[TABLE]

The singular component $\tilde{S}_{\mathrm{CL}}$ and the absolute continuous component $\tilde{A}_{\mathrm{CL}}=\tilde{C}_{\mathrm{CL}}-\tilde{S}_{\mathrm{CL}}$ weights are $\sum_{k=1}^{K}\min(p_{k},q_{k})$ and $1-\sum_{k=1}^{K}\min(p_{k},q_{k})$, respectively.

A key property of the comonotonic-based Liebscher copula (17) is the presence of multiple singular components lying on the curves $v=u^{r_{k}}$ with associated weights $\min(p_{k},q_{k})$, $k\in\{1,\dots,K\}$. As an illustrative example, let us consider the bivariate comonotonic-based Liebscher copula defined with $p_{1}=1-p_{2}=1/3$, $q_{1}=1-q_{2}=3/4$,

[TABLE]

which entails $r_{0}=0,r_{1}=4/9,r_{2}=8/3$ and $r_{3}=\infty$. The moon shaped subsets of the partition of the unit square $[0,1]^{2}$ are represented on Fig. 1, and the expressions of $\tilde{C}_{\mathrm{CL}}$ and the singular component $\tilde{S}_{\mathrm{CL}}$ are as follows:

[TABLE]

Appropriate choices of pairs $(p_{k},q_{k})_{k=1,\dots,K}$ may lead to a number of singular components ranging from [math] to $K$. The independence copula (no singular component) is obtained for instance with $p_{i}=b_{j}=1$ for a given pair $(i,j)$, $i\neq j$, the comonotonic copula (one singular component) is obtained by choosing $p_{k}=q_{k},\,\forall k\in\{1,\ldots,K\}$, and a copula with exactly $K$ singular components can be obtained provided that $0<r_{1}<r_{2}<\dots<r_{K}<\infty$. See Fig. 2 and Fig. 3 for illustrations. As a comparison, Cuadras and Augé [6] copula given by $C(u,v)=(uv)^{1-\theta}\min(u,v)^{\theta}$, $\theta\in[0,1]$ is limited to a single singular component, necessarily on the diagonal $v=u$. Similarly, Marshall and Olkin [26] copula is defined by $C(u,v)=\min(u^{1-\alpha}v,uv^{1-\beta})$, $(\alpha,\beta)\in[0,1]^{2}$ and has only one singular component located on the curve $v=u^{\alpha/\beta}$. The proposal by Lauterbach and Pfeifer [20] based on singular mixture copulas includes comonotonic-based Liebscher copula (17) in the particular case when $K=2$ but is limited to two singular components. Finally, Sibuya copulas [15] is a very general family of copulas: Let us point out that, in the bivariate case, a non-homogeneous Poisson Sibuya copula allows for only one singular component, this singular component being supported by a curve with very flexible shape (see Remark 4.2 in the previously referenced work for further details).

3.2 Dependence and association properties of $\tilde{C}_{\mathrm{CL}}$

We consider here several measures of dependence and measures of association between the components of the bivariate comonotonic-based Liebscher copula (17). Some of these measures are already dealt with in great generality in Section 2, while some others seem to be tractable only in the comonotonic-based Liebscher copula case: Blomqvist’s medial correlation coefficient, Kendall’s $\tau$ and Spearman’s $\rho$.

Tail dependence

Recall that for the comonotonic copula $C_{\text{C}}$, it holds $\Lambda_{U}(C_{\text{C}};\cdot,\cdot)=\min(\cdot,\cdot)$. Then, Corollary 1 yields

[TABLE]

In other words, the lower tail dependence coefficient is non zero only in the case when $p_{k}=q_{k}$ for all $k\in\{1,\ldots,K\}$, where $\tilde{C}_{\mathrm{CL}}$ boils down to the comonotonic copula, while the upper tail dependence coefficient coincides with the weight of the singular component $\tilde{S}_{\mathrm{CL}}$ (see Proposition 6). These results were also established in [24, Lemma 3], in the particular case where $K=2$.

Dependence

It is well-known that the comonotonic copula $C_{\mathrm{F}}$ fulfills the following positive dependence properties defined in Section 2.2, namely it is $\mathrm{TP}_{2}$, $\mathrm{PQD}$, $\mathrm{LTD}$ and $\mathrm{SI}$ [28]. According to Proposition 2, we can thus conclude that the comonotonic-based Liebscher copula $\tilde{C}_{\mathrm{CL}}$ also satisfies these positive dependence properties.

Stability properties

It is easily seen that comonotonic-based Liebscher copula (17) is max-stable. Proposition 4 thus entails that this copula is stable with respect to Liebscher’s construction. Another consequence is that comonotonic-based Liebscher construction (17) can be interpreted as a possible cdf for modelling bivariate maxima.

Dependence coefficients

The $\beta$-Blomqvist’s medial correlation coefficient ([28], Paragraph 5.1.4) defined by \beta(C)=4C\big{(}\frac{1}{2},\frac{1}{2}\big{)}-1, Kendall’s $\tau$ ([28], Paragraph 5.1.1) and Spearman’s $\rho$ ([28], Paragraph 5.1.2) defined by

[TABLE]

are provided in next proposition.

Proposition 7**.**

Blomqvist’s medial correlation coefficient, Kendall’s $\tau$ and Spearman’s $\rho$ for the comonotonic-based Liebscher copula (17) are respectively given by

[TABLE]

where $\bar{x}_{k}=x_{1}+\ldots+x_{k}$.

It appears that $\beta$-Blomqvist medial correlation coefficient is closely related to the upper tail dependence coefficient: $\beta(\tilde{C}_{\mathrm{CL}})=2^{\lambda_{U}(\tilde{C}_{\mathrm{CL}})}-1$. Besides, in the particular case where $K=2$, these results coincide the ones of Lemma 2 of [24]: the Kendall’s $\tau$ can be simplified as $\tau(\tilde{C}_{\mathrm{CL}})=p_{1}+q_{2}=\lambda_{U}(\tilde{C}_{\mathrm{CL}})$. No similar simplification seems to be possible for Spearman’s $\rho$.

For the special case of copula (21), we have $\lambda_{L}(\tilde{C}_{\mathrm{CL}})=0$, $\lambda_{U}(\tilde{C}_{\mathrm{CL}})=\tau(\tilde{C}_{\mathrm{CL}})=p_{1}+q_{2}=\frac{1}{3}+\frac{1}{4}=\frac{7}{12}\approx 0.583$, $\beta(\tilde{C}_{\mathrm{CL}})=2^{7/12}-1\approx 0.498$ and $\rho(\tilde{C}_{\mathrm{CL}})\approx 0.298$.

3.3 Iterative construction for $\tilde{C}_{\mathrm{CL}}$

Algorithm 1 can be simplified when specified to comonotonic-based Liebscher setting since: (i) sampling from the comonotonic copula is straightforward and only requires sampling from the uniform distribution $\mathcal{U}(0,1)$, and (ii) power functions benefit from an explicit inverse. The specific sampling procedure for this construction is described in detail as Algorithm 2.

4 Bayesian inference

In this section, we provide a simple strategy to make Bayesian inference on any Liebscher copula based on an Approximate Bayesian computation algorithm (ABC, see for instance [25, 31, 18] for reviews). ABC is a “likelihood-free” method usually employed for inference of models with intractable likelihood: it enables to perform approximate Bayesian analysis on any statistical model from which it is possible to sample new data, without the need to explicitly evaluate the likelihood function.

Let $\bm{X}_{\mathrm{obs}}=\{\bm{X}_{\mathrm{obs},1},\ldots,\bm{X}_{\mathrm{obs},n}\}$ be the observed data, where $\bm{X}_{\mathrm{obs},i}=(X_{\mathrm{obs},i,1},\ldots,X_{\mathrm{obs},i,d})$, $i\in\{1,\ldots,n\}$, and assume that the statistical model for $\bm{X}_{\mathrm{obs}}$ is described by a likelihood function $\mathcal{L}_{\theta}$ with parameter $\theta$ which is to be inferred. The basic scheme of one step of ABC is the following:

Sample $\theta$ from the prior distribution $\pi(\theta)$; 2. 2.

Given $\theta$, sample $\bm{X}_{1},\ldots,\bm{X}_{n}$ from $\mathcal{L}_{\theta}$, and set $\bm{X}=\{\bm{X}_{1},\ldots,\bm{X}_{n}\}$; 3. 3.

If $\bm{X}$ is too different from $\bm{X}_{\mathrm{obs}}$, discard $\theta$, otherwise, keep $\theta$.

The outcome of the ABC algorithm is a sample of values of the parameter $\theta$ approximately distributed according to its posterior distribution. The basic (rejection) ABC approach in point 3. amounts to a priori specifying a tolerance level $\epsilon>0$, and then keeping $\theta$ if $d(\bm{X},\bm{X}_{\mathrm{obs}})<\epsilon$ for some distance $d(\,\cdot\,,\,\cdot\,)$ between samples. Another common approach employed in this paper consists in selecting the tolerance level $\epsilon$ as a fixed quantile of the distances $d(\bm{X},\bm{X}_{\mathrm{obs}})$. More specifically, Steps 1. to 3. are repeated $M^{\prime}$ times, out of which $M$ are retained, yielding a quantile of order $M/M^{\prime}$ [31]. In other words, the $M$ retained parameters are those associated with the smallest values of the distance $d(\bm{X},\bm{X}_{\mathrm{obs}})$. In this paper, we choose as distance between samples the Hilbert distance introduced by [3], henceforth denoted by $d_{H}(\,\cdot\,,\,\cdot\,)$. The Hilbert distance is an approximation of the Wasserstein distance between empirical probability distributions which preserves the desirable properties of the latter in the context of ABC, while being considerably faster to compute in multivariate data settings. More precisely, given two samples, $y_{1:n}$ and $z_{1:n}$, the Hilbert distance of order 1 associated with the Euclidean distance (henceforth simply referred to as the Hilbert distance) is defined as $d_{H}(y_{1:n},z_{1:n})=\frac{1}{n}\sum_{i=1}^{n}||y_{i}-z_{\sigma(i)}||$, where $\sigma(i)=\sigma_{z}\circ\sigma^{-1}_{y}(i)$ for all $i=1,\dots,n$, and $\sigma_{y}$ and $\sigma_{z}$ are the permutations obtained by mapping the vectors $y_{1:n}$ and $z_{1:n}$ through their projection via the Hilbert space-filling curve [33] and sorting the resulting vectors in increasing order (see [3] for details).

The choice for ABC is motivated by two main reasons. First, it is nontrivial in general to derive the likelihood of copulas, especially for Liebscher copulas which involve differentiating a product of $K$ terms. All the more, the specific case of comonotonic-based Liebscher copulas induces up to $K$ singular components which precludes a general evaluation of the likelihood. Second, sampling new data (step 2. above) from a Liebscher copula is straightforward and fast thanks to the iterative procedure of Algorithm 1 (Section 2.4).

Section 4.1 introduces the ABC procedure in the case of the Liebscher copula (1) and describes the prior distributions on the model parameters. The methodology is then illustrated on two data generating distributions: Section 4.2.1 focuses on the well-specified setting where the data are sampled from the comonotonic-based Liebscher copula; we show that the estimation procedure performs well in this case. Then, Section 4.2.2 investigates the misspecified setting where the data are sampled from a noisy version of the model; we show that the estimation procedure still performs well, but the estimation accuracy may deteriorate for too large values of the noise. Finally, Section 4.3 compares our proposed ABC approach to a likelihood-based technique.

4.1 ABC inference for Liebscher copulas

The description of the ABC procedure is first completed by specifying the prior distributions on the model parameters. For simplicity, we here focus on the case of the $d$-dimensional Liebscher copula with the functions $f_{j}^{(k)}(\cdot)$, $k\in\{1,\ldots,K\}$, $j\in\{1,\ldots,d\}$ of Algorithm 1 chosen as the power functions (15) introduced in Example 3. The parameters of the $f_{j}^{(k)}(\cdot)$ functions are collected in a $K\times d$ matrix ${A}=[a_{j}^{(k)}]_{k,j}$, where ${a}^{(1)}=(a_{1}^{(1)},\ldots,a_{d}^{(1)})=(1,\ldots,1)$. Since the $(K-1)d$ free parameters are constrained to $a_{j}^{(k)}\in(0,1)$ for $2\leq k\leq K$, we simply choose, by symmetry, independent and uniform distributions $a_{j}^{(k)}\stackrel{{\scriptstyle\mbox{\scriptsize{iid}}}}{{\sim}}\mathcal{U}(0,1)$. More flexible distributions like the Beta distribution could be thought of in order to reflect some prior knowledge on these parameters. Additionally, note that different functions $f_{j}^{(k)}(\cdot)$ would simply require setting prior distributions adapted to the parameters used.

The number of iterative steps $K$ is also considered as a parameter of the model. Independently of parameters $A$, $K$ is assigned a Zipf distribution, $K\sim\text{Zipf}(\xi)+1$, for $\xi>1$. Such a distribution is supported on integers $k\geq 2$ and has probability mass function $\mathbb{P}(K=k)$ proportional to $(k-1)^{-\xi}$. We further choose the parameter $\xi$ to be equal to 2, which insures that $90\%$ of the prior mass for $K$ is supported on most realistic values $2\leq K\leq 6$. This can be changed depending on applications at hand. Another option, useful in case where some prior information is available on $K$, is to adopt a Binomial distribution (translated, such that $K\geq 2$). The choice of the two hyper-parameters of the Binomial density could then be set as a function of prior knowledge, such as prior mode and confidence, that one may be able to elicit thanks to expert knowledge or previous studies.

We are now ready to state the main ABC inference procedure as Algorithm 3.

In general, the sequence of copulas $(C_{k})_{k}$ depends on some sequence of parameters $(\gamma_{k})_{k}$. In such a case, Algorithm 3 can be easily amended by adding a step consisting in sampling $\gamma_{k}$ parameters from some prior distribution to be set based on available prior information or expert knowledge. In the following section, we focus on the case of comonotonic-based Liebscher copulas, which do not depend on any additional parameter. Thus, the sampling step for new data $\bm{X}^{(s)}$ in Algorithm 3 is performed with the iterative construction of Algorithm 2 tailored to comonotonic-based Liebscher copulas.

4.2 Numerical illustrations with comonotonic-based Liebscher copulas

This section provides two illustrations of the inferential procedures described so far. The first investigates a setting where data are sampled from the comonotonic-based Liebscher model, while the second is concerned with observations from a noisy version of it. The code, implemented in R using the copula package [40] and winference package [2] for the Hilbert distance implementation [3], is available at the following link:

https://sites.google.com/site/crispinostat/research?authuser=0.

4.2.1 Well-specified setting: data from comonotonic-based Liebscher copula

We generate $n=500$ data points from a $2$-dimensional comonotonic-based Liebscher copula (17) with $d=2$, varying values of $K\geq 2$ and of the parameters in the matrix ${A}$, using Algorithm 2. The estimation is then performed with the ABC procedure summarized in Algorithm 3.

Our method provides the full (approximate) posterior distribution of the parameters of interest, making possible to select any strategy to summarize them, possibly driven by the application at hand. One can for instance compute the posterior distribution of the Spearman’s $\rho$, and, thanks to the retained samples, any other quantity of interest. Here, the performance of the estimation procedure is assessed basing on the following three summary statistics: (i) Kendall’s distribution function $\mathcal{K}(t)=\Pr_{C}({C}(U,V)<t)$, (ii) Spearman’s $\rho$ index of association, introduced in Section 3.2, and for which an explicit closed form is obtained for comonotonic-based Liebscher copula in Proposition 7, and (iii) an asymmetry measure, since it is a central motivation of the present work. More specifically, the Cramér-von Mises test statistics ${\rm E}_{C}[({C}(U,V)-{C}(V,U))^{2}]$ defined in Genest et al. [12] has been selected since it emerged as a powerful statistic to test the symmetry of a copula. Following the strategy of Genest et al. [12], the approximate p-values associated with the symmetry test, performed both on the observed sample $\bm{X}_{\text{obs}}$ and on the retained samples $\bm{X}_{m},m\in\{1,\ldots,M\}$, are computed on the basis of 250 bootstrap replicates.

The results obtained on a single simulation experiment are displayed on Fig. 4, where $n=500$ data points were simulated from the comonotonic-based Liebscher copula (17) with $d=2$ and $K=3$ (top-left panel). The number of ABC iterations was set to $M^{\prime}=10^{4}$, of which $M=300$ were retained (resulting in a quantile of order $3\%$). The empirical Kendall’s distribution functions of the observed and retained samples are compared on the top-right panel; The posterior distribution of $\rho$ is compared to the empirical Spearman’s $\rho$ of the observed sample on the bottom-left panel; Finally, the posterior distribution of the approximate p-values is displayed on the bottom-right panel. Let us highlight that the estimating procedure provides distributions around the true values in the three considered cases.

We then vary the generating number of iterative steps $K$ and compute the average relative errors $\eta_{\mathcal{K}}$ and $\eta_{\rho}$ for $\mathcal{K}$ and $\rho$ between the values computed on the observed sample and on the $M$ samples retained by ABC:

[TABLE]

where $\|\cdot\|_{1}$ denotes the $\ell_{1}$-norm. In order to take care of the randomness involved in sampling the parameters in the matrix $A$, the previous procedure has been replicated 20 times based on 20 independent data samples repetitions. The average relative errors $\eta_{\mathcal{K}}$ and $\eta_{\rho}$ in (22) were therefore averaged over the 20 independent samples, and reported as $\bar{\eta}_{\mathcal{K}}$ and $\bar{\eta}_{\rho}$ in Table 1 (first two rows), along with standard deviations in parentheses. As for the asymmetry test, we computed for each of the 20 data replications the fraction of times (out of $M$) that the same decision is taken (‘reject’ vs ‘do not reject’) at the 5% level, based on the approximate p-values computed on $\bm{X}_{m}$ and $\bm{X}_{\text{obs}}$. The obtained values were averaged over the 20 independent data samples repetitions and reported in the third row of Table 1 as $\bar{\text{f}}_{\text{test}}$.

Table 1 suggests a general trend: the larger $K$ is, the more difficult the estimation is. However, the estimation procedure yields satisfactory results for all cases considered.

4.2.2 Misspecified setting: data from a noisy comonotonic-based Liebscher copula

In this section, we generate data from a noisy version of the comonotonic-based Liebscher copula, and demonstrate that our inference procedure works well even if the data are not sampled from the exact model (so-called misspecified setting). In order to sample data from such a noisy model, a slightly changed version of Algorithm 2 is used in which the parameters ${A}$ are not fixed. Instead, they are sampled from a beta distribution with given variance $\sigma_{a}^{2}$ (interpreted as the error variance) around some fixed value corresponding to the zero noise version. The latter is illustrated on Fig. 5: a sample of $n=10^{3}$ data points from a 2-dimensional comonotonic-based Liebscher copula is depicted on the top-left panel with $K=2$ iterative steps, and with the two parameters of the power functions set to $a_{1}^{(2)}=0.4$ and $a_{2}^{(2)}=0.8$ (recall that $a_{1}^{(1)}=a_{2}^{(1)}=1$). The remaining five panels correspond to samples from comonotonic-based Liebscher copula with increasing noise variance $\sigma_{a}^{2}=10^{-5},10^{-4},10^{-3},10^{-2},10^{-1}$.

The inference results obtained from Algorithm 3 (ran with $M^{\prime}=10^{4}$ ABC iterations of which $M=300$ were retained) are reported in Table 2. Unsurprisingly, it appears that the larger the noise variance is, the more difficult is the estimation. Again, the estimation procedure yields good results for all cases considered. Let us note that the results reported in Table 2 are not averaged over 20 independent replications like in the previous section. The reason is that, here, the interest is in illustrating how the procedure deteriorates with the increasing noise in the observed data. Therefore, it is sufficient to illustrate the results of the analysis on a single dataset.

4.3 Comparison of ABC with likelihood-based estimation

In this section, the inference based on our ABC procedure is compared with the likelihood-based method provided in the copula R package [40].

As a matter of fact, in some specific cases it is possible to derive the density of Liebscher copulas, and hence, to perform likelihood-based inference on them. We work here with the following bivariate Liebscher copula obtained by combining Clayton copula $C_{1}(u,v)=\left(u^{-\theta}+v^{-\theta}-1\right)^{-1/\theta}$ with the independence copula $C_{2}(u,v)=uv$, and by using power functions (Example 3),

[TABLE]

In this illustration we set $\theta=5$, $p=0.3$, and $q=0.8$, and estimate these parameters on 100 data sets independently sampled from copula (23), with both our ABC procedure and the optimization procedure implemented in the fitCopula function of the copula package. The above procedure is repeated for samples of size $n=500$ and $n=10\,000$. In order to compare the results of our Bayesian procedure with the likelihood-based one, which provides only point estimates (maximum likelihood estimator, MLE), the posterior distribution of the model parameters is summarized into two point estimates: the posterior mean and the posterior median.

The results are plotted in Fig. 6, where each boxplot is made out of the 100 estimated values. In each panel, three quantities are plotted: MLE, which refers to the likelihood-based estimation, Post.median, which refers to the posterior median ABC estimation, and Post.mean, which refers to the posterior mean ABC estimation. It appears from the three plots that, as expected, the MLE is asymptotically unbiased. This can be seen by noticing the convergence to the true value as $n$ increases. On the other hand, both the mean and the median of the posterior distribution obtained with our ABC procedure show some bias. Unexpectedly, the results from MLE show a larger variance than ABC in some cases, particularly when $n=500$. This could be due to difficulties in the optimisation procedure of the fitCopula R function which, for small sample sizes, could have problems converging.

5 Conclusion

In this paper, we have studied the class of asymmetric copulas first introduced by Khoudraji [19], and developed in its general form by Liebscher [22]. Some new theoretical properties of these copulas were provived, including novel closed form expressions for its tail dependence coefficients, thus complementing the partial results of Liebscher [22] and Liebscher [23]. An iterative procedure is also introduced to flexibly sample from these copulas, which makes it easy to apply an Approximate Bayesian computation procedure to make inference on them.

Acknowledgments

The authors would like to thank two referees and the Associate Editor for their valuable suggestions, which have significantly improved the paper. The work of the third author benefited from the support of the Chair Stress Test, RISK Management and Financial Steering, led by the French École polytechnique and its Foundation and sponsored by BNP Paribas.

Appendix A Proofs of main results

Proof of Proposition 1

(i) Recall that $d=2$ and let $g^{(k)}\coloneqq g_{1}^{(k)}=g_{2}^{(k)}$ for all $k\in\{1,\dots,K\}$. For all $x>0$, we have

[TABLE]

since the product of the $g^{(k)}$ functions is the identity (a) and by definition of regular variation (b). It follows that $\sum_{k=1}^{K}\gamma^{(k)}=1$. Besides,

[TABLE]

by the regular variation property (a), by Lipschitz property for copulas and because of the continuity of $g^{(k)}$ at the origin (c). The result is thus proved.

(ii) Let $x\leq y$, the proof being similar when $x>y$. We have

[TABLE]

since, for any copula $C$ and any $(u,v)\in[0,1]^{2}$, $C(u,v)\leq\min(u,v)$. By assumption here, there exists $k_{0}\in\{1,\dots,K\}$ such that $\eta(t)\coloneqq g_{1}^{(k_{0})}(t)/g_{2}^{(k_{0})}(t)$, where $\eta(t)\to 0$ as $t\to 0$. Taking into account that $g_{1}^{(k_{0})}$ is increasing, this entails

[TABLE]

Recalling that $\prod_{k=1}^{K}g_{2}^{(k)}$ is the identity function, it follows that $\tilde{C}(\varepsilon x,\varepsilon y)\leq\eta(\varepsilon y)\varepsilon y$ and consequently $\Lambda_{L}(\tilde{C};x,y)=0$.

(iii) Differentiating the product of the $g_{j}^{(k)}$ functions, we obtain for $j\in\{1,2\}$:

[TABLE]

Now, since $g_{j}^{(\ell)}(1)=1$ for all $\ell\in\{1,\ldots,K\}$, and $j=1,2$, we obtain $\sum_{k=1}^{K}\text{d}_{j}^{(k)}=1$. Turning to the upper dependence function, let us write

[TABLE]

where $\stackrel{{\scriptstyle[}}{{\varepsilon}}\to 0]{}{\sim}$ denotes the asymptotic equivalence as $\varepsilon\to 0$. By assumption, all $g_{j}^{(k)}$ are differentiable at 1, with derivative denoted by $\text{d}_{j}^{(k)}$, and satisfy $g_{j}^{(k)}(1)=1$. Hence, the first order Taylor expansion is $g_{j}^{(k)}(1-\varepsilon)=1-\text{d}_{j}^{(k)}\varepsilon+o(\varepsilon)$. By the Lipschitz property of copulas,

[TABLE]

Plugging in to (24) yields

[TABLE]

Finally, in view of $\sum_{k=1}^{K}\text{d}_{j}^{(k)}=1$ for $j=1,2$ it follows that

[TABLE]

Taking the limit $\varepsilon\to 0$ yields the result:

[TABLE]

(iv) The case $\text{d}_{1}^{(k)}=\text{d}_{2}^{(k)}$ for all $k\in\{1,\dots,K\}$ is then a simple consequence of the homogeneity property of the upper tail dependence function. ∎

Proof of Corollary 1

This is a direct consequence of Proposition 1. ∎

Proof of Proposition 2

Liebscher [22, Proposition 2.2] proves the $\mathrm{TP}_{2}$ and LTD properties. The LTI property can be proven in the same way as LTD. For PQD, it suffices to remark that for any $u,v\in[0,1]$,

[TABLE]

while NQD works similarly with a reversed inequality in $(a)$.

Let us prove the SI part. According to Equation (6), $\mathrm{SI}$ is a property of $u\mapsto\tilde{C}(u,v)$ and $v\mapsto\tilde{C}(u,v)$ functions. Focusing without loss of generality on the former function, and omitting the $v$ variable for notational simplicity, $\tilde{C}$ can be written as

[TABLE]

with $g^{(1)}(u)\ldots g^{(K)}(u)=u$ for all $u\in[0,1]$. Differentiating this function twice yields $\tilde{C}^{\prime\prime}(u)=(T_{1}(u)+T_{2}(u))\tilde{C}(u)$ where

[TABLE]

By assumption, $C_{1},\ldots,C_{K}$ are SI, hence they are concave, thus $C_{1}^{\prime\prime}\leq 0$,…, $C_{K}^{\prime\prime}\leq 0$ and therefore $T_{1}\leq 0$. In view of Theorem 5.2.12 and Corollary 5.2.6 in Nelsen [28], $u\mapsto C_{\ell}^{\prime}(u)-{C_{\ell}(u)}/{u}$ is a negative function for all $\ell\in\{1,\ldots,K\}$. As a consequence, $\tau_{k}$ can be upper bounded as follows:

[TABLE]

where the equality is due to the fact that the product of all $g^{(\ell)}$ functions is the identity $\mathrm{Id}$, thus the derivative of product logarithm is $1/\mathrm{Id}$. Additionally, since $g^{(k)}$ is concave and $g^{(k)}(0)=0$, Theorem 5 by [5] implies that $-g^{(k)}$ is star-shaped, i.e., $-g^{(k)}/Id$ is increasing. This, in turn, proves that the right-hand side of (25) can be further upper bounded by zero, and therefore $T_{2}\leq 0$. As a conclusion, $\tilde{C}^{\prime\prime}<0$ and, in virtue of Equation (6), $\tilde{C}$ is also SI. The proof for SD follows similar lines. ∎

Proof of Proposition 3

By definition of the copula $\tilde{C}$,

[TABLE]

and the result is proved. ∎

Proof of Proposition 4

Let us remark that Example 2 shows that $C$ is max-stable implies $\tilde{C}^{(K)}=C$ for all $K\geq 1$ and for all sequence $(p^{(k)})_{k}\subset(0,1)$. Conversely, assume that $\tilde{C}^{(K)}=C$ for all $K\geq 1$ and for all sequence $(p^{(k)})_{k}\subset(0,1)$. From (1) with $K=2$, it follows that

[TABLE]

for all $\bm{u}\in[0,1]^{d}$. Introducing $\varphi:{\mathbb{R}}_{-}^{d}\to{\mathbb{R}}_{-}$ the continuous function defined by $\varphi=\ln\mathop{\circ}C\mathop{\circ}\exp$, we thus have $\varphi(\bm{v})=\varphi(p^{(1)}\bm{v})+\varphi((1-p^{(1)})\bm{v})$, for all $\bm{v}\in{\mathbb{R}}_{-}^{d}$ and $p^{(1)}\in(0,1)$. Lemma 1 in B entails that $\varphi$ is homogeneous of degree 1 or equivalently that $C$ is max-stable. ∎

Proof of Proposition 5

Let us first show that the copula $\tilde{C}^{(K)}$, $K\geq 1$ defined iteratively by (9), (10) is a Liebscher copula. The proof is done by induction on $K$. First, it is clear that

[TABLE]

Secondly, assume that

[TABLE]

and

[TABLE]

where the $\bigodot$ notation is defined in (11). From (10), it follows that

[TABLE]

by letting for all $j\in\{1,\dots,d\}$:

[TABLE]

As a first result, (28) proves (12) while (30) proves (13). Letting $\ell=K-k+1$ and in view of (27), equation (29) can be rewritten for all $2\leq\ell\leq K$ as

[TABLE]

which proves (14) for $\ell\in\{2,\dots,K\}$. The case $\ell=1$ is straightforward in view of (30).

Conversely, let us prove that any Liebscher copula can be constructed iteratively by (9), (10). Let $g_{j}^{(k)}$, $j\in\{1,\dots,d\}$, $k\in\{1,\dots,K\}$ be the set of functions associated with Liebscher copula (1). Our goal is to find a set of functions $f_{j}^{(k)}$, $j\in\{1,\dots,d\}$, $k\in\{1,\dots,K\}$ verifying the set of equations (9), (10). For all $j\in\{1,\dots,d\}$, define $f_{j}^{(K)}\coloneqq\mathrm{Id}/g_{j}^{(1)}$ and, for all $k\in\{K-1,\dots,1\}$,

[TABLE]

Let $j\in\{1,\dots,k\}$. The first part of the proof consists in establishing that $F_{j}^{(k+1)}$ is strictly increasing $[0,1]\to[0,1]$ and thus that its inverse $(F_{j}^{(k+1)})^{-1}$ is well-defined. To this end, remark that $F_{j}^{(K)}=\mathrm{Id}$ by definition and $F_{j}^{(k)}=f_{j}^{(k)}\mathop{\circ}F_{j}^{(k+1)}=F_{j}^{(k+1)}/g_{j}^{(K+1-k)}$, in view of (14). Iterating, it follows that

[TABLE]

in view of (2). It is then clear that $F_{j}^{(k+1)}$ is strictly increasing $[0,1]\to[0,1]$. The goal of the second part of the proof is to show that $f_{j}^{(k)}\in{\cal F}$ for all $k\in\{1,\dots,K\}$. It is clear that $\mathrm{Id}/f_{j}^{(k)}$ is strictly increasing for all $k\in\{1,\dots,K\}$ from (31). Besides, as already noticed, $f_{j}^{(k)}\mathop{\circ}F_{j}^{(k+1)}=F_{j}^{(k)}$, and thus $f_{j}^{(k)}$ is strictly increasing as the composition of strictly increasing functions. This concludes the proof that $f_{j}^{(k)}\in{\cal F}$. The third part of the proof consists in showing that $f_{j}^{(k)}$, $k\in\{1,\dots,K\}$ is solution of the set of equations (13), (14). It is readily seen that (13) holds and

[TABLE]

which proves (14). ∎

Proof of Proposition 6

The expression (19) of the copula $\tilde{C}_{\mathrm{CL}}$ is a simple consequence of the definition of the partition $\mathcal{A}_{0},\ldots,\mathcal{A}_{K}$. We only derive the singular component expression in the case of a product of two terms, $K=2$; the general case follows similar lines. Denote $p_{1}$ and $q_{1}$ by $p$ and $q$, then $\tilde{C}_{\mathrm{CL}}(u,v)$ can be respectively written on $\mathcal{A}_{0}$, $\mathcal{A}_{1}$, and $\mathcal{A}_{2}$ by $u$, $u^{1-p}v^{q}$ and $v$. The cross derivative $\frac{\partial^{2}\tilde{C}_{\mathrm{CL}}}{\partial u\partial v}(x,y)$ then vanishes on $\mathcal{A}_{0}$ and $\mathcal{A}_{2}$, and is equal to $(1-p)qx^{-p}y^{q-1}$ on $\mathcal{A}_{1}$. Using the formula [28, Eq. (2.4.1)],

[TABLE]

and dividing the double integral above into the three sets of $[0,u]\times[0,v]$ intersected with $\mathcal{A}_{0}$, $\mathcal{A}_{1}$, and $\mathcal{A}_{2}$ yields

[TABLE]

and therefore routine calculations yield the following expressions:

[TABLE]

Taking the difference $\tilde{S}_{\mathrm{CL}}=\tilde{C}_{\mathrm{CL}}-\tilde{A}_{\mathrm{CL}}$ yields a simple expression for the singular component

[TABLE]

Given that $r_{1}\leq r_{2}$ yields $p\leq q$ in the case $K=2$, thus leading to expression (20) when $K=2$. The general case $K>2$ is derived similarly. ∎

Proof of Proposition 7

(i) Taking account of the identity $\min(x,y)+\max(x,y)=x+y$ yields

[TABLE]

and consequently $\beta(\tilde{C}_{\mathrm{CL}})=2^{\sum_{k=1}^{K}\min(p_{k},q_{k})}-1$.

(ii) We use the following expression of Kendall’s $\tau$, which is convenient since the copula $\tilde{C}_{\mathrm{CL}}$ may have non null singular component ([28], Eq. 5.1.12):

[TABLE]

Let $k\in\{0,\dots,K\}$. Then, for all $(u,v)\in{\cal A}_{k}$, one has

[TABLE]

Besides, remarking that $\bar{q}_{0}=0$ and $\bar{p}_{K}=1$ shows that both terms $k=0$ and $k=K$ do not contribute to the sum in $\tau(\tilde{C}_{\mathrm{CL}})$. The result

[TABLE]

then follows.

(iii) Recall that

[TABLE]

and the result is then established. ∎

Appendix B Auxiliary results

Lemma 1**.**

Let $\varphi:{\mathbb{R}}_{-}^{d}\to{\mathbb{R}}_{-}$ be a continuous function such that

[TABLE]

for all $a\in(0,1)$ and $x\in{\mathbb{R}}_{-}^{d}$. Then, necessarily, $\varphi$ is homogeneous of degree 1.

Proof

Firstly, let us prove by induction the property $(P_{n})$: $\varphi(x/n)=\varphi(x)/n$ for all $n\in{\mathbb{N}}\setminus\{0\}$ and $x\in{\mathbb{R}}_{-}^{d}$. $(P_{1})$ is straightfowardly true. Assume $(P_{n})$ holds. Then,

[TABLE]

and (33) entails

[TABLE]

which proves $(P_{n+1})$.

Second, for all $m\in{\mathbb{N}}\setminus\{0\}$, $(P_{m})$ shows that $\varphi(x)=m\varphi(x/m)$ and thus, letting $y=x/m$, $\varphi(my)=m\varphi(y)$ for all $y\in{\mathbb{R}}_{-}^{d}$. This property can be extended to $m=0$ since letting $a\to 0$ in (33) yields $\varphi(0)=0$.

Third, let $q\in{\mathbb{Q}}_{+}$. There exists $(m,n)\in{\mathbb{N}}\times{\mathbb{N}}\setminus\{0\}$ such that $q=m/n$. From the first two points, $\varphi(qx)=\varphi(mx/n)=m\varphi(x/n)=m\varphi(x)/n=q\varphi(x)$.

Finally, the continuity of $\varphi$ and the density of ${\mathbb{Q}}_{+}$ in ${\mathbb{R}}_{+}$ imply that $\varphi(tx)=t\varphi(x)$ for all $t\in{\mathbb{R}}_{+}$ and $x\in{\mathbb{R}}_{-}^{d}$. The result is thus proved. ∎

Lemma 2**.**

For all $k\geq 1$ let $C_{k}$ be a $d$-variate copula and $f^{(k)}_{j}\in{\cal F}$ for all $j\in\{1,\dots,d\}$, with the assumption $f^{(1)}_{j}(t)=1$ for all $t\in[0,1]$. The sequence $(\tilde{C}^{(k)})_{k\geq 1}$ defined iteratively by (9) and (10) is a sequence of $d$-variate copulas.

Proof

The proof is done by induction on $k$. Let ${\cal C}$ be the set of all $d$-variate copulas. First, it is clear that $\tilde{C}^{(1)}\in{\cal C}$ from (9). Second, let us assume that $\tilde{C}^{(k-1)}\in{\cal C}$ and prove that $\tilde{C}^{(k)}\in{\cal C}$. for all $k\geq 2$. Let $(Y_{1},\dots,Y_{d})$ and $(Z_{1},\dots,Z_{d})$ be two independent random vectors in $[0,1]^{d}$ drawn respectively from the cdf $\tilde{C}^{(k-1)}(f_{1}^{(k)},\dots,f_{d}^{(k)})$ and ${C}_{K}(\mathrm{Id}/f_{1}^{(k)},\dots,\mathrm{Id}/f_{d}^{(k)})$. For all $j\in\{1,\dots,d\}$ define the random variable $X_{j}=\max(Y_{j},Z_{j})$. For all $\bm{u}\in[0,1]^{d}$, the cdf of $(X_{1},\dots,X_{d})$ is given by

[TABLE]

from (10). This proves that $\tilde{C}^{(k)}$ is a cdf. The margins are uniform by construction. ∎