A comparison of dependence function estimators in multivariate extremes

Abstract

Various nonparametric and parametric estimators of extremal dependence have been proposed in the literature. Nonparametric methods commonly suffer from the curse of dimensionality and have been mostly implemented in extreme-value studies up to three dimensions, whereas parametric models can tackle higher-dimensional settings. In this paper, we assess, through a vast and systematic simulation study, the performance of classical and recently proposed estimators in multivariate settings. In particular, we first investigate the performance of nonparametric methods and then compare them with classical parametric approaches under symmetric and asymmetric dependence structures within the commonly used logistic family. We also explore two different ways to make nonparametric estimators satisfy the necessary dependence function shape constraints, finding a general improvement in estimator performance either (i) by substituting the estimator with its greatest convex minorant, developing a computational tool to implement this method for dimensions $D\ge 2$ or (ii) by projecting the estimator onto a subspace of dependence functions satisfying such constraints and taking advantage of Bernstein–Bézier polynomials. Implementing the convex minorant method leads to better estimator performance as the dimensionality increases.

Appendices

Appendix 1: Nonparametric estimators

If ${\mathbf {Z}}=(Z_1,\ldots ,Z_D)^{\top }$ is a max-stable vector with distribution $G({\mathbf {z}})$ (1) and estimated margins ${\hat{G}}_d$ (3), then the variable

$$\begin{aligned} \kappa \left( \varvec{\omega }\right) =\underset{1\le d\le D}{\text {min}}\left[ \frac{-\text {log}\left\{ {\hat{G}}_{d}(Z_d)\right\} }{\omega _d}\right] , \quad \quad \varvec{\omega }\in S_{D}, \end{aligned}$$

(15)

has survival function

$$\begin{aligned} \Pr \left\{ \kappa \left( \varvec{\omega }\right)>t\right\}= & {} \Pr \left[ -\text {log}\left\{ {\hat{G}}_{1}(Z_1)\right\}> \omega _1t,\ldots ,\right. \\&\left. -\,\text {log}\left\{ {\hat{G}}_{d}(Z_D)\right\} >\omega _D t\right] \\= & {} \exp \left\{ -tA\left( \varvec{\omega }\right) \right\} , \end{aligned}$$

where $t\ge 0$. Thus, the random variable $\kappa \left( \varvec{\omega }\right)$ has standard exponential distribution with

$$\begin{aligned} \text {E}\left\{ \kappa \left( \varvec{\omega }\right) \right\}= & {} A(\varvec{\omega })^{-1}, \quad \quad \text { and } \nonumber \\ \text {E}\left[ \hbox {log}\left\{ \kappa \left( \varvec{\omega }\right) \right\} \right]= & {} -\hbox {log}\left\{ A(\varvec{\omega })\right\} -\eta , \end{aligned}$$

(16)

where $\eta$ is the Euler–Mascheroni constant. The most popular nonparametric estimators of $A(\varvec{\omega })$ are based on the equations in (16). For instance, Pickands (1981) proposed to estimate the dependence function by the reciprocal of the sample mean of $\kappa \left( \varvec{\omega }\right)$, i.e.

$$\begin{aligned} {\hat{A}}_n^{\text {RP}}(\varvec{\omega })=\left\{ \frac{1}{M}\sum _{i=1}^{M}\kappa _i\left( \varvec{\omega }\right) \right\} ^{-1},\quad \quad \varvec{\omega }\in S_{D}, \end{aligned}$$

(17)

where $\kappa _i\left( \varvec{\omega }\right)$ is the empirical counterpart of (15) using sample maxima ${\mathbf {m}}_i$. While ${\hat{A}}_n^{\text {RP}}(\varvec{\omega })$ provides a sensible estimation of $A(\varvec{\omega })$, it does not satisfy any of the conditions A1, A2 and A3. As such, Hall and Tajvidi (2000) proposed a variant which verifies A2 and A3, defined as

$$\begin{aligned} {\hat{A}}_n^{\text {RHT}}(\varvec{\omega })=\left\{ \frac{1}{M}\sum _{i=1}^{M}\underset{1\le d\le D}{\text {min}}\left( \frac{{\tilde{m}}^{-1}_{i;d}/{\bar{m}}^{-1}_{M;d}}{\omega _d}\right) \right\} ^{-1}{,}\quad \varvec{\omega }\in S_{D},\nonumber \\ \end{aligned}$$

(18)

where ${\bar{m}}^{-1}_{M;d}=M^{-1} \sum _{i=1}^{M} {\tilde{m}}^{-1}_{i;d}, \;d=1,\ldots ,D$, and ${\tilde{m}}$ as defined in (4). The second equation in (16) leads to another dependence function estimator

$$\begin{aligned} {\hat{A}}^*\left( \varvec{\omega }\right) =\exp \left[ -\eta -\frac{1}{M}\sum _{i=1}^{M}\hbox {log}\{\kappa _{i}\left( \varvec{\omega }\right) \}\right] ,\quad \quad \varvec{\omega }\in S_{D}.\nonumber \\ \end{aligned}$$

(19)

However, ${\hat{A}}^*\left( \varvec{\omega }\right)$ does not satisfy any of the constraints A1, A2, A3 either. Capéràa et al. (1997) further modified (16) introducing an estimator meeting the constraint A3 for $D=2$; its extension to arbitrary dimension, provided by Zhang et al. (2008), see also Gudendorf and Segers (2011), corresponds to the following estimator

$$\begin{aligned}&{\hat{A}}_n^{\text {RMCFG}}(\varvec{\omega })\nonumber \\&=\exp \left[ \hbox {log}\left\{ A^*(\varvec{\omega })\right\} - \sum _{d=1}^D \delta _d(\varvec{\omega }) \hbox {log}\left\{ A^*({\mathbf {e}}_d)\right\} \right] ,\nonumber \\&\varvec{\omega }\in S_{D}, \end{aligned}$$

(20)

where the $\delta _d, d=1,\ldots ,D$, are continuous functions verifying $\delta _d({\mathbf {e}}_k)=1_{(d=k)}$ for all $d,k=1,\ldots ,D$. In particular, in the present paper we use $\delta _d(\varvec{\omega })=\omega _d$. See Gudendorf and Segers (2010) for a more complete explanation of the construction of the above estimators and Segers (2007) for a review of their properties.

Another type of nonparametric estimator of the dependence function can be based on the madogram (Matheron 1987; Cooley et al. 2006). The multivariate madogram is defined as (see Naveau et al. 2009) the expected value of the difference between the maximum and the mean of the rescaled variables $G({\mathbf {Z}})^{\varvec{\omega }^{-1}}$, i.e.

$$\begin{aligned}&\nu \left( \varvec{\omega }\right) =\text {E}\left[ \underset{1\le d\le D}{\text {max}}\left\{ G_{d}(z_d)^{1/\omega _d}\right\} -\frac{1}{D}\sum _{d=1}^D G_{d}(z_d)^{1/\omega _d}\right] , \nonumber \\&\quad \varvec{\omega }\in S_{D}. \end{aligned}$$

(21)

The estimator suggested by Marcon et al. (2017) consists of a slight modification of the estimator introduced by Naveau et al. (2009) for $D=2$ that satisfies the conditions A2 and A3 and applies to high dimensions. It is defined as

$$\begin{aligned} {\hat{A}}_n^{\text {MMD}}(\varvec{\omega })=\frac{{\hat{\nu }}_{M}(\varvec{\omega })+c(\varvec{\omega })}{1-{\hat{\nu }}_{M}(\varvec{\omega })-c(\varvec{\omega })},\quad \quad \varvec{\omega }\in S_{D}, \end{aligned}$$

(22)

where $c(\varvec{\omega })=D^{-1}\sum ^D_{d=1}\omega _d/(1+\omega _d)$ and ${\hat{\nu }}_{M}$ denotes the following estimator of the multivariate madogram (21):

$$\begin{aligned} {\hat{\nu }}_{M}\left( \varvec{\omega }\right)= & {} \frac{1}{M}\sum _{i=1}^{M}\left[ \underset{1\le d\le D}{\text {max}}\left\{ {\hat{G}}_{d}(m_{i;d})^{1/\omega _d}\right\} \right. \\&\left. -\frac{1}{D}\sum _{d=1}^D {\hat{G}}_{d}(m_{i;d})^{1/\omega _d}\right] . \end{aligned}$$

Finally, Cormier et al. (2014) proposed to estimate the Pickands dependence function through a graphical tool. Given the bivariate copula $C(u_1, u_2)$, consider the transformation

$$\begin{aligned} T_i= & {} \frac{\hbox {log}\left\{ {\hat{G}}(u_{i;1})\right\} }{\hbox {log}\left\{ {\hat{G}}(u_{i;1}){\hat{G}}(u_{i;2})\right\} },\\ Z_i= & {} \frac{\hbox {log}\left[ {\hat{C}}_n\left\{ {\hat{G}}(u_{i;1}) {\hat{G}}(u_{i;2})\right\} \right] }{\hbox {log}\left\{ {\hat{G}}(u_{i;1}){\hat{G}}(u_{i;2})\right\} }. \end{aligned}$$

where for each $i =1,\ldots , n, {\hat{C}}_n$ denote the empirical copula. The Pickands dependence function is approximated by fitting to the points $(t_1, z_1), \ldots , (t_n, z_n)$ a constrained B-spline smoothing of order $m = 3$ for a suitable sequence of k interior knots, $\beta _1, \ldots , \beta _{m+k}$ and B-spline basis $\epsilon _{1;m}, \ldots , \epsilon _{m+k;m}$, i.e.

$$\begin{aligned} {\hat{A}}_n^{\text {COBS}}= \sum _{j=1}^{m+k}{\hat{\beta }}_j\epsilon _{m:j}, \end{aligned}$$

(23)

where ${\hat{\beta }}_j$ is an estimate of $\beta _j$ for each $j =1,\ldots ,m + k$.

Appendix 2: Parametric models for the dependence function

Several parametric models to describe extremal dependence have been proposed in the literature. The logistic model, introduced by Gumbel (1960a, b), is one of the most widely used thanks to its simplicity. It has dependence function of the form

$$\begin{aligned} A(\varvec{\omega })=\left( \sum _{d=1}^D \omega _{d}^{1/\alpha } \right) ^{\alpha }, \quad \varvec{\omega }\in S_{D}, \end{aligned}$$

(24)

with dependence parameter $0<\alpha \le 1$. The strength of dependence decreases as $\alpha$ increases. In particular, the cases of independence and perfect dependence correspond to $\alpha =1$ and ${\alpha \downarrow 0}$, respectively. When $D=2$, the model (24) simplifies to $A(\omega )=\left\{ \left( 1-\omega \right) ^{1/\alpha }+\omega ^{1/\alpha }\right\} ^{\alpha }, \omega \in [0,1]$.

A limitation of the logistic model resides in its lack of flexibility as its dependence structure relies on only one parameter and it does not allow for asymmetry, i.e. the margins are exchangeable. An asymmetric generalisation was suggested by Tawn (1988) for dimension $D=2$ and extended to higher dimensions by Tawn (1990); see also Smith et al. (1990), Coles and Tawn (1991) and Stephenson (2009). The asymmetric logistic model has dependence function

$$\begin{aligned} A(\varvec{\omega })= \sum _{C\in S}\left\{ \sum _{d\in C}(\phi _{Cd}\omega _d)^{1/\alpha _C}\right\} ^{\alpha _C},\quad \varvec{\omega }\in S_{D}, \end{aligned}$$

(25)

with dependence parameters $0<\alpha _C \le 1$ and asymmetry parameters $\phi _{Cd}=0$ if $d\notin C, 0\le \phi _{Cd}\le 1$ for $d=1,\ldots ,D$, with $\sum _{C \in S} \phi _{Cd}=1$, where S is the set of all non-empty subsets of ${\mathcal {D}}=\{1,\ldots ,D\}$. When $D=2$, the model (25) simplifies to $A(\omega )=(\phi _2-\phi _1)\omega +1-\phi _2+\left[ \left\{ \phi _2(1-\omega )\right\} ^{1/\alpha }+(\phi _1\omega )^{1/\alpha }\right] ^{\alpha }, \omega \in [0,1]$, with dependence parameter $0<\alpha \le 1$ and asymmetry parameters $0\le \phi _1, \phi _2\le 1$. The extent of asymmetry in the bivariate dependence structure decreases as the asymmetry parameters $\phi _1$ and $\phi _2$ approach unity. The symmetric case arises when $\phi _1=\phi _2=1$. Without much loss of flexibility, $\phi _1$ or $\phi _2$ could be fixed to 1 (Tawn 1988). Independence corresponds to $\phi _1=\phi _2=1$ and $\alpha =1$, whereas perfect dependence corresponds to $\phi _1=0$ or $\phi _2=0$ or $\alpha \downarrow 0$.

Appendix 3: Proof of equation (11)

If ${\mathbf {U}}=(U_1,\ldots ,U_D)^\top \sim C_1(u_{1},\ldots ,u_{D})$, where $C_1$ denotes the outer power Clayton copula with generator $\varphi (t)=(t^{\alpha }+1)^{-1}$, then the random vector ${\mathbf {X}}=\left\{ -1/\hbox {log}(U_1),\right. \left. \ldots ,-1/\hbox {log}(U_D)\right\} ^\top$ is distributed according to the same copula $C_1$ and its marginal distributions are unit Fréchet. According to Fougères (2004, p. 376), taking the sequences ${\mathbf {a}}_{n}=(n, \ldots ,n)^\top$ and ${\mathbf {b}}_{n}=(0, \ldots ,0)^\top$ defined in Sect. 1, the distribution of ${\mathbf {X}}$ is in the MDA of the logistic distribution (24). This means that, defining

$$\begin{aligned} (M_{n;1},\ldots ,M_{n;D})^\top= & {} \left[ \underset{1\le i\le n}{\text {max}}\left\{ -1/\hbox {log}(U_{i;1})\right\} ,\ldots ,\right. \\&\left. \underset{1\le i\le n}{\text {max}}\left\{ -1/\hbox {log}(U_{i;D})\right\} \right] ^\top , \end{aligned}$$

one has

$$\begin{aligned}&\Pr \left( \frac{M_{n;1}}{n}\le x_1,\ldots ,\frac{M_{n;D}}{n}\le x_D\right) \\&\quad \rightarrow \exp \left\{ -V\left( x_1,\ldots ,x_D\right) \right\} , \quad \quad n\rightarrow \infty , \end{aligned}$$

where $V\left( x_1,\ldots ,x_D\right) =\left( \sum _{d=1}^D x_d^{-1/\alpha }\right) ^\alpha , \alpha \in (0,1]$. Now, notice that any vector $\mathbf {U^*}=(U^*_1,\ldots ,U^*_D)^\top$ distributed as $\mathbf {U^*}\sim C_2({\mathbf {u}}^*)=C_{1}(u_{1}^{\phi _1},\ldots ,u_{D}^{\phi _D})\prod _{d=1}^Du_{d}^{1-\phi _d}$ has uniform margins and may be represented by setting $U_d^*=\max \left\{ U_d^{1/\phi _d},{\tilde{U}}_d^{1/(1-\phi _d)}\right\} , d=1,\ldots ,D$, with $(U_1,\ldots ,U_D)^T\sim C_1(u_1,\ldots ,u_D)$, independent of ${\tilde{U}}_d \underset{\text { i.i.d.}}{\sim } \text {Unif}(0,1)$. Therefore, writing $M_{n;d}^*=\underset{1\le i\le n}{\text {max}}\left\{ -1/\hbox {log}(U_{i;d}^*)\right\}$ and ${\tilde{M}}_{n;d}=\underset{1\le i\le n}{\text {max}}\left\{ -1/\hbox {log}({\tilde{U}}_{i;d})\right\} , d=1,\ldots ,D$, one obtains

$$\begin{aligned}&\Pr \left( \frac{M^*_{n;1}}{n}\le x_1,\ldots ,\frac{M^*_{n;D}}{n}\le x_D\right) \\&\quad =\Pr \left( \frac{M_{n;1}}{n}\le \frac{x_1}{\phi _1},\ldots ,\frac{M_{n;D}}{n}\le \frac{x_D}{\phi _D}\right) \\&\qquad \times \prod _{d=1}^D\Pr \left( \frac{{\tilde{M}}_{n;d}}{n}\le \frac{x_d}{1-\phi _d}\right) ,\\&\quad \rightarrow \exp \left[ -\left\{ V\left( \frac{x_1}{\phi _1},\ldots ,\frac{x_D}{\phi _D}\right) +\sum _{d=1}^D(1-\phi _d)x_d^{-1}\right\} \right] ,\\&\qquad \quad n\rightarrow \infty , \end{aligned}$$

which is a multivariate extreme-value distribution with exponent function

$$\begin{aligned} V\left( x_1,\ldots ,x_D\right) =\left\{ \sum _{d=1}^D \left( \frac{x_d}{\phi _d }\right) ^{-1/\alpha }\right\} ^\alpha +\sum _{d=1}^D(1-\phi _d)x_d^{-1}. \end{aligned}$$

Applying (2) we obtain

$$\begin{aligned} A(\varvec{\omega })=\left\{ \sum _{d=1}^D(\omega _d \phi _{d})^{1/\alpha }\right\} ^\alpha +\sum _{d=1}^D \omega _d(1-\phi _{d}), \end{aligned}$$

which is a special case of the asymmetric logistic model dependence function in (25).