Exchangeable FGM copulas

Abstract

Copulas provide a powerful and flexible tool for modeling the dependence structure of random vectors, and they have many applications in finance, insurance, engineering, hydrology, and other fields. One well-known class of copulas in two dimensions is the Farlie–Gumbel–Morgenstern (FGM) copula, since its simple analytic shape enables closed-form solutions to many problems in applied probability. However, the classical definition of the high-dimensional FGM copula does not enable a straightforward understanding of the effect of the copula parameters on the dependence, nor a geometric understanding of their admissible range. We circumvent this issue by analyzing the FGM copula from a probabilistic approach based on multivariate Bernoulli distributions. This paper examines high-dimensional exchangeable FGM copulas, a subclass of FGM copulas. We show that the dependence parameters of exchangeable FGM copulas can be expressed as a convex hull of a finite number of extreme points. We also leverage the probabilistic interpretation to develop efficient sampling and estimating procedures and provide a simulation study. Throughout, we discover geometric interpretations of the copula parameters that assist one in decoding the dependence of high-dimensional exchangeable FGM copulas.

1. Introduction

Copulas are a powerful tool for modeling dependence between the components of a random vector. A well-known family of copulas is that of the Farlie–Gumbel–Morgenstern (FGM) copulas, first studied by [17, 18, 26, 42].

FGM copulas are attractive since their simple shape enables exact calculus. Being quadratic in each marginal, FGM copulas allow one to develop closed-form expressions for many quantities of interest. For a given set of dependence parameters, many basic properties of FGM copulas are known; see, for instance, [9, 31, 38], [35, Chapter 5], [34, Section 44.10], [16], or [47]. FGM copulas have been applied in many disciplines, including, for instance, finance [38], actuarial science [5], bioinformatics [32], and hydrology [22]. However, it is not yet clear how to interpret the FGM copula parameters in higher dimensions—in particular, it is not clear how the copula parameters affect the dependence structure and how one may compare dependence constructions in terms of dependence orders.

In [7], the authors establish a one-to-one correspondence between the family of d-variate FGM copulas and the family of d-variate symmetric multivariate Bernoulli distributions. By symmetric Bernoulli distributions, we mean a discrete random variable (RV) taking the value 1 with probability 1/2 and 0 with probability 1/2. One advantage of this representation is that one may construct subfamilies of FGM copulas by selecting subfamilies of multivariate symmetric Bernoulli distributions. Then the subfamily of FGM copulas will share the dependence properties of the symmetric Bernoulli distributions, thus simplifying the parameter space and related operations like sampling and estimation. Another advantage of the stochastic representation of FGM copulas is that multivariate Bernoulli distributions are simpler to understand. While FGM copulas only induce weak dependence, they are the simplest of the Bernstein copulas. A d-dimensional Bernstein copula, introduced by [54], is dense on the hypercube $[0, 1]^d$ , but many dependence parameters are required to specify the copula. Indeed, because of their flexibility, Bernstein copulas are sometimes used as alternatives to the empirical copula, as in [56]. Understanding the stochastic nature of FGM copulas is essential preliminary work before investigating the properties, geometries, and stochastic nature of Bernstein copulas. In light of the new results presented in this paper, we will return to a stochastic representation of an exchangeable Bernstein copula in the conclusion.

The present paper investigates exchangeable FGM (eFGM) copulas. The eFGM copulas are subfamilies of FGM copulas that we construct with exchangeable symmetric multivariate Bernoulli random vectors. It follows that eFGM copulas are the simplest of the exchangeable Bernstein copulas and that understanding the geometry and properties of eFGM copulas lays the groundwork for extending these to exchangeable Bernstein copulas. In turn, one will be better positioned to investigate the general class of Bernstein copulas. The exchangeability assumption is reasonable and useful in some contexts; consider, for instance, the study of litter-mates in laboratory experiments [36], finance [49], reliability theory [45], actuarial science [33], or credit default risk [20, 41]. Exchangeability also plays an important role in Bayesian statistics [55].

Another advantage of studying the class of eFGM copulas is that an FGM copula corresponding to the lower bound under the supermodular order is a special case of eFGM copulas. We study this lower bound in detail in Section 6. We also introduce subfamilies of eFGM copulas that display a specific shape of dependence structure, and one may compare copulas under the supermodular order within the subfamilies. Ordering of random vectors with respect to the supermodular order is important for practical applications. For instance, in applied probability, finance, and actuarial science, one may be interested in the distribution of the sum of the components of a random vector. If one may order two copulas under the supermodular order, then one may order the two aggregate distributions under the stop-loss order, which implies inequalities of certain useful risk measures. See, for instance, Section 8.3 of [44] or Section 6.3 of [13] for details on the supermodular order, the stop-loss order, and aggregate distributions.

The remainder of this paper is organized as follows. In Section 2, we introduce the subclass of eFGM copulas. Section 3 presents construction methods for symmetric exchangeable Bernoulli RVs and their relationship to eFGM copulas. In Section 4, we show that the parameters of all eFGM copulas can be expressed as a convex hull of eFGM copula dependence parameters. We also provide a method for analytically obtaining extreme points corresponding to the copula parameters. One can view an eFGM copula as a finite mixture model. By studying the extreme points of the convex hull of the dependence parameters of d-variate eFGM copulas, we gain a geometric understanding of the class of eFGM copulas. We characterize, in Section 5, the class of d-variate eFGM copulas that can be represented as the first elements of an infinite sequence of RVs. We deal with dependence ordering in Section 6, providing methods to compare d-variate eFGM copulas under the supermodular order. In Section 7, we discuss sampling and estimation for high-dimensional eFGM copulas. In particular, we leverage the stochastic representation of eFGM copulas to propose an efficient stochastic sampling method, and we leverage the finite mixture representation to propose an estimation algorithm. In Section 8, we offer some conclusions and discussions for future research. Certain proofs are deferred to the appendix.

2. Definition

In this section, we introduce the subfamily of copulas studied in the paper. First, recall that copulas are multivariate cumulative distribution functions (CDFs) of RVs with uniform marginals.

Definition 1. A (d-variate) copula is a function $C\,:\, [0, 1]^d \to [0, 1]$ satisfying the following conditions:

1. 1. $C(u_1, \dots, u_d) = 0$ if any $u_j = 0$ , $j \in \{1, \dots, d\}$ .

2. 2. $C(u_1, \dots, u_d) = u_j$ if $u_k = 1$ for all $k \in \{1, \dots, d\}$ and $k \neq j$ .

3. 3. C is d-increasing on $[0, 1]^d$ ; that is,

   $$\sum_{i_1 = 1}^2\dots \sum_{i_d = 1}^2 ({-}1)^{i_1 + \dots + i_d}C(u_{1i_1}, \dots, u_{di_d}) \geq 0$$
   for all $0 \leq u_{j1} \leq u_{j2} \leq 1$ and $j \in \{1, \dots, d\}.$

An important family of copulas is that of the FGM copulas, first studied by [17, 18, 26, 42]. One may refer to [9, 31], [35, Chapter 5], [34, Section 44.10], or [22] for properties of this family of copulas. A d-variate FGM copula is defined as

(1) \begin{align} C\left( u_1, \dots, u_d\right) = & \left(\prod_{j=1}^d u_j\right) \left( 1+\sum_{k=2}^{d}\sum_{1\leq j_{1}<\dots <j_{k}\leq d}\theta_{j_{1}\dots j_{k}}\overline{u}_{j_{1}}\overline{u}_{j_{2}}\dots \overline{u}_{j_{k}}\right)\!,\nonumber\\ & \qquad\qquad\qquad\qquad\qquad\qquad\qquad (u_1, \dots, u_d)\in [0,1]^{d},\end{align}

where $\overline{u}_{j}=1-{u}_{j}$ , $j \in \{1,\dots,d\}$ . The set of admissible parameters for FGM copulas, derived in [9], is given by

(2) \begin{equation} \left\{(\theta_{12}, \dots, \theta_{1\dots d}) \in \mathbb{R}^{2^d - d - 1}\, :\, 1+\sum_{k=2}^{d}\sum_{1\leq j_{1}<\cdots <j_{k}\leq d}\theta_{j_{1}\ldots j_{k}}\varepsilon _{j_{1}}\varepsilon _{j_{2}}\ldots \varepsilon_{j_{k}}\geq 0\right\}\!,\end{equation}

for all $\{\varepsilon_{j_{1}},\varepsilon _{j_{2}},\dots, \varepsilon_{j_{k}}\} \in \{-1,1\}^d$ . We call the $\binom{d}{k}$ parameters $\theta_{j_{1}\dots j_{k}}$ , for $1 \leq j_1 < \dots < j_k \leq d$ , the k-dependence parameters, $k \in \{2, \dots, d\}$ . A d-variate FGM copula has $2^d - d - 1$ parameters; the large number of parameters becomes impractical for high-dimensional applications of FGM copulas. However, one may rely on a new stochastic representation of FGM copulas, as introduced in the following theorem from [7].

Theorem 1. The copula in (1) has an equivalent representation

(3) \begin{equation} C(u_1, \dots, u_d) = {\mathbb{E}}\left[ \prod_{m = 1}^d u_m \left(1 + ({-}1)^{I_m}\overline{u}_m\right)\right], \end{equation}

for $\left(u_1, \dots, u_d\right) \in [0, 1]^d$ , where $\boldsymbol{I} = (I_1, \dots, I_d)$ is a symmetric multivariate Bernoulli random vector with

(4) \begin{equation} \theta_{j_1 \dots j_k} = ({-}2)^k {\mathbb{E}}\left[\prod_{j = 1}^k\left(I_{j} - \frac{1}{2}\right)\right], \quad k \in \{2, \dots, d\}. \end{equation}

In particular, the probability mass function (PMF) of the underlying random vector associated with an FGM copula is given by

(5) \begin{equation} f_{\boldsymbol{I}}(\boldsymbol{i}) = \frac{1}{2^d}\left(1 + \sum_{k = 2}^d\sum_{1\leq j_{1}<\cdots <j_{k}\leq d}({-}1)^{i_{j_1}+\dots + i_{j_k}}\theta_{j_1\dots j_k}\right)\!, \end{equation}

for $\boldsymbol{i} \in \{0, 1\}^d$ , while the copula parameters associated with the PMF of a symmetric multivariate Bernoulli random vector are

\begin{equation*} \theta_{j_1 \dots j_k} = \sum_{\left(i_{j_1}, \dots, i_{j_k}\right)\in \{0, 1\}^k}({-}1)^{i_{j_1} + \dots + i_{j_k}} f_{I_{j_1}, \dots, I_{j_k}}\left(i_{j_1}, \dots, i_{j_k}\right)\!, \end{equation*}

for $1 \leq j_1 < \dots < j_k \leq d$ and $k \in \{2, \dots, d\}$ .

The current paper studies the subfamily of eFGM copulas that have the shape

(6) \begin{equation} C_d(u_1, \dots, u_d) = \left(\prod_{j = 1}^d u_j\right)\left(1 + \sum_{k = 2}^d \sum_{1\leq j_1< \dots < j_k\leq d} \theta_k \overline{u}_{j_1} \dots \overline{u}_{j_k}\right)\!, \quad (u_1, \dots, u_d) \in [0, 1]^d,\end{equation}

for $d \geq 2$ . For $k \in \{2, \dots, d\}$ , this class of FGM copulas sets each of the $\binom{d}{k}$ parameters $\theta_{j_{1}\dots j_{k}} = \theta_k$ for all $1 \leq j_1 < \dots < j_k \leq d$ ; that is, all k-dependence parameters are equal. By symmetry of the bivariate FGM copula, it is obvious that with $d = 2$ , each admissible parameter $\theta_2 \in [-1, 1]$ corresponds to an exchangeable bivariate FGM copula parameter; that is, the entire class of bivariate FGM copulas are also eFGM copulas.

A d-variate eFGM copula is specified by a vector of $d - 1$ parameters $(\theta_2, \dots, \theta_d) \in \mathcal{T}_d$ (as opposed to $2^d - d - 1$ for the complete class of d-variate FGM copulas) where, for all $\{\varepsilon_{1},\dots, \varepsilon_{d}\} \in \{-1,1\}^d$ , the set of admissible parameters is

(7) \begin{equation} \mathcal{T}_d=\left\{(\theta_2, \dots, \theta_d) \in \mathbb{R}^{d-1}\, :\, 1+\sum_{k=2}^{d}\sum_{1\leq j_{1}<\cdots <j_{k}\leq d}\theta_{k}\varepsilon _{j_{1}}\dots \varepsilon _{j_{k}}\geq 0\right\}.\end{equation}

As the dimension d increases, satisfying the $2^{d-1}$ constraints for the parameters $(\theta_2, \dots, \theta_d)$ in (7) becomes tedious (it is a computation in exponential time). A preferable approach to studying eFGM copulas is by constructing a stochastic representation, along the same lines as in Theorem 1.

Corollary 1. Let C be a copula as in (6). Then C also admits a stochastic representation as in (3) if and only if $\boldsymbol{I}$ is an exchangeable symmetric Bernoulli random vector.

Proof. We first show that if $\boldsymbol{I}$ is an exchangeable symmetric Bernoulli random vector, then C is an eFGM copula. For all $1 \leq j_1 < \dots < j_k \leq d$ , $1 \leq j_1^{\prime}<\dots < j_k^{\prime} \leq d$ , and $k \in \{2, \dots, d\}$ , we have from (4) that

$$\theta_{j_1 \dots j_k} = ({-}2)^k {\mathbb{E}}\left[\prod_{l = 1}^k\left(I_{j_l} - \frac{1}{2}\right)\right] \quad \text{and} \quad \theta_{j_1^{\prime} \dots j_k^{\prime}} = ({-}2)^k {\mathbb{E}}\left[\prod_{l = 1}^k\left(I_{j_l}^{\prime} - \frac{1}{2}\right)\right].$$

By exchangeability of $\boldsymbol{I}$ , we have that

$${\mathbb{E}}\left[\prod_{l = 1}^k\left(I_{j_l} - \frac{1}{2}\right)\right] = {\mathbb{E}}\left[\prod_{l = 1}^k\left(I_{j_l}^{\prime} - \frac{1}{2}\right)\right]$$

for all $1 \leq j_1 < \dots < j_k \leq d$ , $1 \leq j_1^{\prime}<\dots < j_k^{\prime} \leq d$ , and $k \in \{2, \dots, d\}$ ; hence $\theta_{j_1 \dots j_k} = \theta_{j_1^{\prime}\dots j_k^{\prime}}$ .

Next, we must show that if C is an eFGM copula, then the underlying symmetric multivariate Bernoulli random vector $\boldsymbol{I}$ is exchangeable. Replacing $\theta_{j_1\dots j_k}$ by $\theta_k$ in (5) for all $1\leq j_1 < \dots < j_k \leq d$ and $k \in \{2, \dots, d\}$ in (5), we see that $f_{\boldsymbol{I}}(i_1, \dots, i_d) = f_{\boldsymbol{I}}(i_1^{\prime}, \dots, i_d^{\prime})$ whenever $i_1 + \dots + i_d = i_1^{\prime} + \dots + i_d^{\prime}$ ; it follows by the associative property of addition that $\boldsymbol{I}$ is exchangeable.

We offer an interpretation of the dependence structure within eFGM copulas in the following. Let $V_1$ and $V_2$ be a pair of independent standard uniform RVs. Define $U_{[1]} \overset{\mathcal{D}}{=} \min(V_1, V_2)$ and $U_{[2]} \overset{\mathcal{D}}{=} \max(V_1, V_2)$ , where $\overset{\mathcal{D}}{=}$ means equality in distribution. We have that $U_{[1]}$ is beta distributed with CDF $F_{U_{[1]}}(u) = u(2-u)$ , and $U_{[2]}$ is beta distributed with CDF $F_{U_{[2]}}(u) = u^2$ , for $0\leq u\leq 1$ . Let $\boldsymbol{U}_{[j]}$ be a vector of independent RVs where each component has CDF $F_{U_{[j]}}$ , for $j \in \{1, 2\}$ . Let $\boldsymbol{U}$ be a random vector whose joint CDF corresponds to an eFGM copula. Then we have from Theorem 1 that there exists a random vector $\boldsymbol{I}$ such that

(8) \begin{equation} \boldsymbol{U} = (\boldsymbol{1} - \boldsymbol{I}) \boldsymbol{U}_{[1]} + \boldsymbol{I}\boldsymbol{U}_{[2]},\end{equation}

where $\boldsymbol{1}$ is a vector of ones. Within the context of this paper, we require that $\boldsymbol{I}$ is a vector of exchangeable RVs.

3. Construction methods and examples

We first present a few methods of constructing exchangeable symmetric Bernoulli RVs. Alternating between different construction methods, along with the stochastic representation of eFGM copulas in Corollary 1, will enable us to study the properties of eFGM copulas.

3.1. Construction based on the sum of Bernoulli RVs

We first define the (univariate) RV $N_d = \sum_{j = 1}^{d} I_j$ , with support $\{0, 1, \dots, d\}$ , representing the sum of d exchangeable Bernoulli RVs. The relationship between the PMF of $(I_1, \dots, I_d)$ and $N_d$ is

\begin{align*} {\mathbb{P}}(N_d = k) &= \sum_{\substack{\{i_1, \dots, i_d\} \in \{0, 1\}^d \\ i_\bullet = k}} {\mathbb{P}}\left(I_1 = i_1, \dots, I_d = i_d\right)\\ &= \binom{d}{k} {\mathbb{P}}(I_1 = 1, \dots, I_k = 1, I_{k+1} = 0, \dots, I_d = 0),\end{align*}

where $i_\bullet = \sum_{j = 1}^d i_j$ ; the second equality follows by exchangeability of $(I_1, \dots, I_d)$ .

Let $\mathcal{N}_d$ represent the class of PMFs for univariate RVs with support $\{0, \dots, d\}$ with mean $d/2$ . In [20, Section 3.2], the authors provide a one-to-one correspondence between the class of PMFs for d-variate exchangeable Bernoulli random vectors and $\mathcal{N}_d$ . This construction is useful since it identifies the joint PMF of $(I_1, \dots, I_d)$ only through the PMF of $N_d$ .

3.2. Construction based on a vector of probabilities

One can specify the multivariate distribution of $(I_1, \dots, I_d)$ by the vector of probabilities $(\zeta_0, \dots, \zeta_d)$ , where $\zeta_0 = 1$ and $\zeta_k = {\mathbb{P}}(I_1 = 1, \dots, I_k = 1)$ , $k \in \{1,\dots,d\}$ and $d \in \{1, 2, \dots\}$ . For eFGM copulas, we require $\zeta_1 = 1/2$ . We recall [37, Theorem 1], which provides a sufficient condition for the values of $\zeta_k$ , $k \in \{1, \dots, d\}$ .

Theorem 2. Let $d\in \{2, 3, \dots\}$ be fixed, and let $\psi(t)$ be a completely monotone function for $t \geq 0$ . If $\zeta_k = \psi(k)$ , then ${\mathbb{P}}(N_d = k) \geq 0$ , for $k \in \{0, 1, \dots, d\}$ .

The relationship between the values $(\zeta_0, \dots, \zeta_d)$ , characterizing the multivariate distribution of a vector of d exchangeable RVs $(I_1, \dots, I_d)$ , and the values of the components of the vector of dependence parameters $(\theta_2, \dots, \theta_d)$ of the eFGM copula is established in the next result.

Corollary 2. Let $d \in \{2, 3, \dots \}$ be fixed. For a given vector $(\zeta_0, \dots, \zeta_d)$ satisfying the conditions of Theorem 2 with $\zeta_0 = 1$ , $\zeta_1 = 1/2$ , we have

(9) \begin{align} \theta_{k} &= ({-}2)^{k} \sum_{l = 0}^{k}\binom{k}{l} \zeta_{l}\left({-}\frac{1}{2}\right)^{k-l} = \sum_{l = 0}^{k}\binom{k}{l} \zeta_{l}\left({-}2\right)^{l}, \quad k\in \{2,\ldots,d\}. \end{align}

Proof. Expanding the product in (4) yields

\begin{align*} \theta_k &= ({-}2)^k {\mathbb{E}}\left[\left({-}\frac{1}{2}\right)^k + \sum_{l = 1}^k I_l\left({-}\frac{1}{2}\right)^{k-1} + \sum_{l = 2}^k \sum_{1 \leq j_1 < \dots < j_l \leq k} I_{j_1} \dots I_{j_l} \left({-}\frac{1}{2}\right)^{k - l} \right]. \end{align*}

Since $\left(I_1, \dots, I_d\right)$ are exchangeable random vectors, one has ${\mathbb{E}}\left[I_{j_1}\dots I_{j_k}\right] = \zeta_k$ for all k-dimensional vectors $(j_1,\dots,j_k)$ such that $1 \leq j_1 < \dots < j_k \leq d$ and $k \in \{2,\dots,d\}$ .

Since (9) does not depend on d, the first k-dependence parameters from Corollary 2 are

(10) \begin{equation} \theta_{k} = \begin{cases} 4\zeta_2 - 1,& k = 2,\ d \geq 2,\\[3pt] -8\zeta_3 + 12 \zeta_2 - 2, & k = 3,\ d\geq 3,\\[3pt] 16\zeta_4 - 32 \zeta_3 + 24\zeta_2 - 3, & k = 4,\ d \geq 4,\\[3pt] -32\zeta_5 + 80\zeta_4 - 80\zeta_3 + 40\zeta_2 - 4, & k = 5,\ d \geq 5. \end{cases}\end{equation}

Example 1. (Model 3 of [37].) Madsen [37] considered the model $\zeta_k = \beta + (1 - \beta)\alpha^k$ , for $(\alpha, \beta) \in [0, 1]^2$ and $k \in \{0, 1, \dots\}$ . With the constraint that $\zeta_1 = 1/2$ , we have $\alpha = (1/2-\beta)/(1 - \beta)$ for $\beta \in [0, 0.5]$ , and we have one free parameter $\beta$ , meaning that the parameters $(\theta_2, \dots, \theta_d)$ are entirely determined by $\beta$ . Inserting these probabilities in (9) yields $\theta_k = \beta({-}1)^k + (1 - \beta)(1 - (1-2\beta)/(1 - \beta))^k$ , for $k \in \{2, 3, \dots\}$ . The case $\beta = 0$ yields the independence copula, and $\beta = 0.5$ yields the extreme positive dependence copula $C^{EPD}$ that we will describe in Theorem 8.

A more convenient way of specifying the values of $\zeta_k$ for $k \in \{0, 1, \dots\}$ uses Laplace–Stieltjes transforms. Let us first recall Bernstein’s theorem, originally from [6]; see also Theorem 1a in [19, Section XIII.4].

Theorem 3. If $\psi(0) = 1$ and $\psi$ is completely monotone, then $\psi$ is the Laplace–Stieltjes transform of a strictly positive RV Y; that is, $\psi(t) = \mathcal{L}_Y(t) = {\mathbb{E}}\left[e^{-Yt}\right]$ .

Corollary 3. Setting $\zeta_k$ , $k = 0, 1, \dots, d$ , to $\zeta_k = \mathcal{L}_Y(rk)$ , for $r > 0$ , will generate probability values which satisfy the conditions of Theorem 2. For a symmetric multivariate exchangeable Bernoulli random vector, we have $\zeta_1 = \mathcal{L}_Y(r) = 1/2$ , implying $\zeta_k = \mathcal{L}_Y\left(k \times \mathcal{L}_Y^{-1}\left(1/2\right)\right)$ , for $k = 1, \dots, d$ .

Remark 1. As shown in [24], the constructions based on $N_d$ from Subsection 3.1 and on the vector of probabilities $(\zeta_0, \dots, \zeta_d)$ of the current subsection are equivalent and related through the relationship

$${\mathbb{P}}(N_d = k) = \binom{d}{k} \sum_{l = 0}^{d-k} ({-}1)^l \binom{d - k}{l}\zeta_{k+l}, \quad k \in \{0, 1, \dots, d\}.$$

3.3. Construction with mixtures

One can also construct exchangeable Bernoulli distributions using mixtures, that is,

(11) \begin{equation} \mathbb{P}(I_1 = i_1, \dots, I_d = i_d) = \int_{0}^{1} \mathbb{P}(I_1 = i_1 \vert \Lambda = \lambda) \times \dots \times \mathbb{P}(I_d = i_d \vert \Lambda = \lambda) \mathrm{d}F_{\Lambda}(\lambda),\end{equation}

where $\Lambda$ is a mixing RV defined on [0,1]. According to (11), conditional on the mixing RV $\Lambda$ , $\left(I_1, \dots, I_d\right)$ are conditionally independent. One must select a distribution for $\Lambda$ such that ${\mathbb{E}}\left(\Lambda\right) = 1/2$ . From (11), it follows that

\begin{equation*} \zeta_k = {\mathbb{P}}(I_1 = 1, \dots, I_k = 1) = \int_0^1 \lambda^k {\, \mathrm{d}} F_{\Lambda}(\lambda) = {\mathbb{E}}\left[\Lambda^k\right],\end{equation*}

for $k \in \{0, \dots, d\}$ , and

(12) \begin{align} f_{I_1, \dots, I_d}(i_1, \dots, i_d) = \int_0^1 \lambda^{i_\bullet}(1 - \lambda)^{d - i_\bullet} {\, \mathrm{d}} F_{\Lambda}(\lambda) = {\mathbb{E}}\left[\Lambda^{i_\bullet}\left(1-\Lambda\right)^{d-i_\bullet}\right],\end{align}

for $(i_1, \dots, i_d) \in \{0, 1\}^d$ . Furthermore, for $k \in \{2,\ldots,d\}$ , combining (4) and (11), the parameters of the copula are defined in terms of the central mixed moments of $\Lambda$ as follows:

(13) \begin{equation} \theta_k = ({-}2)^k {\mathbb{E}}_{\Lambda}\left[{\mathbb{E}}\left\{\left.\prod_{j = 1}^k \left(I_j - \frac{1}{2}\right)\right\vert \Lambda\right\} \right] = ({-}2)^k {\mathbb{E}}\left[\left(\Lambda - \frac{1}{2}\right)^k\right].\end{equation}

Remark 2. Let $(I_1, \dots, I_d, I_{d + 1}, \dots\! )$ be an infinite sequence of exchangeable symmetric Bernoulli RVs. Let $(I_1, \dots, I_d)$ be the first d RVs from that sequence. The famous result from de Finetti [12] states that there exists an RV $\Lambda$ such that (11) holds. On the other hand, if we have (11), then $(I_1, \dots, I_d)$ can be extended to higher dimensions. We will study the extendability of eFGM copulas in Section 5.

Remark 3. Using the mixture construction and (13), one can interpret the copula parameters from the mixing RV. We have $\theta_{2}\propto Var(\Lambda)$ , which implies that the variance of the mixing RV induces the 2-dependence parameters. Then, since $\theta_{3} \propto - {\mathbb{E}}[\{\Lambda - {\mathbb{E}}(\Lambda)\}^3]$ , we interpret $\theta_{3}$ as proportional to the negative of the skewness. If the density function of $\Lambda$ is symmetric about the mean (skewness of 0), then $\theta_{k} = 0$ when k is odd. When k is even, we have ${\mathbb{E}}[\{\Lambda - {\mathbb{E}}(\Lambda)\}^k] \geq 0$ , implying $\theta_{k} \geq 0$ , so the mixture construction induces positive dependence.

A family of distributions for $\Lambda$ will generate a specific family of eFGM copulas. The following example presents the beta-eFGM family of copulas with $\Lambda$ set to follow a beta distribution.

Example 2. Let $\Lambda \sim Beta(\alpha, \alpha)$ for $\alpha > 0$ with probability density function

$$f_{\Lambda}(\lambda) = \frac{[\lambda(1-\lambda)]^{\alpha - 1}}{B(\alpha, \alpha)}, \quad 0\leq \lambda \leq 1,$$

in which case $\zeta_1 = E\left(\Lambda\right) = 1/2$ . The only parameter within this example is $\alpha$ ; hence it acts as a genuine dependence parameter, meaning that the vector $(\theta_2, \dots, \theta_d)$ is entirely determined by $\alpha$ . The PMF in (11) becomes

(14) \begin{equation} f_{I_1, \dots, I_d}(i_1, \dots, i_d) = \frac{\Gamma(2\alpha)}{\Gamma(\alpha)^2} \frac{\Gamma(2\alpha + d)}{\Gamma\left(\alpha + i_\bullet\right)\Gamma\left(\alpha + d - i_\bullet\right)} = \frac{B\left(\alpha + i_\bullet, \alpha + d - i_\bullet\right)}{B(\alpha, \alpha)}, \end{equation}

where B(a, b) is the beta function $\Gamma(a)\Gamma(b)/\Gamma(a+b)$ ; see Chapter 7 of [28] for details. Note that the distribution of $N_d$ , when $(I_1, \dots, I_d)$ has PMF (14), is called the beta-binomial distribution. From (13), the dependence parameters are $\theta_k = B(\alpha+1/2, (k+1)/2)/B(\alpha+(k+1)/2, 1/2)$ for $k = 2, 4, 6, \dots$ and $\theta_k = 0$ for $k = 3, 5, 7, \dots$ . We provide a proof in Appendix A. The beta-eFGM copula is

(15) \begin{equation} C\left(u_1, \dots, u_d\right) = \prod_{j = 1}^{d} u_j \left(1 + \sum_{l = 1}^{\left\lfloor \frac{d}{2} \right\rfloor}\sum_{1\leq j_{1}<\cdots <j_{2l}\leq d} \frac{B(\alpha+1/2, l + 1/2)}{B(\alpha+l+1/2, 1/2)}\overline{u}_{j_1}\cdots \overline{u}_{j_{2l}}\right)\!, \end{equation}

for $\left(u_1, \dots, u_d\right) \in [0,1]^{d}$ , where $\lfloor y \rfloor$ is the floor function returning the greatest integer smaller than or equal to y. We have the following representations for the dependence parameters when k is even:

(16) \begin{align} \theta_{k} &= \begin{cases} \frac{1}{2\alpha + 1},& k = 2,\ d \geq 2,\\[5pt] \frac{3}{(2\alpha + 1)(2\alpha + 3)}, & k = 4,\ d \geq 4,\\[5pt] \frac{15}{(2\alpha + 1)(2\alpha + 3)(2\alpha + 5)}, & k = 6,\ d \geq 6,\\ \end{cases} \nonumber \\[6pt] \theta_{k} &= \prod_{l = 1}^{k/2} \frac{2l - 1}{2\alpha + 2l - 1}, \quad k \in \{2, 4, 6, \dots\},\ d \geq k. \end{align}

For $\alpha \downarrow 0$ , we have that $\Lambda \sim Beta(\alpha, \alpha)$ converges to a discrete distribution with ${\mathbb{P}}(\Lambda = 0) = {\mathbb{P}}(\Lambda = 1) = 1/2$ . The resulting FGM copula becomes the extreme positive dependence copula $C^{EPD}$ that we will describe in Theorem 8. When $\alpha \uparrow \infty$ , we have that $\Lambda \sim Beta(\alpha, \alpha)$ converges to a discrete distribution with ${\mathbb{P}}(\Lambda = 1/2) = 1$ . Then ${\mathbb{E}}\left[(\Lambda - 1/2)^k\right] = 0$ for $k \in \{2, 3, \dots\}$ , and the corresponding FGM copula is the independence copula. One could define the beta-eFGM copula with the parametrization $\Lambda \sim Beta(1/\alpha, 1/\alpha)$ such that the dependence is monotone increasing with $\alpha$ , leading to a more intuitive ordering of the effect of $\alpha$ on the overall dependence. However, we opt to keep the current definition since the expressions for the copula parameters are more convenient.

4. Extreme points of eFGM copulas

In this section, we show that the parameters of any eFGM copula can be expressed as convex combinations of linearly independent parameters of eFGM copulas. More precisely, in every dimension $d \in \{2, 3, \dots\}$ , we find the convex hull for the parameters for the class of d-variate eFGM copulas. We call each vertex in the hull an extreme point, and we seek the extreme points of $\mathcal{T}_d$ , which are the extreme points from the set of inequalities in (7). See, for example, [52, Section 18] for the relationship between convex hulls and extreme points, and [59] for details on extremal rays of convex cones.

Theorem 4. Mixtures of eFGM copulas are also eFGM copulas.

Proof. Consider a vector of probabilities $(\lambda_1, \dots, \lambda_{n})$ such that $\lambda_j \geq 0$ for $j \in \{1, \dots, n\}$ and $\lambda_1 + \dots + \lambda_{n} = 1$ , with the notation that $\theta_{k, j}$ is the k-dependence parameter for the jth copula in the mixture of eFGM copulas. Consider the parameters $(\theta_{2, 1}, \dots, \theta_{d, 1}), \dots, (\theta_{2, n}, \dots, \theta_{d, n})$ . A convex combination of eFGM copulas has parameters $\theta_k = \sum_{j = 1}^{n} \lambda_j \theta_{k, j}$ , $k \in \{2, \dots, d\}$ . One must then verify that the constraints in (7) remain satisfied; indeed,

\begin{align*} 1+\sum_{k=2}^{d}\sum_{1\leq j_{1}<\dots <j_{k}\leq d}\theta_k\varepsilon _{j_{1}}\dots \varepsilon _{j_{k}}&= 1+\sum_{k=2}^{d}\sum_{1\leq j_{1}<\dots <j_{k}\leq d}\sum_{m = 1}^{n} \lambda_m\theta_{ k, m}\varepsilon _{j_{1}}\dots \varepsilon _{j_{k}}\nonumber\\ &=\sum_{j = 1}^{n} \lambda_j \left(1+\sum_{k=2}^{d}\sum_{1\leq j_{1}<\dots <j_{k}\leq d}\theta_{k, m}\varepsilon _{j_{1}}\dots \varepsilon _{j_{k}}\right) \end{align*}

for $\{\varepsilon_{1},\dots, \varepsilon_{d}\} \in \{-1,1\}^d$ . Since $\lambda_j \geq 0$ for $j \in \{1, \dots, n\}$ , every summand above satisfies (7), implying that $(\theta_2, \dots, \theta_d)$ also satisfies (7).

Remark 4. Theorem 4 can be understood probabilistically as follows. Let $\boldsymbol{U}_j$ be a random vector whose dependence structure can be expressed as an eFGM copula with parameters $(\theta_{2, 1}, \dots, \theta_{d, 1})$ and $(\theta_{2, 2}, \dots, \theta_{d, 2})$ . Construct a new random vector $\boldsymbol{U}_3$ with $F_{\boldsymbol{U}_3}(u) = (1 - \alpha) F_{\boldsymbol{U}_1}(u) + \alpha F_{\boldsymbol{U}_2}(u)$ , for $0\leq \alpha \leq 1$ and $u \in [0, 1]$ . Then the dependence structure of $\boldsymbol{U}_3$ is given by an eFGM copula with parameters $\theta_{k, 3} = (1 - \alpha) \theta_{k, 1} + \alpha \theta_{k, 2}$ , for $k \in \{2, \dots, d\}$ .

The authors of [20] show that the class $\mathcal{N}_d$ is a convex polytope generated from a finite number of extremal points. They also provide expressions for these extremal points, but let us first set up some notation. Denote by $n_d$ the number of extremal points in $\mathcal{N}_d$ . In Corollary 4.6 of [20], the authors show that

$$n_d = \begin{cases} (d + 1)^2/4, & d \text{ is odd},\\ d^2/4 + 1, & d \text{ is even}.\end{cases}$$

Let $(j_1^\wedge, j_2^\vee) = ((d - 1)/2, (d + 1)/2)$ if d is odd, and $(d/2 - 1, d/2 + 1)$ if d is even. Furthermore, define the one-to-one correspondence between the index $j\in \{1, \dots, n_d\}$ and every combination of the pairs $(j_1, j_2) \in \{0, 1, \dots, j_1^\wedge\}\times \{j_2^\vee, j_2^\vee + 1, \dots, d\}$ with

(17) \begin{equation} j = \begin{cases} 1 + j_1 + (j_1^\wedge + 1) (j_2 - j_2^{\vee}), & \substack{j_1 \in \{0, \dots, j_1^\wedge\}, \\ j_2 \in \{j_2^\wedge, \dots, d\},}\\[4pt] d^2/4 + 1, & d \text{ is even}. \end{cases}\end{equation}

The following proposition provides the expressions for the PMFs for the degenerate or two-point distributions that make up the extremal points of $\mathcal{N}_d$ . In particular, (17) provides a correspondence between the index $j\in\{1, \dots, n_d\}$ and the pair $(j_1, j_2) \in \{0, 1, \dots, j_1^\wedge\}\times \{j_2^\vee, j_2^\vee + 1, \dots, d\}$ .

Proposition 1. The PMFs corresponding to the extreme points of $\mathcal{N}_d$ are given by

(18) \begin{equation} {\mathbb{P}}(N_{jd} = k) = \begin{cases} \frac{j_2 - d/2}{j_2 - j_1}, & k = j_1,\\[4pt] \frac{d/2 - j_1}{j_2 - j_1}, & k = j_2,\\[4pt] 0, & {otherwise,} \end{cases} \end{equation}

for every combination of $j_1 \in \{0, 1, \dots, j_1^\wedge\}$ and $j_2 \in \{j_2^\vee, j_2^\vee + 1, \dots, d\}$ . If d is even, the set of extreme points also contains the one-point distribution at $k = d/2$ .

Proof. See [20, Proposition 4.5].

It follows that any PMF in $\mathcal{N}_d$ can be expressed as a convex combination of the $n_d$ extreme points, that is,

(19) \begin{equation} {\mathbb{P}}\left(N_d = k\right) = \sum_{j = 1}^{n_d} \lambda_j {\mathbb{P}}\left(N_{jd} = k\right)\!, \quad k = 0, \dots, d,\end{equation}

where $N_{jd}$ corresponds to the RV associated with the jth extreme point of $\mathcal{N}_d$ , $j \in \{1, \dots, n_d\}$ , and $(\lambda_1, \dots, \lambda_{n_d})$ is a vector such that $\lambda_m \geq 0$ , $m \in \{1, \dots, n_d\}$ , and $\lambda_1 + \dots + \lambda_{n_d} = 1$ .

From each extreme point in (18), one can extract an associated extreme point for the dependence parameters of eFGM copulas; the solutions solve the dual problem of finding the extreme points in the set of inequalities in (7). In other words, there is a one-to-one correspondence between the extreme points of $\mathcal{N}_d$ and the extreme points of $\mathcal{T}_d$ . For $d = 2$ , the eFGM copula parameters corresponding to the extreme points of $\mathcal{T}_2$ are $-1$ and 1.

Figure 1(a) presents the coordinates $(p_0, p_1, p_2)$ of the extremal points of $\mathcal{N}_3$ ; the last value is not free since $p_3 = 1 - p_0 - p_1 - p_2$ . The convex hull forms a geometric kite on a plane. In Figure 1(b), we present the kite from (a) along with the kite’s inscribed circle, which has a radius of $1/3$ . In Figure 1(c) we present the extreme points of $\mathcal{T}_3$ , which also form a kite, but not a scaled version of the kite in (b). The coordinates associated to the extreme point $(\theta_2, \theta_3) = (0, 1)$ are presented in bold.

Figure 1. Convex hull of admissible eFGM copula parameters for three dimensions.

Figure 2(a) presents the convex hull of PMFs $(p_0, \dots, p_4)$ generated by the extreme points of $\mathcal{N}_4$ . Each coordinate represents a point $(p_0, p_1, p_2, p_3)$ since $p_4$ is not free; we have $p_4 = 1 - p_1 - p_2 - p_3$ . We represent the coordinates inside a tesseract defined by the Cartesian product $[0, 1]^4$ , to represent the coordinates in four dimensions. (We use a tesseract designed by Claude Bragdon; see [53] for details.) The convex polytope generated by the extreme points of $\mathcal{N}_4$ corresponds to a pyramid with a kite base, called a kite pyramid, with the most negative dependence case at the apex. Figure 2(b) presents the extremal points of $\mathcal{T}_4$ within the cube $[-1, 1]^3$ , given with dotted lines for scale. The shape of $\mathcal{T}_4$ is also a kite pyramid. Finally, we present the dependence parameters associated to the extreme points of $\mathcal{N}_{10}$ in Table 1. Rows 22 and 27 (in bold) represent respectively the extremal positive and negative dependence within the class of 10-variate FGM copulas, as we will show in (24).

Table 1. Extremal points of the set of parameters $\mathcal{T}_{10}$ associated to $\mathcal{N}_{10}$ .

Figure 2. Convex hull of admissible eFGM copula parameters for four dimensions.

Remark 5. The representation as a convex combination of extreme points is not unique, meaning there may be different sets of weights $\lambda_1, \dots, \lambda_{n_d}$ which yield the same eFGM copula parameters. For example, let $d = 3$ ; then one recovers the independence parameter vector (0, 0) with $2^{-1} (0, 1) + 2^{-1}(0, -1)$ and $3/4({-}1/3, 0) +4^{-1}(1, 0)$ .

Exact methods of computing the extremal points of eFGM copula parameters not only are useful for understanding the properties of eFGM copulas, but also provide methods for solving problems within applications. For instance, we have from (7) that eFGM copula parameters must satisfy a set of $2^d$ inequalities. The computational cost of verifying these inequalities can be prohibitive; indeed, the complexity of doing so is exponential in time. However, one may compute the extremal points of the eFGM copula in quadratic time. Then verifying whether a set of parameters lies within the convex hull generated by the extremal points of eFGM copula parameters is a linear feasibility problem, which can be solved via linear programming (see, for instance, [10, Chapter 29] for an introduction to linear programming). In particular, if one can find a vector $(\lambda_1, \dots, \lambda_{n_d})$ such that $\theta_k = \sum_{j = 1}^{n_d}\lambda_j \theta_{k, j}$ for $k \in \{2, \dots, d\}$ , then $(\theta_2, \dots, \theta_d) \in \mathcal{T}_d$ . Hence, the extreme points of eFGM copula parameters enable a membership testing algorithm in polynomial time, as opposed to exponential time.

5. Extendability of eFGM copulas

We now address the question of the extendability of eFGM copulas. We start by studying the class of trivariate eFGM copula parameters that can be extended to k-variate eFGM copulas for $k > 3$ . In that case, we may visualize the admissible set of parameters $(\theta_2, \theta_3)$ graphically. We then present a characterization of infinitely extendable eFGM copulas.

Let $\mathcal{T}_{d, k}$ be the subset of $\mathcal{T}_d$ that can be extended to a valid element of $\mathcal{T}_{k}$ , where $d \leq k$ (in our analysis of extendability, we will always assume that $d \leq k$ ). In other words, for any vector $(\theta_2, \dots, \theta_d) \in \mathcal{T}_{d, k}$ , there exists $(\theta_{d+1}, \dots, \theta_k)$ such that $(\theta_2, \dots, \theta_d, \theta_{d+1}, \dots, \theta_k) \in \mathcal{T}_k$ . We have that

(20) \begin{equation} \mathcal{T}_{d, k} \subset \mathcal{T}_{d, k-1} \subset \dots \subset \mathcal{T}_{d, d+1}\subset \mathcal{T}_{d, d} = \mathcal{T}_d,\end{equation}

which follows from the observation that a k-variate FGM copula evaluated at $C(u_1, \dots, u_d, 1, \dots, 1)$ is a d-variate FGM copula. Our interest lies in studying the relationship between $\mathcal{T}_{d, k}$ and $\mathcal{T}_{k}$ —in particular, in studying the set of d-variate copula parameters that can be extended to a k-variate copula but not to a $(k+1)$ -variate copula, i.e. the set $\mathcal{T}_{d, k} \setminus \mathcal{T}_{d, k+1}$ .

Let us examine what we mean by the previous statement, by observing $\mathcal{T}_{3, k}$ for $k \in \{3, 4, \dots, 8\}$ . Note that one can find $\mathcal{T}_{3}$ in Figure 1(c). The convex hull generated by $\mathcal{T}_{3, 4}$ corresponds to Figure 2(b) but ignoring the coordinate $\theta_4$ . The subset $\mathcal{T}_{3} \setminus \mathcal{T}_{3, 4}$ , represented in diagonal red lines in Figure 3, corresponds to parameters $(\theta_2, \theta_3) \in \mathcal{T}_3$ , but for which there does not exist a $\theta_4$ such that $(\theta_2, \theta_3, \theta_4) \in \mathcal{T}_4$ . Equivalently, the orange checkerboard area corresponds to pairs $(\theta_2, \theta_3)$ that can be extended to 4-variate eFGM copulas but not to 5-variate eFGM copulas. In general, we obtain the set $\mathcal{T}_{3, d}$ from the convex hull generated by the parameters $(\theta_2, \theta_3)$ from the extremal points in $\mathcal{T}_{d}$ .

Figure 3. Extendability of trivariate eFGM copulas.

The set $\mathcal{T}_{d, \infty}$ corresponds to the d-variate eFGM copula parameters that are infinitely extendable. That is, if $(\theta_2, \dots, \theta_d) \in \mathcal{T}_{d, \infty}$ , then, for any dimension $k > d$ , there exist parameters $(\theta_{d+1}, \dots, \theta_k)$ such that $(\theta_2, \dots, \theta_d, \theta_{d + 1}, \dots, \theta_k) \in \mathcal{T}_{k, \infty} \subset \mathcal{T}_k$ . Since infinitely extendable eFGM copulas admit a de Finetti representation, it is possible to specify $\mathcal{T}_{3, \infty}$ . It is proved in Appendix B that the area inscribed between the functions $\theta_3 = \pm \theta_2 (1 - \theta_2)$ for $0\leq \theta_2 \leq 1$ generates the pairs of admissible $\mathcal{T}_{3, \infty}$ ; we also represent this area in violet dots in Figure 3.

We have characterized the set of copula parameters $(\theta_2, \theta_3)$ that are infinitely extendable. This brings up the more general problem of infinite extendability; see [39] for an overview. In particular, the following theorem characterizes the class of infinitely extendable eFGM copulas.

Theorem 5. Let $\boldsymbol{U}$ be a d-variate infinitely extendable random vector with eFGM dependence. Then there exists an RV $\Lambda$ such that $C = F_{\boldsymbol{U}}$ admits the representation

(21) \begin{equation} C(\boldsymbol{u}) = \int_{0}^1 \prod_{k = 1}^d \left\{(1-\lambda) u_k(2-u_k) + \lambda u_k^2\right\}{\, \mathrm{d}} F_{\Lambda}(\lambda). \end{equation}

Proof. From the construction of eFGM random vectors in (8), we have that for $\boldsymbol{U}$ to be infinitely extendable, it is necessary for $\boldsymbol{I}$ to also be infinitely extendable. It follows that $\boldsymbol{I}$ admits a de Finetti representation as in (11); hence one may write

$$C(\boldsymbol{u}) = {\mathbb{E}}\left[ \prod_{k = 1}^d u_k \left(1 + ({-}1)^{I_k} \overline{u}_k\right)\right] = \int_{0}^1 {\mathbb{E}}\left[ \prod_{k = 1}^d u_k \left(1 + ({-}1)^{I_k} \overline{u}_k\right) \vert \Lambda = \lambda\right] {\, \mathrm{d}} F_{\Lambda}(\lambda).$$

From the conditional independence of $\boldsymbol{I}$ on $\Lambda$ and the linearity of expectation, the copula becomes

$$C(\boldsymbol{u}) = \int_{0}^1 \prod_{k = 1}^d \left(u_k \left(1 + {\mathbb{E}}\left[({-}1)^{I_k}\vert \Lambda = \lambda\right] \overline{u}_k\right)\right){\, \mathrm{d}} F_{\Lambda}(\lambda).$$

Evaluating the (univariate) expectation and simplifying yields the desired result.

Let us interpret the results of Theorem 5. Conditionally on $\Lambda$ , the copula C can be expressed as the product of univariate CDFs, which is not surprising thanks to de Finetti’s representation theorem. Each univariate CDF corresponds to a mixture of the CDFs of $U_{[1]}$ and $U_{[2]}$ . For a given value of $\Lambda = \lambda$ , we have that if $\lambda < 0.5$ ( $\lambda > 0.5$ ), then the CDF of each margin is more likely to be that of $U_{[1]}$ (of $U_{[2]}$ ). The random mixture in the CDF of $\boldsymbol{U}$ in (21) is the same as the one in the PMF of $\boldsymbol{I}$ in (11), emphasizing the fact that the random vector $\boldsymbol{I}$ determines the dependence structure of the underlying FGM copula.

6. Dependence ordering

6.1. Supermodular order

In this section, we aim to compare the strength of dependence between two random vectors $\boldsymbol{V}$ and $\boldsymbol{V}'$ whose multivariate CDF is an eFGM copula. We will do so with dependence stochastic orders. The supermodular order is a valuable tool for comparing d-variate vectors of RVs with respect to the degree (strength) of dependence among their components. The supermodular order has been applied in a wide spectrum of fields, such as economics, applied probability, operations research, and statistics, in particular, to compare the dependence among risks within a portfolio of insurance policies or among assets in a financial institution’s portfolio. It is also used in the analysis of dependence properties and the estimation of copulas. In Sections 3.8 and 3.9 of [44], the authors provide an excellent introduction to the supermodular order. See also Section 6.3 of [13] and Section 9.A of [58]. Applications for ordering of actuarial risks can be found in Section 8.3 of [44] and Section 6.3 of [13]. The supermodular order allows one to identify the extremal positive dependence structures in the supermodular sense within a Fréchet class with fixed marginals, leading to the most positive dependence among random vectors. Also, one may derive results about other important orders from the supermodular order. Over the last two decades, there has been increasing interest in supermodular games in the fields of economics and operations research (see, for example, [2, 3, 60]). Supermodular functions and their mirror images, submodular functions, also appear in various branches of discrete mathematics and have numerous applications in computer science and optimization. As examples of results closer to those of the present section, the author of [62] gives sufficient and necessary conditions for the supermodular order of multivariate elliptical random vectors. In [4], the authors extend the ordering conditions to elliptical distributions that characterize the stronger supermodular ordering established for Gaussian distributions by [43]. As an application, they obtain several results on risk bounds in elliptical classes of risk models with fixed marginals. See also [8, 57, 61] for characterizations of the supermodular order for multivariate Marshall–Olkin exponential distributions, Gaussian copulas, and the family of Archimedean copulas.

The supermodular order is defined in terms of supermodular functions. A function $\phi\, :\, \mathbb{R}^{d}\rightarrow \mathbb{R}$ is said to be supermodular if

\begin{eqnarray*} &&\phi (x_{1},\ldots,x_{i}+\varepsilon ,\ldots,x_{j}+\delta ,\ldots,x_{d})-\phi (x_{1},\ldots,x_{i}+\varepsilon ,\ldots,x_{j},\ldots,x_{d}) \\ &\geq\ &\phi (x_{1},\ldots,x_{i},\ldots,x_{j}+\delta ,\ldots,x_{d})-\phi (x_{1},\ldots,x_{i},\ldots,x_{j},\ldots,x_{d})\end{eqnarray*}

holds for all $(x_1, \dots, x_d)\in \mathbb{R}^{d}$ , $1\leq i < j\leq d$ , and all $\varepsilon$ , $\delta >0$ . Examples of supermodular functions are $\phi(x_1,\dots,x_n) = x_1 + \dots + x_n$ , $\phi(x_1,\dots,x_n) = \min(x_1, \dots , x_n)$ , and $\phi(x_1,\dots,x_n) = h(x_1 + \dots + x_n)$ , where h is a convex function. Other examples are provided in Section 6.D of [40]. In economics, operations research, and machine learning, one is interested in optimizing a submodular function $-\phi$ , where $\phi$ is a supermodular function (see, for example, [3]).

Definition 2. (Supermodular order.) We say $(V_1, \dots, V_d)$ is smaller than $(V_1^{\prime}, \dots, V_d^{\prime})$ under the supermodular order, and we write $(V_1, \dots, V_d)\preceq _{sm}(V_1^{\prime}, \dots, V_d^{\prime})$ , if ${\mathbb{E}}\left\{\phi (V_1, \dots, V_d)\right\} \leq {\mathbb{E}}\left\{ \phi (V_1^{\prime}, \dots, V_d^{\prime})\right\} $ for all supermodular functions $\phi $ , given that the expectations exist.

The supermodular order satisfies the nine desired properties for dependence orders as mentioned in Section 3.8 of [44]. Ordering random vectors according to the supermodular order is desirable since it implies stochastic ordering results for the sum or functions of the components of those vectors of RVs.

6.2. Supermodular ordering within eFGM copulas

The following theorem from [7] uses the one-to-one correspondence between the family of d-variate FGM copulas and the family of d-variate symmetric Bernoulli distributions to characterize the supermodular order within the family of d-variate FGM copulas with the stochastic representation.

Theorem 6. Let $(I_1, \dots, I_d)$ and $(I_1^{\prime}, \dots, I_d^{\prime})$ be two symmetric multivariate Bernoulli distributed random vectors. Let $(U_1, \dots, U_d)$ and $(U_1^{\prime}, \dots, U_d^{\prime})$ be two random vectors constructed with the representation in Corollary 1, respectively using the random vectors $(I_1, \dots, I_d)$ and $(I_1^{\prime}, \dots, I_d^{\prime})$ . If $(I_1, \dots, I_d) \preceq_{sm} (I_1^{\prime}, \dots, I_d^{\prime})$ , then $(U_1, \dots, U_d) \preceq_{sm} (U_1^{\prime}, \dots, U_d^{\prime})$ .

From Theorem 6, if one wants to derive the supermodular ordering between $(U_1, \dots, U_d)$ and $(U_1^{\prime}, \dots, U_d^{\prime})$ , it suffices to establish the supermodular ordering between $(I_1, \dots, I_d)$ and $(I_1^{\prime}, \dots, I_d^{\prime})$ . We need to recall the definition of the convex order before establishing the supermodular ordering between the two vectors of symmetric Bernoulli exchangeable RVs $(I_1, \dots, I_d)$ and $(I_1^{\prime}, \dots, I_d^{\prime})$ .

Definition 3. (Convex order.) Let $X$ and $X^{\prime}$ be RVs with finite means. We say that $X$ is smaller than $X^{\prime}$ under the convex order, and we write $X\preceq_{cx}X'$ , if ${\mathbb{E}}[\varphi(X)] \leq {\mathbb{E}}[\varphi(X')]$ for all real convex functions $\varphi$ such that the expectations exist.

For the special case of constructions with mixtures presented in Section 3.3, we have the following ordering property.

Proposition 2. Consider two random vectors $(I_1, \dots, I_d)$ and $(I_1^{\prime}, \dots, I_d^{\prime})$ with PMFs constructed with mixtures as in (11) with respective mixing RVs $\Lambda$ and $\Lambda'$ . If $\Lambda \preceq_{cx} \Lambda'$ , it follows that $(I_1, \dots, I_d) \preceq_{sm} (I_1^{\prime}, \dots, I_d^{\prime})$ and $(U_1, \dots, U_d) \preceq_{sm} (U_1^{\prime}, \dots, U_d^{\prime})$ .

Proof. Given the representation in (11), we can use either Proposition 4.1.iii of [14] or Theorem 2.11 of [11] to deduce that $\Lambda \preceq_{cx} \Lambda'$ implies $(I_1, \dots, I_d) \preceq_{sm} (I_1^{\prime}, \dots, I_d^{\prime})$ . Then the second inequality follows from Theorem 6.

Example 3. Let $\Lambda \sim Beta(\alpha, \alpha)$ and $\Lambda' \sim Beta(\alpha', \alpha')$ ; then we have that $\Lambda \preceq_{cx}\Lambda'$ when $0 < \alpha' < \alpha < \infty$ (see Table 1.1 of [44]). When the representation in (11) is used with a beta RV, the dependence between the components of $(U_1, \dots, U_d)$ increases as $\alpha$ decreases and tends to 0.

Example 4. Let $Y \sim Gamma(1/\alpha, 1/\alpha)$ and $Y' \sim Gamma(1/\alpha', 1/\alpha')$ , with $0< \alpha < \alpha' < \infty$ . Let also $\Lambda = \exp({-}Y k)$ and $\Lambda' = \exp({-}Y' k')$ , where $k = \mathcal{L}_{Y}^{-1}(0.5)$ and $k' = \mathcal{L}_{Y'}^{-1}(0.5)$ . Then we have that $\Lambda \preceq_{cx} \Lambda'$ . Constructing $(I_1, \dots, I_d)$ and $(I_1^{\prime}, \dots, I_d^{\prime})$ with the representation in (11) and respective mixing RVs $\Lambda$ and $\Lambda'$ , we have $(I_1, \dots, I_d) \preceq_{sm} (I_1^{\prime}, \dots, I_d^{\prime})$ and $(U_1, \dots, U_d) \preceq_{sm} (U_1^{\prime}, \dots, U_d^{\prime})$ .

Since we have not seen the proof that $\Lambda \preceq_{cx} \Lambda'$ within the context of Example 4, we provide one in Appendix D. In the remainder of this section, we aim to identify the extremal negative and positive structures in the sense of the supermodular order within the family of eFGM copulas.

Theorem 7. Let $(I_1^-, \dots, I_d^-)$ be a vector of symmetric Bernoulli RVs with PMF

(22) \begin{equation} {\mathbb{P}}(I_1^- = i_1, \dots, I_d^- = i_d) = \begin{cases} (r + 1 - d/2) \binom{d}{r}^{-1}, & i_\bullet = r,\\ (d/2 - r) \binom{d}{r + 1}^{-1}, & i_\bullet = r + 1,\\ 0 & \text{otherwise}, \end{cases} \end{equation}

where $r \leq d/2 \leq r + 1$ . Also define $(I_1^+, \dots, I_d^+)$ as the vector of symmetric Bernoulli RVs with PMF

(23) \begin{equation} {\mathbb{P}}(I_1^+ = i_1, \dots, I_d^+ = i_d) = \begin{cases} 1/2, & i_\bullet \in \{0, d\},\\ 0 & \text{otherwise.} \end{cases} \end{equation}

For all vectors of exchangeable symmetric Bernoulli RVs $(I_1, \dots, I_d)$ , we have

$$(I_1^-, \dots, I_d^-) \preceq_{sm} (I_1, \dots, I_d) \preceq_{sm} (I_1^+, \dots, I_d^+).$$

Proof. The relationship in (22) follows from (7.22) of [28] with $\pi = 1/2$ , which identifies this PMF as the most negative dependence for exchangeable Bernoulli RVs. Theorem 7 of [21] shows that (22) is also the lower bound under the supermodular order. The joint PMF in (23) corresponds to the joint PMF when the exchangeable and symmetric Bernoulli RVs are comonotonic, which also coincides with the PMF derived from the Fréchet–Hoeffding upper bound with symmetric Bernoulli marginals.

We note respectively the extreme negative dependence (END) and extreme positive dependence (EPD) eFGM copulas such that the following holds:

(24) \begin{equation} \left(U_1^{END}, \dots, U_d^{END}\right) \preceq_{sm} \left(U_1, \dots, U_d\right) \preceq_{sm} \left(U_1^{EPD}, \dots, U_d^{EPD}\right)\!,\end{equation}

for all random vectors $(U_1, \dots, U_d)$ with CDF $F_{U_1, \dots, U_d} = C$ in the eFGM family of copulas. The EPD eFGM copula, provided in [7, Theorem 5], is recalled in the following theorem.

Theorem 8. The FGM copula associated with the vector of comonotonic RVs $(I_1^+, \dots, I_d^+)$ is the EPD eFGM copula $C^{EPD}$ . The dependence parameters are $\theta_k = (1 + ({-}1)^k)/2$ ; that is, $\theta_k = 1$ when k is even, and $\theta_k = 0$ when k is odd. The expression for $C^{EPD}$ is given by

(25) \begin{equation} C^{EPD}\left(u_1, \dots, u_d\right) = \prod_{j = 1}^{d} u_j \left(1 + \sum_{k = 1}^{\left\lfloor \frac{d}{2} \right\rfloor}\sum_{1\leq j_{1}<\cdots <j_{2k}\leq d} \overline{u}_{j_1}\cdots \overline{u}_{j_{2k}}\right)\!, \quad \left(u_1, \dots, u_d\right) \in [0,1]^{d}. \end{equation}

The case $j_1 = 0$ and $j_2 = d$ in (18) leads to the EPD FGM copula. The lower bound within the family of d-variate eFGM copulas under the supermodular order is defined in the following theorem; the proof is provided in Appendix C.

Theorem 9. The copula constructed with the vector of RVs $\left(I_1^-, \dots, I_d^-\right)$ is the END copula, denoted by $C^{END}$ . The dependence parameters $(\theta_2, \dots, \theta_d)$ for the END copula $C^{END}$ are given by

(26) \begin{equation} \theta_k = {}_2F_1\left({-}\left\lfloor \frac{d + 1}{2}\right\rfloor, -k, 2\left\lfloor \frac{d + 1}{2}\right\rfloor, 2\right) = \frac{(1 + ({-}1)^k)}{2}\frac{\Gamma(k+1)\Gamma\left(\frac{1}{2} -\left\lfloor \frac{d + 1}{2}\right\rfloor\right)}{2^k \Gamma\left(\frac{k}{2}+1\right) \Gamma\left(\frac{k+1}{2} - \left\lfloor \frac{d + 1}{2}\right\rfloor\right)}. \end{equation}

Corollary 4. An alternate representation for the k-dependence parameters in Theorem 9 for d odd is

$$\theta_k = \begin{cases} -\frac{1}{d}, & k = 2,\ d \geq 3,\\ \frac{3}{(d-2)d}, & k = 4,\ d \geq 7,\\ -\frac{15}{(d-4)(d-2)d}, & k = 6,\ d \geq 11, \end{cases} \qquad \theta_k = \prod_{l = 1}^{k/2} \frac{1 - 2l}{d-2l+2}, \quad d \geq 2k - 1,$$

while for d even the k-dependence parameters are

$$\theta_k = \begin{cases} -\frac{1}{d-1}, & k = 2,\ d \geq 4,\\[5pt] \frac{3}{(d-3)(d-1)}, & k = 4,\ d \geq 8,\\[5pt] -\frac{15}{(d-5)(d-3)(d-1)}, & k = 6,\ d \geq 12, \end{cases} \qquad \theta_k = \prod_{l = 1}^{k/2} \frac{1 - 2l}{d-2l+1}, \quad d \geq 2k.$$

Remark 6. The dependence structure within the random vector $(I_1^-, \dots, I_d^-)$ introduced in Theorem 9 corresponds to complete mixability; see [51] for details.

Remark 7. Table 2 presents the values of $(\theta_2, \dots, \theta_d)$ of the END eFGM copula, $C^{END}$ , for $d \in \{2, \dots, 12\}$ . Although the parameters follow some pattern, it is not obvious from the values that this corresponds to the dependence parameters inducing the most negative dependence within the family of eFGM copulas. We offer a few observations on the patterns exhibited by the parameters of the END eFGM copula. First, (26) gives the same values for consecutive values of $d \in \{3, 5, 7, \dots\}$ and $(d + 1) \in \{4, 6, 8, \dots\}$ . Then, we always have $\theta_k = 0$ for k odd. One also notices alternate signs for k even; that is, $\theta_k$ is negative for $k/2 \in \{1, 3, 5, \dots\}$ and positive for $k/2 \in \{2, 4, 6, \dots\}$ . Since $\theta_0 = 1$ and $\theta_1 = 0$ for every FGM copula, one notices that the magnitude of the END dependence parameters $\theta_k$ is symmetric, decreasing for $k < d/2$ and increasing again for $k > d/2$ . Also, $\theta_k \to 0$ as $d \to \infty$ for $k \neq d$ , and the k-dependence parameters $\theta_k$ depend on d.

Table 2. Extreme negative dependence copula parameters.

Remark 8. The term $(1 + ({-}1)^k)/2$ in (26) implies that $\theta_k = 0$ for $k \in \{3, 5, 7, \dots\}$ , which is also the case for the EPD eFGM copula. As noted in [7], the dependence parameters for odd indices k do not contribute to the overall strength of dependence.

7. Sampling and estimation

7.1. Sampling

In [7], an efficient stochastic sampling method is proposed based on the stochastic representation of FGM copulas. In Algorithm 1, we leverage the representation based on the convex set $\mathcal{N}_d$ from Subsection 3.1 to sample observations from eFGM copulas efficiently. Note that when the PMF of $N_d$ is an extreme point of $\mathcal{N}_d$ , sampling is faster since the vector of probabilities $(p_0, \dots, p_d)$ will have at most two non-zero values. Also, for subfamilies of eFGM copulas based on mixtures as in (11), one may sample $\tilde{N}_d$ from line 1 of Algorithm 1 by first sampling $\tilde{\Lambda}$ , then sampling $\tilde{N}_d$ from a binomial distribution with d trials and success probability $\tilde{\Lambda}$ .

Algorithm 1. Stochastic sampling method for eFGM copulas

7.2. Estimation difficulties with the FGM family of copulas

The main difficulty in estimating the parameters of an FGM copula is that they must respect the constraints in (7). For this reason, the method of moments is unlikely to provide a set of parameters that satisfy the $2^d$ constraints. The paper [48] makes an attempt to estimate the parameters of FGM copulas by estimating the parameters one at a time and using the simplex algorithm to constrain the valid parameter set after each parameter is estimated. However, this method does not scale well to high dimensions and will provide different parameter values if the order of parameter estimation changes. In the following subsection, we provide an algorithm that guarantees that the resulting parameters satisfy (7).

7.3. Maximum likelihood estimation

The likelihood of a set of $m_{obs}$ independent observations of identically distributed random vectors $\left(u_{1m}, \dots, u_{dm}\right)\!,$ for $m \in \{1, \dots, m_{obs}\}$ , for the stochastic representation of eFGM copulas is

(27) \begin{equation} L(\theta_2, \dots, \theta_d) = \prod_{m = 1}^{m_{obs}}\sum_{\{i_1, \dots, i_d\}\in \{0,1\}^d} f_{I_1, \dots, I_d}(i_1, \dots, i_d) \prod\limits_{l=1}^{d} \left[1 + ({-}1)^{i_l}(1 - 2u_{ml})\right].\end{equation}

Maximizing (27) is feasible but is computationally inconvenient since one needs to apply the system of constraints in (7). Another approach involves using the representation from Section 3.1, estimating the parameters $p_k$ , $k \in \{0, \dots, d\}$ , under the constraints $\sum_{k = 0}^d p_k = 1$ , $\sum_{k = 0}^d kp_k = d/2$ , and $p_k \geq 0$ , $k \in \{0, \dots, d\}$ . With this representation, one estimates $d-1$ parameters and the procedure admits a unique solution, but we have not found an efficient algorithm to perform this optimization.

Instead, we use the construction based on Section 4, which defines eFGM copula parameters as convex combinations of parameters from extreme points of the PMFs in $\mathcal{N}_d$ . The main advantage of this construction is that the likelihood is expressed as a finite mixture of $n_d$ points. We can use an expectation-maximization approach to optimize the likelihood, which lets us estimate parameters in higher dimensions than we have observed in the literature with FGM copulas. The disadvantage of this approach is that the solution using the convex combinations of extreme points is not unique, as stated in Remark 5. This non-identifiability is not an issue (from a modeling perspective) if we convert the estimated parameters back to the values of $(\theta_2, \dots, \theta_d)$ , although such a method may not be statistically efficient. Using the eFGM copula representation from Section 3.1, we find that the likelihood is

(28) \begin{equation} L(\theta_2, \dots, \theta_d) = \prod_{m = 1}^{m_{obs}}\sum_{k = 0}^{d}{\mathbb{P}}(N_d = k)\frac{1}{\binom{d}{k}} \sum_{\substack{\{i_1, \dots, i_d\} \in \{0, 1\}^d \\ i_\bullet = k}}\prod\limits_{l=1}^{d} \left[1 + ({-}1)^{i_l}(1 - 2u_{ml})\right].\end{equation}

If we replace (19) into (28), the likelihood becomes

(29) \begin{equation} \prod_{m = 1}^{m_{obs}}\sum_{k = 0}^{d}\sum_{j = 1}^{n_d} \lambda_j {\mathbb{P}}\left(N_{jd} = k\right)\frac{1}{\binom{d}{k}} \sum_{\substack{\{i_1, \dots, i_d\} \in \{0, 1\}^d \\ i_\bullet = k}}\prod\limits_{l=1}^{d} \left[1 + ({-}1)^{i_l}(1 - 2u_{ml})\right].\end{equation}

Rearranging (29) yields

(30) \begin{equation} \prod_{m = 1}^{m_{obs}}\sum_{j = 1}^{n_d} \lambda_j \sum_{k = 0}^{d} {\mathbb{P}}\left(N_{jd} = k\right)\frac{1}{\binom{d}{k}} \sum_{\substack{\{i_1, \dots, i_d\} \in \{0, 1\}^d \\ i_\bullet = k}}\prod\limits_{l=1}^{d} \left[1 + ({-}1)^{i_l}(1 - 2u_{ml})\right] = \prod_{m = 1}^{m_{obs}}\sum_{j = 1}^{n_d} \lambda_j \xi_{mj},\end{equation}

where

(31) \begin{equation} \xi_{mj} = \sum_{k = 0}^{d} {\mathbb{P}}\left(N_{jd} = k\right)\frac{1}{\binom{d}{k}} \sum_{\substack{\{i_1, \dots, i_d\} \in \{0, 1\}^d \\ i_\bullet = k}}\prod\limits_{l=1}^{d} \left[1 + ({-}1)^{i_l}(1 - 2u_{ml})\right];\end{equation}

this does not depend on the parameters $\lambda_j$ , $j \in \{1, \dots, n_d\}$ , so it can be computed once at the beginning of the optimization procedure. Using Lagrange multipliers to impose constraints on the parameters $\lambda_{j}$ , $j = \{1, \dots, n_d\}$ , the log-likelihood to maximize is

(32) \begin{equation} \mathcal{J}(\lambda_1, \dots, \lambda_{n_d}, \mu) = \sum_{m = 1}^{m_{obs}}\ln \left(\sum_{j = 1}^{n_d} \lambda_j \xi_{mj}\right) + \mu \left(\sum_{j = 1}^{n_d}\lambda_j - 1\right).\end{equation}

We find the Lagrange multiplier $\mu = -m_{obs}$ and

(33) \begin{equation} \sum_{j = 1}^{m_{obs}} \frac{\xi_{jt}}{\sum_{l = 1}^{n_d}\hat{\lambda}_j \xi_{jl}} = m_{obs} \Longrightarrow \hat{\lambda}_{t} = \frac{\sum_{j = 1}^{m_{obs}}\frac{\hat{\lambda}_t \xi_{jt}}{\sum_{l = 1}^{n_d} \hat{\lambda}_l \xi_{jl}}}{m_{obs}}, \quad t \in \{0, \dots, n_d\}.\end{equation}

In Algorithm 2, we propose an iterative algorithm to estimate the weights $\hat{\lambda}_1, \dots, \hat{\lambda}_{n_d}$ .

Algorithm 2. MLE estimation as a combination of extreme points

7.4. Simulation study: random parameters

To illustrate the estimation procedure, we perform a simulation study and attempt to estimate the corresponding parameters of eFGM copulas. We use Algorithm 1 to sample observations and Algorithm 2 to estimate the parameters $\lambda_j$ , $j \in \{1, \dots, n_d\}$ . However, we compare the resulting values of $\theta_k$ , $k \in \{2, \dots, d\}$ , to identify unique parameters.

In this study, we consider estimation based on known uniform margins; in cases with unknown margins, one should compute pseudo-observations based on the ranks of the empirical distribution function (using the semiparametric method of [23] or information from margins of [29]). We consider dimension $d = 10$ and sample multivariate observations $(u_{1m}, \dots, u_{10m})$ for $m \in \{1, \dots, 10\ 000\}$ . We then estimate the parameters $\theta_k$ , $k \in \{2, \dots, d\}$ . To the best of our knowledge, this is the first study estimating the parameters of FGM copulas for 10 dimensions, since both the stochastic representation and the exchangeability assumption simplify the parameter space.

We repeat the simulation and estimation 100 times with the same set of randomly generated parameters (but satisfying (7)) and present the results in Figure 4. The coordinate represents the true parameter, and the box-plot presents the range of estimates across the 100 replications. We present the main estimation diagnostic statistics in Table 3. Even when the true value of the parameters is close to zero, there is little variation in the parameter estimates. For example, the value of $\theta_2$ is $0.0667$ , which induces weak dependence (Spearman’s correlation coefficient between each pair of marginals, that is, $\rho_{S}(U_{j_1}, U_{j_2})$ for $1 \leq j_1 < j_2 \leq d$ , is only $0.0667/3$ ), but the interquartile range is only 0.006, making the estimates significantly different from 0, on an empirical basis. Only the parameter $\theta_{10}$ has a real value outside of the interquartile range of estimated values. This is not surprising: estimation of $\theta_5$ is based on $10!/5!/5! = 252$ different 5-tuples for each observation, while $\theta_{10}$ is based on a single 10-tuple. However, as discussed in [7], the k-dependence parameters for k close to d have less impact on the overall dependence: for the multivariate extensions of Spearman’s rho presented in [46], the contribution of 10-dependence parameters is $1/3^{10}$ , while that of 2-dependence parameters is $1/3^2$ .

Table 3. Estimation statistics for the simulation study.

Figure 4. Box-plot of estimates for the simulation study.

One can obtain a more accurate estimate of the parameter $\theta_{10}$ by increasing the dimension. For instance, if we consider $d = 13$ , then we have $13!/10!/3!$ different 10-tuples for each observation; hence the estimate of $\theta_{10}$ is more accurate. We present the box-plot of the estimation for $d = 13$ in Figure 5. Note that $\theta_{10}$ is accurately estimated, as are the first values $(\theta_2, \dots, \theta_7)$ . We conclude that the most important parameters are adequately estimated.

Figure 5. Box-plot of estimates for the simulation study with $d = 13$ .

7.5. Simulation study: extremal points

In the first simulation study, we considered an arbitrary set of parameters for $(\theta_2, \dots, \theta_d)$ . In this study, we will consider a set of parameters corresponding to an extremal point of $\mathcal{T}_d$ . In this context, the parameters $(\lambda_1, \dots, \lambda_{n_d})$ are identifiable from Algorithm 2. We let the dimension d vary in $\{5, 10, 15, 20\}$ and the number of observations $m_{obs}$ vary in $\{100, 500, 1000\}$ . For a fixed d, we sample observations from the parameters generated by selecting $j = 1$ from the extremal points in (17). We are interested in the ability of our algorithm to recover the correct extremal point. In Figure 6, we present the estimated parameter associated with the extremal point $j = 1$ . We note that for $d = 20$ , we have $n_d = 101$ ; the panel $m_{obs} =100$ , $d = 20$ , is therefore overdetermined, yet the algorithm recovers non-zero values of $\widehat{\lambda}_1$ for over 75% of the replications. As expected, increasing the number of observations generally increases the frequency of identifying the correct extremal point.

Figure 6. Histograms of the predicted $\widehat{\lambda}_1$ within the simulation study.

One drawback of the algorithm is that convergence of (29) may take many steps. A faster algorithm would attempt to formulate the problem as a convex optimization problem to leverage algorithms constructed for that purpose. Future works will involve studying the asymptotic or non-asymptotic properties of the estimator generated by Algorithm 2 and scaling the algorithm to higher dimensions.

8. Conclusion

In this paper, we have considered the class of eFGM copulas, including their constructions and properties. Thanks to the one-to-one correspondence between FGM copulas and symmetric multivariate Bernoulli RVs, we can leverage the extensive literature on symmetric exchangeable Bernoulli RVs to study eFGM copulas. We obtain extreme points of eFGM copulas, then study the extendability of eFGM copulas. We compare random eFGM vectors under the supermodular order, which has important implications for practical applications of copulas.

As mentioned in the introduction, FGM copulas are the simplest case of Bernstein copulas, the latter in d dimensions having the expression

(34) \begin{equation} C(\boldsymbol{u}) = \sum_{j_1 = 0}^{m_1} \dots \sum_{j_d = 0}^{m_d}\alpha\left(\frac{j_1}{m_1}, \dots, \frac{j_d}{m_d}\right) P_{j_1, m_1}(u_1) \dots P_{j_d, m_d}(u_d),\end{equation}

where $P_{v, m}(u) = \binom{m}{v}u^v(1-u)^{m-v}$ and $\alpha$ is a d-variate copula, for $(m_1, \dots, m_d) \in \mathbb{N}^{d}$ . When $m_1 = \dots = m_d = 1$ , we obtain the FGM copula. Bernstein copulas being much more flexible than FGM copulas, it would be interesting to extend the results from the current paper to the family of exchangeable Bernstein copulas. The subfamily of exchangeable Bernstein copulas has the additional constraint that $\alpha$ be an exchangeable copula and $m_1 = \dots = m_d$ . That analysis will require more background; in particular, one would need to extend the results of [20] to exchangeable multinomial distributions (see [25] for a construction). For this reason, we defer this analysis to future works. However, we can identify one extremal point of exchangeable Bernstein copulas corresponding to the EPD. The analogue to the EPD FGM copula within the class of Bernstein copulas occurs when $\alpha$ is the comonotonic copula, leading to the EPD Bernstein copula

\begin{align*} C(\boldsymbol{u}) &= \sum_{j_1 = 0}^{m} \dots \sum_{j_d = 0}^{m}\min\left(\frac{j_1}{m}, \dots, \frac{j_d}{m}\right) P_{j_1, m}(u_1) \dots P_{j_d, m}(u_d)\\ &= \frac{1}{m + 1} \sum_{k = 0}^m \prod_{l = 1}^d I_{u_l}(k + 1, m + 1 - k),\end{align*}

where $I_{x}(a, b)$ is the regularized incomplete beta function.