Parametric dependence between random vectors via copula-based divergence measures

TL;DR

This paper introduces a copula-based framework for quantifying dependence among multiple random vectors, providing new estimators, theoretical properties, and practical applications in finance.

Contribution

It develops a novel dependence measure using copula-based Phi-divergences, including estimation procedures and asymptotic analysis for both parametric and semi-parametric models.

Findings

Explicit dependence coefficients for Gaussian copulas.

Asymptotic properties of estimators established.

Application demonstrated on financial data.

Abstract

This article proposes copula-based dependence quantification between multiple groups of random variables of possibly different sizes via the family of $Phi$-divergences. An axiomatic framework for this purpose is provided, after which we focus on the absolutely continuous setting assuming copula densities exist. We consider parametric and semi-parametric frameworks, discuss estimation procedures, and report on asymptotic properties of the proposed estimators. In particular, we first concentrate on a Gaussian copula approach yielding explicit and attractive dependence coefficients for specific choices of $Phi$, which are more amenable for estimation. Next, general parametric copula families are considered, with special attention to nested Archimedean copulas, being a natural choice for dependence modelling of random vectors. The results are illustrated by means of examples. Simulations and a real-world application on financial data are provided as well.

Tables

Table 1. Table 1: Common choices for the function Φ Φ \Phi .

$\Phi \left(t\right)$	Name	$\Phi \left(0\right)+{\Phi }^{\ast }\left(0\right)$
$t\log \left(t\right)$	mutual information	$0+\infty$
${\left(t-1\right)}^2$	Pearson distance	$1+\infty$
${\left(\sqrt{t}-1\right)}^2$	Hellinger distance	$1+1$
$\left|t-1\right|$	total variation distance	$1+1$
$-\left(t+1\right)\log \left(\frac{t+1}{2}\right)+t\log \left(t\right)$	Jensen-Shannon distance	$\log \left(2\right)+\log \left(2\right)$

Table 2. Table 2: Sample bias, variance and mean squared error of two estimators for the half Hellinger distance in a bivariate Gumbel copula with parameter θ 0 = 3 subscript 𝜃 0 3 \theta_{0}=3 , based on 1000 1000 1000 replications, a sample size n = 200 𝑛 200 n=200 , and M 200 = 100 , 1000 , 10 ​ 000 subscript 𝑀 200 100 1000 10 000 M_{200}=100,1000,10\hskip 1.42271pt000 . The true value equals 𝒟 ( t − 1 ) 2 / 2 = 0.20528 subscript 𝒟 superscript 𝑡 1 2 2 0.20528 \mathcal{D}_{(\sqrt{t}-1)^{2}}/2=0.20528 .

${𝐌}_{\mathrm{𝟐𝟎𝟎}}$	${\widehat{𝒟}}_{{\left(\sqrt{𝐭}-\mathrm{𝟏}\right)}^{\mathrm{𝟐}},\mathrm{𝟐𝟎𝟎},{𝐌}_{\mathrm{𝟐𝟎𝟎}}}/\mathrm{𝟐}$				${\widetilde{𝒟}}_{{\left(\sqrt{𝐭}-\mathrm{𝟏}\right)}^{\mathrm{𝟐}},\mathrm{𝟐𝟎𝟎},{𝐌}_{\mathrm{𝟐𝟎𝟎}}}/\mathrm{𝟐}$
	bias	var	mse		bias	var	mse
$100$	$-0.0109$	$0.1313$	$0.1314$		$0.0031$	$0.0039$	$0.0039$
$1000$	$-0.0078$	$0.0249$	$0.025$		$0.001$	$0.0008$	$0.0008$
$10000$	$0.0007$	$0.0139$	$0.0139$		$0.0008$	$0.0005$	$0.0005$

Table 3. Table 3: Empirical Monte Carlo variances n ​ Var ^ 𝑛 ^ Var n\widehat{\text{Var}} of estimated dependence coefficients for different sample sizes n 𝑛 n in different settings, based on 1000 1000 1000 replications.

$𝐧$	Mutual information			Half Hellinger distance
	Setting 1	Setting 2		Setting 1	Setting 2
$50$	$0.002$	$0.1311$		$0.0001$	$0.0124$
$200$	$0.0002$	$0.0734$		$1.19\cdot {10}^{-5}$	$0.0073$
$1000$	$1.62\cdot {10}^{-5}$	$0.0322$		$1.10\cdot {10}^{-6}$	$0.0042$
$5000$	$2.01\cdot {10}^{-6}$	$0.0217$		$2.35\cdot {10}^{-7}$	$0.0045$

Table 4. Table 4: P-values ( 20 ) for test of increased linkages between inter-regional equity indexes from pre-crisis to crisis and ( 21 ) for test of decreased linkages from crisis to post-crisis, using the mutual information or half Hellinger distance as dependence measure (NA = North America, SA = South America, EU = Europe, AS = Asia).

Financial contagion $p$-values
	Mutual information			Half Hellinger distance
	$p_{12}$	$p_{23}$		$p_{12}$	$p_{23}$
NA - SA	$1.1657\cdot {10}^{-7}$	$0.0008$		$2.8958\cdot {10}^{-7}$	$0.0008$
NA - EU	$8.1439\cdot {10}^{-5}$	$0.0039$		$0.0001$	$0.0052$
NA - AS	$0.0151$	$0.0672$		$0.0157$	$0.0655$
SA - EU	$1.7442\cdot {10}^{-7}$	$0.0012$		$7.3208\cdot {10}^{-7}$	$0.0016$
SA - AS	$0.0004$	$0.0114$		$0.0006$	$0.0126$
EU - AS	$0.0016$	$0.0388$		$0.0018$	$0.0317$

Equations

$$\mathcal{D}\big{(}T_{1}(\mathbf{X}_{1}),\dots,T_{k}(\mathbf{X}_{k})\big{)}=\mathcal{D}(\mathbf{X}_{1},\dots,\mathbf{X}_{k}),$$

$$\mathcal{D}\big{(}\mathbf{X}_{1},\dots,T_{i}(\mathbf{X}_{i}),\dots,\mathbf{X}_{k}\big{)}=\mathcal{D}(\mathbf{X}_{1},\dots,\mathbf{X}_{k}),$$

$$\underset{n\to \infty }{\lim }{\mathcal{D}}^{d_1,\dots ,d_k}({\mathbf{X}}_n)={\mathcal{D}}^{d_1,\dots ,d_k}(\mathbf{X})$$

$$\Phi (0)=\underset{\begin{array}{c}t\to 0\\ >\end{array}}{\lim }\Phi (t),\text{mbox}and{\Phi }^{\ast }(0)=\underset{t\to \infty }{\lim }\frac{\Phi (t)}{t}.$$

$$\mathcal{D}_{\Phi}=\int\Phi\bigg{(}\frac{d\mu_{C}^{\text{ac}}}{d(\mu_{C_{1}}\times\cdots\times\mu_{C_{k}})}\bigg{)}d(\mu_{C_{1}}\times\cdots\times\mu_{C_{k}})+\Phi^{*}(0)\mu_{C}^{\text{s}}(B),$$

$${\mathcal{D}}_{\Phi }\left({\mathbf{X}}_1,\dots ,{\mathbf{X}}_k\right)={\int }_{{\mathbb{I}}^q}\prod _{i=1}^k{c}_i({\mathbf{u}}_i)\Phi \left(\frac{c(\mathbf{u})}{{\prod }_{i=1}^k{c}_i({\mathbf{u}}_i)}\right)d\mathbf{u},$$

$$N_{\Phi }\circ {\mathcal{D}}_{\Phi }=\{\begin{array}{l}{\mathcal{D}}_{\Phi }/(\Phi (0)+{\Phi }^{\ast }(0))\text{if}\Phi (0)+{\Phi }^{\ast }(0)<\infty \\ N({\mathcal{D}}_{\Phi })\text{if}\Phi (0)+{\Phi }^{\ast }(0)=\infty ,\end{array}$$

$$\begin{split}&\mathcal{D}_{\Phi}(\mathbf{X}_{1},\dots,\mathbf{X}_{k},\mathbf{X}_{k+1})\\ &\hskip 28.45274pt=\int_{\mathbb{I}^{q+d_{k+1}}}\prod_{i=1}^{k+1}c_{i}(\mathbf{u}_{i})\Phi\Bigg{(}\frac{\widetilde{c}(\widetilde{\mathbf{u}})}{\prod_{i=1}^{k+1}c_{i}(\mathbf{u}_{i})}\Bigg{)}d\widetilde{\mathbf{u}}\\ &\hskip 28.45274pt=\int_{\mathbb{I}^{q}}\prod_{i=1}^{k}c_{i}(\mathbf{u}_{i})\int_{\mathbb{I}^{d_{k+1}}}c_{k+1}(\mathbf{u}_{k+1})\Phi\Bigg{(}\frac{\widetilde{c}(\mathbf{u},\mathbf{u}_{k+1})}{c_{k+1}(\mathbf{u}_{k+1})\prod_{i=1}^{k}c_{i}(\mathbf{u}_{i})}\Bigg{)}d\mathbf{u}_{k+1}d\mathbf{u}\\ &\hskip 28.45274pt\geq\int_{\mathbb{I}^{q}}\prod_{i=1}^{k}c_{i}(\mathbf{u}_{i})\Phi\Bigg{(}\int_{\mathbb{I}^{d_{k+1}}}\frac{\widetilde{c}(\mathbf{u},\mathbf{u}_{k+1})}{\prod_{i=1}^{k}c_{i}(\mathbf{u}_{i})}d\mathbf{u}_{k+1}\Bigg{)}d\mathbf{u}\\ &\hskip 28.45274pt=\mathcal{D}_{\Phi}(\mathbf{X}_{1},\dots,\mathbf{X}_{k}),\end{split}$$

$$\begin{split}\text{\scalebox{0.95}{$\mathcal{D}_{\Phi}(T_{1}(\mathbf{X}_{1}),\dots,\mathbf{X}_{k})=\int_{\mathbb{I}^{q}}c_{1}(\widetilde{\mathbf{u}}_{1})\prod_{i=2}^{k}c_{i}(\mathbf{u}_{i})\Phi\Bigg{(}\frac{c(\widetilde{\mathbf{u}}_{1},\mathbf{u}_{2},\dots,\mathbf{u}_{k})}{c_{1}(\widetilde{\mathbf{u}}_{1})\prod_{i=2}^{k}c_{i}(\mathbf{u}_{i})}\Bigg{)}d\widetilde{\mathbf{u}}=\mathcal{D}_{\Phi}(\mathbf{X}_{1},\dots,\mathbf{X}_{k})$}},\end{split}$$

$$\left(\prod _{i=1}^k{c}_{ni}\right)\Phi \left(\frac{c_n}{{\prod }_{i=1}^k{c}_{ni}}\right)\to \left(\prod _{i=1}^k{c}_i\right)\Phi \left(\frac{c}{{\prod }_{i=1}^k{c}_i}\right)$$

$$\mathbb{P}(X_1\le x_1)=\frac{2\arcsin (x_1)+\pi }{2\pi }\text{for}{x}_1\in [-1,1],$$

$$C(u_{1},u_{2})=\begin{cases}\frac{2(u_{1}+u_{2})-1}{4}+\frac{2(u_{1}+u_{2})-3}{4}\mathds{1}\Big{(}u_{1}+u_{2}>\frac{3}{2}\Big{)}&\mbox{if }(u_{1},u_{2})\in\Big{[}\frac{1}{2},1\Big{]}^{2}\\ \frac{2(u_{1}+u_{2})-1}{4}\mathds{1}\Big{(}u_{1}+u_{2}>\frac{1}{2}\Big{)}&\mbox{if }(u_{1},u_{2})\in\Big{[}0,\frac{1}{2}\Big{]}^{2}\\ \frac{1+u_{1}-\max\{1-u_{1},\frac{3}{2}-u_{2}\}}{2}&\mbox{if }(u_{1},u_{2})\in\Big{[}0,\frac{1}{2}\Big{]}\times\Big{[}\frac{1}{2},1\Big{]}\\ \frac{\min\{1+u_{1},\frac{3}{2}+u_{2}\}+u_{2}-\frac{3}{2}}{2}&\mbox{if }(u_{1},u_{2})\in\Big{[}\frac{1}{2},1\Big{]}\times\Big{[}0,\frac{1}{2}\Big{]},\end{cases}$$

$$C(u_{1},u_{2},u_{3},u_{4})=\exp\bigg{(}-\Big{(}\big{(}-\log(\min\{u_{1},u_{3}\})\big{)}^{\theta}+\big{(}-\log(\min\{u_{2},u_{4}\})\big{)}^{\theta}\Big{)}^{\frac{1}{\theta}}\bigg{)}$$

$$\mathbf{\Sigma }=\left(\begin{array}{cccc}{\mathbf{\Sigma }}_{11}& {\mathbf{\Sigma }}_{12}& \dots & {\mathbf{\Sigma }}_{1k}\\ {\mathbf{\Sigma }}_{12}^{\text{T}}& {\mathbf{\Sigma }}_{22}& \dots & {\mathbf{\Sigma }}_{2k}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{\Sigma }}_{1k}^{\text{T}}& {\mathbf{\Sigma }}_{2k}^{\text{T}}& \dots & {\mathbf{\Sigma }}_{kk}\end{array}\right)\in {\mathbb{R}}^{q\times q},$$

$${\mathcal{D}}_{\Phi }^{\mathcal{N}}(\mathbf{\Sigma })={\int }_{{\mathbb{R}}^q}\frac{\exp \left(-\frac{1}{2}{\mathbf{x}}^{\text{T}}{\mathbf{\Sigma }}_0^{-1}\mathbf{x}\right)}{(2\pi {)}^{q/2}{\prod }_{i=1}^k{\left|{\mathbf{\Sigma }}_{ii}\right|}^{1/2}}\Phi \left(\frac{{\prod }_{i=1}^k{\left|{\mathbf{\Sigma }}_{ii}\right|}^{1/2}}{{\left|\mathbf{\Sigma }\right|}^{1/2}}\exp \left(-\frac{1}{2}{\mathbf{x}}^{\text{T}}{\mathbf{B}}^{-1}\mathbf{x}\right)\right)d\mathbf{x},$$

$${\mathbf{\Sigma }}_0=\left(\begin{array}{cccc}{\mathbf{\Sigma }}_{11}& 0_{12}& \dots & 0_{1k}\\ 0_{12}^{\text{T}}& {\mathbf{\Sigma }}_{22}& \dots & 0_{2k}\\ ⋮& ⋮& \ddots & ⋮\\ 0_{1k}^{\text{T}}& 0_{2k}^{\text{T}}& \dots & {\mathbf{\Sigma }}_{kk}\end{array}\right)$$

$${\mathcal{D}}_{t\log (t)}^{\mathcal{N}}(\mathbf{\Sigma })=-\frac{1}{2}\log \left(\frac{\left|\mathbf{\Sigma }\right|}{{\prod }_{i=1}^k\left|{\mathbf{\Sigma }}_{ii}\right|}\right),$$

$${\mathcal{D}}_{(\sqrt{t}-1{)}^2}^{\mathcal{N}}(\mathbf{\Sigma })=2-2\frac{2^{q/2}{\left|\mathbf{\Sigma }\right|}^{1/4}}{{\left|{\mathbb{I}}_q+{\mathbf{\Sigma }}_0^{-1}\mathbf{\Sigma }\right|}^{1/2}{\prod }_{i=1}^k{\left|{\mathbf{\Sigma }}_{ii}\right|}^{1/4}}.$$

$$\begin{array}{rl}\left|{\mathbb{I}}_q+{\mathbf{\Sigma }}_0^{-1}\mathbf{\Sigma }\right|& =\left|{\mathbf{D}}_{\mathbf{\Sigma }}^{1/2}\right|\left|{\mathbb{I}}_q+{\mathbf{\Sigma }}_0^{-1}\mathbf{\Sigma }\right|\left|{\mathbf{D}}_{\mathbf{\Sigma }}^{-1/2}\right|=\left|{\mathbb{I}}_q+{\mathbf{D}}_{\mathbf{\Sigma }}^{1/2}{\mathbf{\Sigma }}_0^{-1}{\mathbf{D}}_{\mathbf{\Sigma }}^{1/2}\mathbf{R}{\mathbf{D}}_{\mathbf{\Sigma }}^{1/2}{\mathbf{D}}_{\mathbf{\Sigma }}^{-1/2}\right|\\ & =\left|{\mathbb{I}}_q+{\mathbf{D}}_{\mathbf{\Sigma }}^{1/2}{\mathbf{\Sigma }}_0^{-1}{\mathbf{D}}_{\mathbf{\Sigma }}^{1/2}\mathbf{R}\right|=\left|{\mathbb{I}}_q+{\mathbf{R}}_0^{-1}\mathbf{R}\right|,\end{array}$$

$$({\mathbf{R}}_{im}{)}_{jt}={\rho }_{ij,mt}=\text{Corr}\left(({\varphi }^{-1}\circ F_{ij})(X_{ij}),({\varphi }^{-1}\circ F_{mt})(X_{mt})\right),$$

$$\widehat{\mathbf{R}}_{n}=\begin{pmatrix}\widehat{\mathbf{R}}_{11}&\widehat{\mathbf{R}}_{12}&\cdots&\widehat{\mathbf{R}}_{1k}\\ \widehat{\mathbf{R}}_{12}^{\text{T}}&\widehat{\mathbf{R}}_{22}&\cdots&\widehat{\mathbf{R}}_{2k}\\ \vdots&\vdots&\ddots&\vdots\\ \widehat{\mathbf{R}}_{1k}^{\text{T}}&\widehat{\mathbf{R}}_{2k}^{\text{T}}&\cdots&\widehat{\mathbf{R}}_{kk}\end{pmatrix}\hskip 5.69046pt\text{with}\hskip 5.69046pt(\widehat{\mathbf{R}}_{im})_{jt}=\widehat{\rho}_{ij,mt}=\frac{\frac{1}{n}\sum_{\ell=1}^{n}\widehat{Z}_{ij}^{(\ell)}\widehat{Z}_{mt}^{(\ell)}}{\frac{1}{n}\sum_{\ell=1}^{n}\Big{(}\phi^{-1}\big{(}\frac{\ell}{n+1}\big{)}\Big{)}^{2}},$$

$$\widehat{Z}_{ij}^{(\ell)}=\phi^{-1}\left(\frac{n}{n+1}\widehat{F}_{ij}\Big{(}X_{ij}^{(\ell)}\Big{)}\right)$$

$$\begin{split}\frac{1}{n}\sum_{\ell=1}^{n}\widehat{Z}_{ij}^{(\ell)}&=\frac{1}{n}\sum_{\ell=1}^{n}\phi^{-1}\Big{(}\frac{\ell}{n+1}\Big{)}=0\\ \frac{1}{n}\sum_{\ell=1}^{n}\Big{(}\widehat{Z}_{ij}^{(\ell)}\Big{)}^{2}&=\frac{1}{n}\sum_{\ell=1}^{n}\bigg{(}\phi^{-1}\Big{(}\frac{\ell}{n+1}\Big{)}\bigg{)}^{2},\end{split}$$

$$({\mathbb{S}}^s,||\cdot |{|}_{\text{F}})\to (\mathbb{R},|\cdot |):\mathbf{\Sigma }\mapsto ({\mathcal{D}}_{\Phi }^{\mathcal{N}}\circ \phi )(\mathbf{\Sigma })$$

$$\sqrt{n}\left({\hat{\mathcal{D}}}_{\Phi ,n}^{\mathcal{N}}-{\mathcal{D}}_{\Phi }^{\mathcal{N}}(\mathbf{R})\right)\stackrel{d}{\to }\mathcal{N}(0,{\zeta }_{\Phi }^2)$$

$$\zeta_{\Phi}^{2}=2\text{Tr}\Big{(}\big{(}\mathbf{R}(\mathbf{M}_{\Phi}-\mathbf{D}_{\mathbf{M}_{\Phi}\mathbf{R}})\big{)}^{2}\Big{)},$$

$$\begin{split}\mathbf{M}_{\Phi}=\frac{1}{2}\Big{(}\mathbf{F}_{1}-\mathcal{D}_{\Phi}^{\mathcal{N}}(\mathbf{R})\mathbf{R}_{0}^{-1}&-\mathbb{E}_{\mathcal{N}(\mathbf{0},\mathbf{R})}\left[\alpha^{\prime}(\mathbf{X})\right]\left(\mathbf{R}^{-1}-\mathbf{R}_{0}^{-1}\right)\\ &+\mathbf{R}^{-1}\mathbb{E}_{\mathcal{N}(\mathbf{0},\mathbf{R})}\left[\alpha^{\prime}(\mathbf{X})\mathbf{X}\mathbf{X}^{\text{T}}\right]\mathbf{R}^{-1}-\mathbf{F}_{2}\Big{)},\end{split}$$

$$\begin{array}{rl}{\mathbf{F}}_1& ={\mathbf{R}}_0^{-1}\text{diag}\left({\mathbb{E}}_{\mathcal{N}(0,{\mathbf{R}}_0)}\left[\alpha (\mathbf{X}){\mathbf{X}}_1{\mathbf{X}}_1^{\text{T}}\right],\cdots ,{\mathbb{E}}_{\mathcal{N}(0,{\mathbf{R}}_0)}\left[\alpha (\mathbf{X}){\mathbf{X}}_k{\mathbf{X}}_k^{\text{T}}\right]\right){\mathbf{R}}_0^{-1},\\ {\mathbf{F}}_2& ={\mathbf{R}}_0^{-1}\text{diag}\left({\mathbb{E}}_{\mathcal{N}(0,\mathbf{R})}\left[{\alpha }^{\prime }(\mathbf{X}){\mathbf{X}}_1{\mathbf{X}}_1^{\text{T}}\right],\cdots ,{\mathbb{E}}_{\mathcal{N}(0,\mathbf{R})}\left[{\alpha }^{\prime }(\mathbf{X}){\mathbf{X}}_k{\mathbf{X}}_k^{\text{T}}\right]\right){\mathbf{R}}_0^{-1},\end{array}$$

$$k(\mathbf{X})=\frac{{\prod }_{i=1}^k{\left|{\mathbf{\Sigma }}_{ii}\right|}^{1/2}}{{\left|\mathbf{\Sigma }\right|}^{1/2}}\exp \left(-\frac{1}{2}{\mathbf{X}}^{\text{T}}\left({\mathbf{R}}^{-1}-{\mathbf{R}}_0^{-1}\right)\mathbf{X}\right).$$

$${\mathbf{M}}_{t\log (t)}=-\frac{1}{2}\left({\mathbf{R}}^{-1}-{\mathbf{R}}_0^{-1}\right)$$

Full text

For Property (A4), suppose that $(\mathbf{X}_{1},\dots,\mathbf{X}_{k},\mathbf{X}_{k+1})$ has copula density $\widetilde{c}$ with marginal copula density $c_{k+1}$ of $\mathbf{X}_{k+1}$. Put $\widetilde{\mathbf{u}}=(\mathbf{u},\mathbf{u}_{k+1})$ Then,

[TABLE]

where we used Jensen’s inequality. If $\mathbf{X}_{k+1}$ is independent from $(\mathbf{X}_{1},\dots,\mathbf{X}_{k})$, we have $\widetilde{c}(\mathbf{u},\mathbf{u}_{k+1})=c(\mathbf{u})c_{k+1}(\mathbf{u}_{k+1})$ and the equality holds. If $\Phi$ is strictly convex, the equality holds if and only if $\widetilde{c}(\mathbf{u},\mathbf{u}_{k+1})=A(\mathbf{u})c_{k+1}(\mathbf{u}_{k+1})\prod_{i=1}^{k}c_{i}(\mathbf{u}_{i})$, where $A(\mathbf{u})$ is a function of $\mathbf{u}$ not depending on $\mathbf{u}_{k+1}$ almost surely for almost every $\mathbf{u}$. Integrating this equality with respect to $\mathbf{u}_{k+1}$ gives $A(\mathbf{u})=c(\mathbf{u})/\prod_{i=1}^{k}c_{i}(\mathbf{u}_{i})$, i.e. $\widetilde{c}(\mathbf{u},\mathbf{u}_{k+1})=c(\mathbf{u})c_{k+1}(\mathbf{u}_{k+1})$.

Obviously, by definition, Property (A5) and hence also (A6) are fulfilled. Next, in the context of Property (A7), assume without loss of generality that $X_{11}$ gets transformed to $T_{11}(X_{11})$ for a strictly decreasing transformation $T_{11}$ and let $\widetilde{c}$ be the copula density of $(T_{1}(\mathbf{X}_{1}),\mathbf{X}_{2},\dots,\mathbf{X}_{k})$ with $T_{1}(\mathbf{X}_{1})=(T_{11}(X_{11}),X_{12},\dots,X_{1d_{1}})$. Then, $\widetilde{c}(\mathbf{u}_{1},\mathbf{u}_{2},\dots,\mathbf{u}_{k})=c(\widetilde{\mathbf{u}}_{1},\mathbf{u}_{2},\dots,\mathbf{u}_{k})$ with $\widetilde{\mathbf{u}}_{1}=(1-u_{11},u_{12},\dots,u_{1d_{1}})$ and $\widetilde{c}_{1}(\mathbf{u}_{1})=c_{1}(\widetilde{\mathbf{u}}_{1})$ with $\widetilde{c}_{1}$ the copula density of $T_{1}(\mathbf{X}_{1})$. Hence,

[TABLE]

by simply doing a substitution $t_{11}=1-u_{11}$.

Finally, given there exist $m,M>0$ such that $0<m\leq c(\mathbf{u})\leq M$ for all $\mathbf{u}\in\mathbb{I}^{q}$, and $c_{n}\to c$ uniformly on $\mathbb{I}^{q}$, some basic analysis implies that

[TABLE]

uniformly on $\mathbb{I}^{q}$ as $n\to\infty$, where $c_{ni}$ are the marginal copula densities of $c_{n}$ for $i=1,\dots,k$. Property (A8) in terms of copula densities is then satisfied by the Lebesgue dominated convergence theorem. ∎

Remark 1. For showing (A3), we assumed that $\Phi$ is strictly convex at $1$. While convexity is usually defined as a global property of a function, we use Definition 1 of local strict convexity of Liese and Vajda (2006), i.e. $\Phi$ is strictly convex at $1$ if it is convex and not locally linear at $1$.

Remark 2. Consistency results for dependence measures based on the copula cdf are typically based on the weak uniform convergence of the empirical copula process. In copula density terms, the stated conditions in Proposition 1 for fulfilling (A8) are rather stringent. Indeed, it is known that many of the common copula families (e.g. normal, Student, Clayton, Gumbel) have densities that explode to infinity near some boundaries points, see e.g. Omelka et al. (2009). This means that consistency arguments relying on uniform convergence are typically limited to compact subsets of $\mathbb{I}^{q}$. However, there are theoretical properties that favour copula density based dependence measures, see Remark 3.

Remark 3. When using a dependence measure that compares the true copula cdf to the one under independence, like the Hoeffding’s $\Phi^{2}$ of Geißer et al. (2010) (with $d_{1}=\cdots=d_{k}=1$) using the $L_{2}$-distance, the independence characterization Axiom (A3) is still satisfied, but stays rather ambiguous. The reason is that such dependence measures can be made arbitrarily small, while maintaining an exact deterministic relationship (singularity) between all the variables. An explicit proof of this follows from Theorem 3.2.2 of Nelsen (2006), telling us that we can approximate the independence copula arbitrarily and uniformly closely by copulas (shuffles of Min) exhibiting a perfect deterministic relationship (‘complete dependence’ in the sense of Lancaster (1963)). The copula density based $\Phi$-dependence measures are more alert to such singularities.

It is also interesting to think about the meaning of maximal $\Phi$-dependence. Such maximal dependence occurs if there is a certain singularity (and hence the copula density does not exist everywhere). First, we give an overview of popular choices for $\Phi$, the corresponding name, and its maximum value $\Phi(0)+\Phi^{*}(0)$, see Table 1.

Note that all the $\Phi$-functions in Table 1 are strictly convex on $(0,\infty)$, except for the total variation distance, which is only strictly convex at $1$. The mutual information is a prominent quantity in information theory, see e.g. Cover and Thomas (2006). Differential Shannon entropy quantifies the average amount of uncertainty and mutual information equals the difference between the differential entropy under independence and under the true model. A general family is $\Phi(t)=|t-1|^{\alpha}$ with $\alpha\geq 1$, for which $\Phi(0)=1$ and $\Phi^{*}(0)=\infty$ for $\alpha>1$ (for $\alpha=1$, this is the total variation distance).

We refer to Liese and Vajda (2006) and references therein for further statistical applications, as well as for other choices of $\Phi$. For the Jensen-Shannon divergence, we refer to Österreicher and Vajda (2003).

3.2. Perfect dependence

Next to independence, there is some kind of maximal dependence at the opposite end of the spectrum, to which we will refer as perfect dependence. Perfect dependence is inherent to the dependence measure and occurs if and only if the measure in question reaches its maximum value.

Two random variables $X_{1},X_{2}$ are often seen as maximally dependent if their copula is the Fréchet upper or lower bound copula, that is if $F_{1}(X_{1})=U_{1}=U_{2}=F_{2}(X_{2})$ almost surely, with $F_{i}$ the cdf of $X_{i}$ for $i=1,2$, or if $U_{1}=1-U_{2}$ almost surely. In this case, concordance measures like Kendall’s tau and Spearman’s rho are maximal. This however focuses on monotonic dependence, and perfect co- or counter-monotonicity are only particular cases of ‘strict dependence’ as in Rényi (1959), telling that $X_{2}=\Psi_{1}(X_{1})$ for some function $\Psi_{1}$, or $X_{1}=\Psi_{2}(X_{2})$ for some function $\Psi_{2}$ (deterministic predictability of one variable through the other). When $\Psi_{1}$ (or $\Psi_{2}$) is invertible, we get the ‘complete dependence’ of Lancaster (1963).

More general is the ‘pure dependence’ of Geenens and Lafaye de Micheaux (2022), being the existence of a function $\Psi:\mathbb{R}\rightarrow\mathbb{R}^{2}$ such that $(X_{1},X_{2})=\Psi(U)$ for a certain $U\sim\mathcal{U}[0,1]$, and with $\mathbb{P}_{X_{1}\perp X_{2}}((X_{1},X_{2})\in\Psi(\mathbb{I}))=0$, where $\Psi(\mathbb{I})$ is the image of $\mathbb{I}$ under $\Psi$. We provide an example of pure dependence.

Example 1. Consider $X_{1}\stackrel{{\scriptstyle d}}{{=}}X_{2}$ (equality in distribution) with

[TABLE]

and interconnected by the copula

[TABLE]

where $\mathds{1}$ is the indicator function. Figure 1 shows scatter plots of $(U_{1},U_{2})$ and $(X_{1},X_{2})$ based on a random sample of size $200$. This is clearly an example of pure dependence, with function $\Psi$ given by $\Psi(U)=$ $(\cos(2\pi U),\sin(2\pi U))$. Note however that this is not an example of strict dependence nor of complete dependence.

We can easily extent this notion of pure dependence to $q$ random variables $X_{1},\dots,X_{q}$, and formulate it in terms of $F_{1}(X_{1})=U_{1},\dots,F_{q}(X_{q})=U_{q}$, as the existence of a $\Psi:\mathbb{I}^{q-1}\rightarrow\mathbb{I}^{q}$ such that $(U_{1},\dots,U_{q})=\Psi(V_{1},\dots,V_{q-1})$ for certain $V_{1},\dots,V_{q-1}\sim\mathcal{U}[0,1]$, and with $\lambda^{q}(\Psi(\mathbb{I}^{q-1}))=0$. The intuition associated with pure dependence is akin to understanding perfect dependence inherent in the $\Phi$-dependence measures.

Theorem 1 of De Keyser and Gijbels (2023) suggests that perfect $\Phi$-dependence ($\mathcal{D}_{\Phi}$ maximal) between random vectors $\mathbf{X}_{1},\dots,\mathbf{X}_{k}$ arises when there exists a $B\in\mathcal{B}(\mathbb{I}^{q})$ such that $\mu_{C}(B)=1$ and $(\mu_{C_{1}}\times\cdots\times\mu_{C_{k}})(B)=0$ (i.e. $\mu_{C}\perp\mu_{C_{1}}\times\cdots\times\mu_{C_{k}}$, hence $\mu_{C}=\mu_{C}^{\text{s}})$, where $\mathcal{B}(\mathbb{I}^{q})$ is the Borel sigma-algebra on $\mathbb{I}^{q}$. Restricting to $\Phi$-functions that are strictly convex at $1$, and satisfy $\Phi(0)+\Phi^{*}(0)<\infty$, this is the only case in which $\mathcal{D}_{\Phi}$ is maximal (and hence we obtain a characterization of perfect dependence). We illustrate this singularity between random vectors by means of another example.

Example 2. Consider $(X_{1},X_{2},X_{3},X_{4})$ having copula

[TABLE]

for $\theta\in[1,\infty)$. Then $(X_{1},X_{2})$ and $(X_{3},X_{4})$ have a Gumbel$(\theta)$ copula and $(X_{1},X_{3})$ and $(X_{2},X_{4})$ the comonotonicity copula. We have $\mu_{C}\perp\lambda^{4}$, since $\mu_{C}$ is concentrated on $\{(u_{1},u_{2},u_{3},u_{4})\in\mathbb{I}^{4}:u_{1}=u_{3},u_{2}=u_{4}\}$. Also, \mathcal{D}_{\Phi}\big{(}(X_{1},X_{2});(X_{3},X_{4})\big{)}=\infty, since $\mu_{C}\perp\mu_{C_{1}}\times\mu_{C_{2}}$ ($C_{1}$ the copula of $(X_{1},X_{2})$ and $C_{2}$ the copula of $(X_{3},X_{4}))$.

On the other hand, we have \mathcal{D}_{\Phi}\big{(}(X_{1},X_{3});(X_{2},X_{4})\big{)}<\infty since $\mu_{C_{1}}\perp\lambda^{2}$ and $\mu_{C_{2}}\perp\lambda^{2}$ (now defining $C_{1}$ as the copula of $(X_{1},X_{3})$ and $C_{2}$ as the copula of $(X_{2},X_{4})$). In fact, this is an example of $\mu_{C}\perp\lambda^{4}$, $\mu_{C_{1}}\times\mu_{C_{2}}\perp\lambda^{4}$, but $\mu_{C}\ll\mu_{C_{1}}\times\mu_{C_{2}}$. In Example 6, we illustrate that \mathcal{D}_{\Phi}\big{(}(X_{1},X_{3});(X_{2},X_{4})\big{)} (for one particular choice of $\Phi$) is actually nothing more than $\mathcal{D}_{\Phi}$ of a two dimensional Gumbel$(\theta)$ copula, which is pretty intuitive, since $(X_{1},X_{3})$ and $(X_{2},X_{4})$ both have the comonotonicity copula and thus $X_{1}$ and $X_{3}$ can be seen as one, as well as $X_{2}$ and $X_{4}$, and $(X_{1},X_{2})$ (just as $(X_{3},X_{4})$) has a Gumbel$(\theta)$ copula.

In short, $\Phi$-dependence is maximal when there is a singularity among random variables belonging to different random vectors, indifferent to plausible singularities between random variables within one and the same random vector.

4. A Gaussian copula approach

In this section, we assume a Gaussian copula model for $\mathbf{X}=(\mathbf{X}_{1},\dots,\mathbf{X}_{k})$. Although a restricted framework, it leads to interesting statistical inference results. In fact, we only assume finite second moments and the existence of the covariance matrix of $\mathbf{X}$, say

[TABLE]

where $\mathbf{\Sigma}_{ij}\in\mathbb{R}^{d_{i}\times d_{j}}$ contains the covariances between the components of $\mathbf{X}_{i}$ and $\mathbf{X}_{j}$ for $i\neq j$, and $\mathbf{\Sigma}_{ii}\in\mathbb{R}^{d_{i}\times d_{i}}$ the within covariances of $\mathbf{X}_{i}$ for $i=1,\dots,k$. If $\mathbf{X}$ has a normal $\mathcal{N}(\mathbf{0}_{q},\mathbf{\Sigma})$ distribution, a straightforward calculation shows that (1) reduces to

[TABLE]

with

[TABLE]

the covariance matrix under mutual independence of $\mathbf{X}_{1},\dots,\mathbf{X}_{k}$, and $\mathbf{B}^{-1}=\mathbf{\Sigma}^{-1}-\mathbf{\Sigma}_{0}^{-1}$. In (4), we used a superscript $\mathcal{N}$ referring to the Gaussian assumption and emphasize the dependence on merely the covariance matrix $\mathbf{\Sigma}$. For certain specific choices of the function $\Phi$, the above integral allows a closed form solution in terms of the covariance matrix $\mathbf{\Sigma}$. We shall go deeper into the cases $\Phi(t)=t\log(t)$ and $\Phi(t)=(\sqrt{t}-1)^{2}$, and denote the resulting dependence measures as $\mathcal{D}_{t\log(t)}^{\mathcal{N}}$ and $\mathcal{D}_{(\sqrt{t}-1)^{2}}^{\mathcal{N}}$ (for the other choices of $\Phi$ listed in Table 1, there is no such elegant closed form solution to (4)). It is a straightforward calculation to show that

[TABLE]

and, denoting $\mathbb{I}_{q}$ for the $q\times q$ identity matrix,

[TABLE]

It is also easily checked that, if $k=2$ and $d_{1}=d_{2}=2$, i.e. in the case of two univariate random variables, the expression (6) reduces to what was found in Section 4 of Geenens and Lafaye de Micheaux (2022).

Notice that (5) and (6) are formulated in terms of covariance matrices. However, since our dependence measures are copula-based, they should be scale-invariant. Let $\mathbf{D}_{\mathbf{\Sigma}}$ be the diagonal matrix with the variances on the diagonal. Then, the correlation matrix $\mathbf{R}$ corresponding to $\mathbf{\Sigma}$ satisfies $\mathbf{R}=\mathbf{D}_{\mathbf{\Sigma}}^{-1/2}\mathbf{\Sigma}\mathbf{D}_{\mathbf{\Sigma}}^{-1/2}$, such that $|\mathbf{\Sigma}|$ is just a rescaling of $|\mathbf{R}|$ by multiplying the latter with the product of all variances. From this, we easily see that $\mathcal{D}_{t\log(t)}^{\mathcal{N}}(\mathbf{\Sigma})=\mathcal{D}_{t\log(t)}^{\mathcal{N}}(\mathbf{R})$. Moreover,

[TABLE]

with $\mathbf{R}_{0}$ the correlation matrix corresponding to $\mathbf{\Sigma}_{0}$, implying that $\mathcal{D}_{(\sqrt{t}-1)^{2}}^{\mathcal{N}}(\mathbf{\Sigma})=\mathcal{D}_{(\sqrt{t}-1)^{2}}^{\mathcal{N}}(\mathbf{R})$. For general $\Phi$, we obtain this relation by performing the substitution $\mathbf{y}=\mathbf{D}_{\bm{\Sigma}}^{-1/2}\mathbf{x}$ in (4).

Also, since we do not care about univariate marginals, $\mathbf{R}$ needs not to be the traditional Pearson correlation matrix of a multivariate normal distribution, but can be the margin free Gaussian copula correlation matrix of normal scores, that is

[TABLE]

for $i,m=1,\dots,k$, $j=1,\dots,d_{i}$ and $t=1,\dots,d_{m}$, where Corr stands for the traditional Pearson correlation, and $\phi^{-1}$ for the standard normal quantile function.

Statistical inference

Based on a sample $\mathbf{X}^{(\ell)}=(\mathbf{X}_{1}^{(\ell)},\dots,\mathbf{X}_{k}^{(\ell)})$ for $\ell=1,\dots,n$ from $\mathbf{X}$, with $\mathbf{X}_{i}^{(\ell)}=(X_{i1}^{(\ell)},\dots,X_{id_{i}}^{(\ell)})$ for $\ell=1,\dots,n$ a sample from $\mathbf{X}_{i}$ for $i=1,\dots,k$, the sample version of (7) is known as the matrix of normal scores rank correlation coefficients (Hájek and Šidák (1967)),

[TABLE]

defined through normal scores

[TABLE]

obtained from the univariate empirical cdf $\widehat{F}_{ij}(x_{ij})=\frac{1}{n}\sum_{\ell=1}^{n}\mathds{1}\{X_{ij}^{(\ell)}\leq x_{ij}\}$ for $i=1,\dots,k$ and $j=1,\dots,d_{i}$. The quantity $\widehat{\rho}_{ij,mt}$ is computed as the conventional Pearson correlation of the bivariate sample of scores \big{(}(\widehat{Z}_{ij}^{(1)},\widehat{Z}_{mt}^{(1)}),\dots,(\widehat{Z}_{ij}^{(n)},\widehat{Z}_{mt}^{(n)})\big{)} and by noting that

[TABLE]

which follows from the fact that $\phi^{-1}(\alpha)=-\phi^{-1}(1-\alpha)$ for $\alpha\in[0,1]$ and $n\widehat{F}_{ij}(X_{ij}^{(\ell)})$ is the rank of $X_{ij}^{(\ell)}$ in the sample $X_{ij}^{(1)},\dots,X_{ij}^{(n)}$. Being rank-based, makes the variance of the normal scores independent of the data.

A next natural step in estimating $\mathcal{D}_{\Phi}^{\mathcal{N}}(\mathbf{R})$ is to just plug in $\widehat{\mathbf{R}}_{n}$ instead of the unknown matrix $\mathbf{R}$. Let $\mathbb{S}^{q}$ be the set of all $q\times q$ covariance matrices and $\mathbb{S}^{q}_{>}\subset\mathbb{S}^{q}$ the set of all positive definite ones. Let $\varphi$ be the map defined by $\varphi(\mathbf{\Sigma})=\mathbf{D}_{\mathbf{\Sigma}}^{-1/2}\mathbf{\Sigma}\mathbf{D}_{\mathbf{\Sigma}}^{-1/2}$ for $\mathbf{\Sigma}\in\mathbb{S}^{q}$, and $||\cdot||_{\text{F}}$ the Frobenius matrix norm, i.e. $||\mathbf{\Sigma}||_{\text{F}}^{2}=\text{Tr}(\mathbf{\Sigma}^{\text{T}}\mathbf{\Sigma})$. If we can show the Fréchet differentiability of the mapping

[TABLE]

on $\mathbb{S}^{q}_{>}$, then the delta method turns an asymptotic normality result for $\widehat{\mathbf{R}}_{n}$ into an asymptotic normality result for $\mathcal{D}_{\Phi}^{\mathcal{N}}(\widehat{\mathbf{R}}_{n})$. For general $\Phi$, interchanging Fréchet differentiation and Lebesgue integration in (4) would be useful, and holds if there exists an integrable function on $\mathbb{R}^{q}$ dominating the Fréchet derivative of the integrand of (4) uniformly on $\mathbb{S}^{q}_{>}$. Here, special attention is again devoted to $\Phi(t)$ either $t\log(t)$ or $(\sqrt{t}-1)^{2}$, because we have the more explicit expressions (5) and (6). Theorem 1 states formally the asymptotic normality result for $\mathcal{D}_{\Phi}^{\mathcal{N}}(\widehat{\mathbf{R}}_{n})$.

Theorem 1. Let $\mathbf{X}$ have a Gaussian copula with correlation matrix $\mathbf{R}\in\mathbb{S}^{q}_{>}$, and let $\widehat{\mathbf{R}}_{n}$ be given by (8), based on which the plug-in estimator $\widehat{\mathcal{D}}^{\mathcal{N}}_{\Phi,n}=\mathcal{D}_{\Phi}^{\mathcal{N}}(\widehat{\mathbf{R}}_{n})$ is constructed. If the mapping defined in (9) is Fréchet differentiable, then the estimator $\widehat{\mathcal{D}}_{\Phi,n}^{\mathcal{N}}$ is asymptotically normal. Moreover, if differentiation and integration in (4) can be interchanged, we have in addition the expression for the asymptotic variance:

[TABLE]

with asymptotic variance

[TABLE]

where $\mathbf{D}_{\mathbf{M}_{\Phi}\mathbf{R}}$ is the diagonal matrix consisting of the diagonal of $\mathbf{M}_{\Phi}\mathbf{R}$, and with

[TABLE]

defining $\alpha(\mathbf{X})=\Phi(k(\mathbf{X}))$, $\alpha^{\prime}(\mathbf{X})=\Phi^{\prime}(k(\mathbf{X}))$, and

[TABLE]

for

[TABLE]

In particular, one can show that

[TABLE]

and

[TABLE]

defining

[TABLE]

The expression for $\mathbf{M}_{t\log(t)}$ in Theorem 1 was obtained in the proof using the explicit formula (5) that we have for the mutual information (similarly for the Hellinger distance). In the following example, we verify, for $\Phi(t)=t\log(t)$, the general formula for $\mathbf{M}_{\Phi}$ relying on interchanging the Fréchet derivative and Lebesgue integral.

Example 3. In case $\Phi(t)=t\log(t)$, the different terms in the general formula for $\text{Tr}(\mathbf{M}_{\Phi}\mathbf{H})$, for a certain $\mathbf{H}\in\mathbb{S}^{q}_{>}$, can be calculated as follows. First, using (A1) in the proof of Theorem 1, we have

[TABLE]

where

[TABLE]

Next, observe that

[TABLE]

Using $\Phi^{\prime}(t)=\log(t)+1$, it is also straightforward to see that

[TABLE]

and

[TABLE]

and finally

[TABLE]

Combining all this, we obtain that

[TABLE]

[TABLE]

where we used the following properties of quadratic forms

[TABLE]

for certain compatible matrices $\mathbf{A},\mathbf{A}_{1},\mathbf{A}_{2}$, the cyclic trace property, and the fact that $\text{Tr}((\mathbf{R}^{-1}-\mathbf{R}_{0}^{-1})\mathbf{R})=0$ and $\text{Tr}\left(\mathbf{D}\mathbf{R}\right)=-\text{Tr}(\mathbf{R}_{0}^{-1}\mathbf{H})$. Hence, we have verified that the general formula for $\mathbf{M}_{\Phi}$ also brings us to

[TABLE]

We end this section by looking at a specific four dimensional Gaussian copula family.

Example 4. Consider a four dimensional random vector $(X_{1},X_{2},X_{3},X_{4})$ having a Gaussian copula with correlation matrix

[TABLE]

Then one can check that

[TABLE]

and

[TABLE]

In (11), we normalized the Hellinger distance, guaranteeing a dependence measurement in $\mathbb{I}$ (recall Proposition 1). For the mutual information, an artificial normalization is required, but we do not do not implement it here. Figure 2 shows how (10) and (11) depend on $\rho_{2}$ for different values of $\rho_{1}$. Some observations are:

• •

$\mathcal{D}^{\mathcal{N}}_{t\log(t)}=\mathcal{D}^{\mathcal{N}}_{(\sqrt{t}-1)^{2}}=0$ iff $\rho_{2}=0$.

• •

For $\rho_{1}\to-1$ (singularity of $\mu_{C_{1}}$ and $\mu_{C_{2}}$ w.r.t. $\lambda^{2}$), we must have $\rho_{2}\to 0$ and see that $\mathcal{D}^{\mathcal{N}}_{t\log(t)}\to 0$ and $\mathcal{D}^{\mathcal{N}}_{(\sqrt{t}-1)^{2}}\to 0$.

• •

If $\rho_{1}=2|\rho_{2}|-1$ (absolute continuity of $\mu_{C_{1}}$ and $\mu_{C_{2}}$ w.r.t. $\lambda^{2}$, but singularity of $\mu_{C}$ w.r.t. $\lambda^{4}$), we get $\mathcal{D}^{\mathcal{N}}_{t\log(t)}=\infty$ and $\mathcal{D}^{\mathcal{N}}_{(\sqrt{t}-1)^{2}}=1$ .

• •

For $\rho_{1}\to 1$ (singularity of $\mu_{C_{1}}$ and $\mu_{C_{2}}$ w.r.t. $\lambda^{2}$), $\mathcal{D}^{\mathcal{N}}_{t\log(t)}\to-1/2\log(1-\rho_{2}^{2})$, and $\mathcal{D}^{\mathcal{N}}_{(\sqrt{t}-1)^{2}}/2\to 1-(2(1-\rho_{2}^{2})^{1/4})/(4-\rho_{2}^{2})^{1/2}$, being the mutual information and half Hellinger distance of a bivariate Gaussian copula with correlation $\rho_{2}$ and is maximal iff $|\rho_{2}|=1$.

Note that the principal components of $(Z_{1},Z_{2})=((\phi^{-1}\circ F_{1})(X_{1}),(\phi^{-1}\circ F_{2})(X_{2}))$ are

[TABLE]

and similarly, those of $(Z_{3},Z_{4})=((\phi^{-1}\circ F_{3})(X_{3}),(\phi^{-1}\circ F_{4})(X_{4}))$ are

[TABLE]

Moreover, $\text{Corr}(PC_{1},PC_{3})=\text{Corr}(PC_{2},PC_{3})=\text{Corr}(PC_{1},PC_{4})=0$, but

[TABLE]

And so, we see that if $\rho_{1}=2|\rho_{2}|-1$, $|\text{Corr}(PC_{2},PC_{4})|=1$, i.e. the principal components $PC_{2}$ and $PC_{4}$ are perfectly correlated. This means that the four dimensional random vector $(X_{1},X_{2},X_{3},X_{4})$ is propagated in a three dimensional subspace (scatterplot of $(Z_{1},Z_{2},Z_{3},Z_{4})$ constitutes a hyperplane), resulting in the singularity of $\mu_{C}$ with respect to $\mu_{C_{1}}\times\mu_{C_{2}}$ and explaining maximal dependence.

Figure 3 shows the asymptotic standard deviation $\zeta_{t\log(t)}$ of the mutual information and $\zeta_{(\sqrt{t}-1)^{2}}/2$ of the half Hellinger distance (as in Theorem 1) as a function of $\rho_{2}$ for different values of $\rho_{1}$, which can be calculated as

[TABLE]

In general, we see that higher degrees of dependence come with higher asymptotic variance. For the Hellinger distance, however, the asymptotic variance goes down to zero when we get close to the previously discussed singularity, i.e. when $\rho_{2}$ gets close to satisfying $|\rho_{2}|=(\rho_{1}+1)/2$. For instance, if $\rho_{1}=0$, the asymptotic variance is maximal at $\left|\rho_{2}\right|\approx 0.45427$, after which it converges to zero for $\left|\rho_{2}\right|\to 0.5$. From its mathematical expression above, and its general form in Theorem 1, we see the factor $\exp((-1/2)\mathcal{D}_{t\log(t)}^{\mathcal{N}})\to 0$ in case of singularity makes the asymptotic variance of the Hellinger distance tend to zero.

5. A maximum likelihood approach

A Gaussian copula model is restricted to monotone dependence structures in terms of correlations. However, in many cases other relationships are present. Think for instance of comovements in the tails of stock returns. Thanks to Proposition 1 (in particular the fulfillment of Axiom (A3)), we know that $\Phi$-dependence measures are able to capture such associations as well. In practice, a sufficiently flexible estimation methodology is required. If one is willing to, for example, assume a Clayton copula model for two random variables $X_{1}$ and $X_{2}$, we can use this assumption to estimate $\mathcal{D}_{\Phi}$ and hence, next to a monotone relationship, possible lower tail dependence is guaranteed to be incorporated as well.

5.1. Maximum likelihood-based inference

In general, in this section, we assume a specified parametric model for the copula density of $\mathbf{X}=(\mathbf{X}_{1},\dots,\mathbf{X}_{k})$, say $\mathcal{P}_{C}=\{c(\hskip 2.84544pt\cdot\hskip 2.84544pt;\bm{\theta}_{C}):\bm{\theta}_{C}\in\bm{\Theta}_{C}\subset\mathbb{R}^{D}\}$, where $D=\text{dim}(\bm{\Theta}_{C})$. Notice that a Gaussian copula is also a parametric copula of this form when we stack all upper triangle elements of the correlation matrix in a vector. The univariate marginals $F_{ij}$ for $i=1,\dots,k$ and $j=1,\dots,d_{i}$ are approached in either a parametric, or non-parametric way, leading to two different estimation procedures.

Case 1: parametric marginals. A first option is to also assume a parametric statistical density model for $F_{ij}$, say $\mathcal{P}_{ij}=\{f_{ij}(\hskip 2.84544pt\cdot\hskip 2.84544pt;\bm{\theta}_{ij}):\bm{\theta}_{ij}\in\mathbf{\Theta}_{ij}\}$. Based on a sample $\mathbf{X}^{(1)},\dots,\mathbf{X}^{(n)}$ from $\mathbf{X}$, the full log-likelihood is

[TABLE]

Let $\bm{\eta}=(\bm{\theta}_{11},\dots,\bm{\theta}_{kd_{k}},\bm{\theta}_{C})$ and $\widehat{\bm{\eta}}_{n}^{\text{MLE}}=(\widehat{\bm{\theta}}_{11,n}^{\text{MLE}},\dots,\widehat{\bm{\theta}}_{kd_{k},n}^{\text{MLE}},\widehat{\bm{\theta}}_{C,n}^{\text{MLE}})$ be the MLE of $\bm{\eta}$, obtained by maximizing (12). It is a very well-known result that under certain regularity conditions (true $\bm{\eta}$ lies in an open subset in which the log-density admits third derivatives w.r.t. the parameters that are uniformly bounded by an integrable function, and interchanging of derivatives and integral for score equations, see e.g. Theorem 4.1 on page 429 in Lehmann (1983)), one has that $\widehat{\bm{\eta}}_{n}^{\text{MLE}}$ is an asymptotically normal estimator. If $\bm{\theta}_{C}$ is $m$-dimensional and $\bm{\theta}_{ij}$ all one dimensional, (12) is a $q+m$-dimensional optimization problem. In case of a multivariate normal distribution, one has for $\bm{\theta}_{C}$ the vector of upper triangle elements of the correlation matrix of the normal copula, and $\bm{\theta}_{ij}=(\mu_{ij},\sigma_{ij}^{2})$ the mean and variance of $X_{ij}$.

Given estimators $\widehat{\bm{\theta}}_{ij,n}^{\text{IFM}}$ of $\bm{\theta}_{ij}$ (e.g. also based on MLE using the sample $X_{ij}^{(1)},\dots,X_{ij}^{(n)}$), the pseudo log-likelihood for estimating $\bm{\theta}_{C}$ is

[TABLE]

If $\bm{\theta}_{C}$ is $m$-dimensional and $\bm{\theta}_{ij}$ all one dimensional, this is an $m$-dimensional optimization problem, after having done $q$ one dimensional optimization problems. This method is known as the inference functions for margins (IFM) method, and extensively studied in Section 10.1 of Joe (1997). Asymptotic normality is known to hold under the same regularity conditions as for the MLE. For a multivariate normal distribution, it is known that $\widehat{\bm{\eta}}_{n}^{\text{IFM}}=\widehat{\bm{\eta}}_{n}^{\text{MLE}}$.

Case 2: non-parametric marginals. If no appropriate parametric models can be proposed for the marginals, an option consists of estimating the marginals via the empirical cdf

[TABLE]

and the copula parameter $\bm{\theta}_{C}$ through maximizing

[TABLE]

resulting in an estimator $\widehat{\bm{\theta}}_{C,n}^{\text{NP}}$ for $\bm{\theta}_{C}$. Again under the same regularity conditions, it is shown in Genest et al. (1995) that $\widehat{\bm{\theta}}^{\text{NP}}_{C,n}$ is asymptotically normal. If $c$ is the Gaussian copula density, there is no known expression for the maximizer of (13) over all correlation matrices, and numerical maximization is often used (see Section 5.5.3 in McNeil et al. (2005) for more details), being quite unfeasible in high dimensions. Of course, we have the results from Section 4 dealing with the Gaussian copula case.

In all the above cases, we have an asymptotically normal estimator $\widehat{\bm{\theta}}_{n}=\widehat{\bm{\theta}}_{C,n}$ for $\bm{\theta}_{C}$. A natural estimator for the copula density $c(\hskip 2.84544pt\cdot\hskip 2.84544pt;\bm{\theta}_{C})$ is then $\widehat{c}(\hskip 2.84544pt\cdot\hskip 2.84544pt)=c(\hskip 2.84544pt\cdot\hskip 2.84544pt;\widehat{\bm{\theta}}_{n})$. If $c$ is unbounded, we will generally not have uniform consistency of $\widehat{c}$ (recall Remark 2). Nevertheless, Fréchet differentiability of the mapping

[TABLE]

suffices to turn the asymptotic normality result of $\widehat{\bm{\theta}}_{n}$ into an asymptotic normality result for the plug-in estimator $\mathcal{D}_{\Phi}(\widehat{c})=\mathcal{D}_{\Phi}(c(\hskip 2.84544pt\cdot\hskip 2.84544pt;\hskip 2.84544pt\widehat{\bm{\theta}}_{n}))=\mathcal{D}_{\Phi}(\widehat{\bm{\theta}}_{n})$. Still, the estimator $\mathcal{D}_{\Phi}(\widehat{c})$ would often require high-dimensional numerical integration, because usually the integral in (14) does not have a closed-form expression in terms of $\bm{\theta}_{C}$. That is why we suggest the following general approach.

Let $M_{n}$ be a user defined sample size for each $n$, and $\widetilde{\mathbf{U}}^{(1)},\dots,\widetilde{\mathbf{U}}^{(M_{n})}$ with $\widetilde{\mathbf{U}}^{(\ell)}=(\widetilde{\mathbf{U}}_{1}^{(\ell)},\dots,\widetilde{\mathbf{U}}_{k}^{(\ell)})$ for $\ell=1,\dots,M_{n}$, a sample from $\widetilde{\mathbf{U}}$ having distribution $C(\hskip 2.84544pt\cdot\hskip 2.84544pt;\widetilde{\bm{\theta}})$ given that $\widehat{\bm{\theta}}_{n}=\widetilde{\bm{\theta}}$, that is

[TABLE]

We then propose

[TABLE]

as estimator for the $\Phi$-dependence $\mathcal{D}_{\Phi}$ given in (1). The rationale is that $\Phi$-dependence measures can be seen as an expectation and the empirical mean in (15) tries to approximate this integral. Since we have an explicit form of the estimated copula, we can take as many samples as we want when $\widehat{\bm{\theta}}_{n}$ is given, and the larger $M_{n}$, the better the approximation will be. It is intuitively clear (yet not trivial to prove) that, for $M_{n}$ large enough in some sense, the asymptotic properties of $\mathcal{D}_{\Phi}(\widehat{c})$ carry over to $\widehat{\mathcal{D}}_{\Phi,n,M_{n}}$. A conditional law of large numbers will give substance to this, see Remark 4. Theorem 2 states the asymptotic normality result for $\widehat{\mathcal{D}}_{\Phi,n,M_{n}}$.

Theorem 2. Let $\widehat{\bm{\theta}}_{n}$ be an asymptotically normal estimator for $\bm{\theta}_{C}$ based on which the estimator (15) is constructed, where the user chosen parameter $M_{n}$ is such that

[TABLE]

for $n\to\infty$. If the mapping defined in (14) is Fréchet differentiable, the estimator $\widehat{\mathcal{D}}_{\Phi,n,M_{n}}$ is asymptotically normal. Moreover, if the derivative can be moved into the integral, we have the expression for the asymptotic variance-covariance:

[TABLE]

as $n\to\infty$, where $\mathbf{V}$ is the asymptotic variance-covariance matrix of $\widehat{\bm{\theta}}_{n}$, and

[TABLE]

defining $D=\text{dim}(\bm{\Theta}_{C})$, the $i$’th component of $\bm{\theta}_{C}$ as $\theta_{C,i}$, and

[TABLE]

Remark 4. A conditional version of Kolmogorov’s law of large numbers (see Theorem 4.2 in Majerek et al. (2005)) yields, for all $n\in\mathbb{N}$,

[TABLE]

as $M_{n}\to\infty$. Formally, this means that

[TABLE]

This implies that the condition on $M_{n}$ in Theorem 2, stating that

[TABLE]

is reasonable. Indeed, take $\epsilon,\delta>0$ arbitrary. Take for example $N=1$. Let $n\geq 1$ arbitrary. For these $\epsilon,\delta$ and $n$, we can take $M_{n}=M_{n}^{(0)}$ since $M_{n}$ is user defined, and hence

[TABLE]

A small value of $N$ corresponds to a good approximation of the integral in $\mathcal{D}_{\Phi}$ already for smaller $n$.

Before looking into another class of parametric copulas, we focus on the Gaussian setting once more.

Example 5. Consider the simple case of estimating, for two univariate random variables $X_{1},X_{2}$ with cdf’s $F_{1},F_{2}$ and Gaussian copula, the parameter

[TABLE]

with $\phi^{-1}$ the standard normal quantile function and Corr the Pearson correlation, and afterwards their mutual information $\mathcal{D}_{t\log(t)}(\rho)=-(1/2)\log(1-\rho^{2})$ via $\mathcal{D}_{t\log(t)}(\widehat{\rho})$ for a certain estimator $\widehat{\rho}$ of $\rho$. The following semi-parametric approaches might be considered.

• •

Approach 1: Assume a Gaussian copula model and estimate the marginals non-parametrically. In particular, take the estimator (8).

• •

Approach 2: Assume a Gaussian copula model, estimate the marginals non-parametrically, and the copula parameter via the pseudo-likelihood (13).

If we look at Approach 1 from a matrix point of view, i.e. estimator $\widehat{\mathbf{R}}_{n}$ in (8), we can apply Theorem 1 with

[TABLE]

such that

[TABLE]

as $n\to\infty$. From a one parameter $\rho$ point of view, with estimator $\widehat{\rho}_{1}$ that is on the off-diagonal of $\widehat{\mathbf{R}}_{n}$, we know that

[TABLE]

as $n\to\infty$, since, as seen in the proof of Theorem 1, the estimator $\widehat{\rho}_{1}$ has the same asymptotic distribution as in the case where the marginals are known, i.e. as the usual sample Pearson correlation in case of a bivariate normal distribution. Hence, noting that $d/d\rho\mathcal{D}_{t\log(t)}(\rho)=\rho/(1-\rho^{2})$, the univariate delta method implies that

[TABLE]

as it should.

As for Approach 2, the resulting estimator, say $\widehat{\rho}_{2}$, has an asymptotic variance that is rather complex to calculate. In Genest et al. (1995), it is shown that this asymptotic variance cannot be smaller than the one of the maximum likelihood estimator in case the marginals are known. Thus, since $\widehat{\rho}_{1}$ has the same asymptotic distribution as the maximum likelihood estimator in a bivariate Gaussian model, $\widehat{\rho}_{2}$ cannot have a smaller asymptotic variance than $\widehat{\rho}_{1}$, and hence the asymptotic variance of $\mathcal{D}_{t\log(t)}(\widehat{\rho}_{2})$ cannot be smaller than the asymptotic variance of $\mathcal{D}_{t\log(t)}(\widehat{\rho}_{1})$. We conclude that Approach 1, in which we have a nice explicit formula for the estimator, always performs at least equally well as Approach 2 in terms of asymptotic variance.

5.2. Nested Archimedean copulas framework

We now turn some attention to a specific parametric family of copulas. The hierarchical models of nested Archimedean copulas extend the frequently used Archimedean copulas and are an instinctive choice for modelling dependence between random vectors. See e.g. Hofert and Pham (2013) for the following definition.

Definition 2. (nested Archimedean copulas) A nested Archimedean copula $C$ with two nesting levels and $k$ child copulas is given by

[TABLE]

where $k$ denotes the dimension of the root copula $C_{0}$, and each child copula $C_{i}$ for $i\in\{0,\dots,k\}$ is an Archimedean copula with a completely monotone generator $\psi_{i}$, that is

[TABLE]

and $\psi_{i}:[0,\infty)\to[0,1]$ is continuous with $\psi_{i}(0)=1,\lim_{t\to\infty}\psi_{i}(t)=0$ and $(-1)^{\ell}\psi_{i}^{(\ell)}(t)\geq 0$ for all $\ell\in\mathbb{N},t\in(0,\infty)$.

Note that we can further nest the child copulas in (16), although this is superfluous for our purposes (in particular because that makes densities excessively complicated). The condition $(-1)^{\ell}\psi_{i}^{(\ell)}(t)\geq 0$ for all $\ell\in\mathbb{N},t\in(0,\infty)$ is called complete monotonicity of the function $\psi_{i}$ and a sufficient condition to guarantee that (16) indeed is a copula, is that $\psi_{0}^{-1}\circ\psi_{i}$ for all $i\in\{1,\dots,k\}$ have completely monotone first order derivatives. The latter condition is often softened (but definitely not equivalent) to the sufficient nesting condition (e.g. Okhrin and Ristig (2014)), telling us that $\psi_{i}$ all being in a same family of Archimedean copulas for $i\in\{0,\dots,k\}$, say with parameter $\theta_{i}$, such that $\theta_{0}\leq\theta_{i}$ for $i\in\{1,\dots,k\}$ suffices to have complete monotonicity of these derivatives. Hofert and Pham (2013) explicitly calculated the copula density of (16) in several settings. Before looking at an example of the behaviour of a $\Phi$-dependence measure in case of a nested Archimedean copula, we make the following remark.

Remark 5. The estimator (15) is quite a general one in terms of a certain copula family and function $\Phi$, motivated from the population version (1). In some cases, it might be better to first simplify (1) and afterwards do the empirical mean approximation. For example, suppose the interest is in the Hellinger distance $\Phi(t)=(\sqrt{t}-1)^{2}$. A straightforward calculation shows that (1) can also be written as

[TABLE]

suggesting

[TABLE]

where $(\mathbf{U}_{1}^{(\ell)},\dots,\mathbf{U}_{k}^{(\ell)})$ for $\ell=1,\dots,M$ is a sample drawn from $c(\hskip 0.28436pt\cdot\hskip 0.28436pt;\bm{\theta}_{C})$, as numerical approximation, assuming here that no estimation is done (if $\bm{\theta}_{C}$ is estimated, we use the notation $\widetilde{\mathcal{D}}_{(\sqrt{t}-1)^{2},n,M_{n}}$). The benefit of $\widetilde{\mathcal{D}}_{(\sqrt{t}-1)^{2}}$ is that

[TABLE]

guaranteeing fast convergence of the law of large numbers, while the variance of the summand in (15) might be infinite. We further illustrate this in the following example.

Example 6. Consider $(X_{1},X_{2},X_{3},X_{4})$ having a four dimensional partially nested Archimedean copula given by

[TABLE]

where $\psi_{i}(t)=\exp(-t^{1/\theta_{i}})$ is the generator of a Gumbel copula with parameter $\theta_{i}\in[1,\infty)$ for $i=0,1,2$ satisfying $\theta_{0}\leq\theta_{1}$ and $\theta_{0}\leq\theta_{2}$ (sufficient nesting condition). We numerically approximate half the Hellinger distance $\mathcal{D}_{(\sqrt{t}-1)^{2}}/2$ using (17) with $M=10\hskip 1.42271pt000$. Note that no estimation of the marginals is involved here. Figure 4 shows half the Hellinger distance as a function of $\theta_{1}=\theta_{2}$ for different values of $\theta_{0}$, and as a function of $\theta_{0}$ for different values of $\theta_{1}=\theta_{2}$.

Note that if $\theta_{0}=1$, we have $\psi_{0}(t)=\exp(-t)$ such that

[TABLE]

meaning that $(X_{1},X_{2})$ and $(X_{3},X_{4})$ are independent. In general, we observe that the strength of dependence between $(X_{1},X_{2})$ and $(X_{3},X_{4})$ is predominantly determined by the parameter $\theta_{0}$. Nested Archimedean copulas allow us to on the one hand control the dependence within each random vector (by parameters $\theta_{1}$ and $\theta_{2}$ here) and on the other hand control what remains to be specified between the random vectors (parameter $\theta_{0}$ here). Notice that if $\theta_{1}\to\infty$ and $\theta_{2}\to\infty$ (i.e. $(X_{1},X_{2})$ and $(X_{3},X_{4})$ both tend to have a comonotonicity copula), the half Hellinger distance $\mathcal{D}_{(\sqrt{t}-1)^{2}}/2$ tends to the half Hellinger distance of a bivariate Gumbel copula with parameter $\theta_{0}$ (red line in the right plot).

Coming back to Remark 5, suppose we would compute the Hellinger distance of this bivariate Gumbel copula with parameter $\theta_{0}$, say $c(u_{1},u_{2};\theta_{0})$, using the non-simplified empirical version (15) (but now with known $\theta_{0}$ and a true sample from $c$ with fixed $M$) instead of $\eqref{eq: helapprox}$, i.e. using

[TABLE]

where $(U_{1}^{(\ell)},U_{2}^{(\ell)})$ for $\ell=1,\dots,M$ is a sample drawn from $c$. Then, the pitfall is that

[TABLE]

[TABLE]

for $\theta_{0}\geq 2$ since

[TABLE]

for a fixed $u_{1}\in(0,1)$. The above can be seen from the fact that for the Gumbel generator $\psi(t)=\exp(-t^{1/\theta_{0}})$, it holds that $\psi^{\prime}(\psi^{-1}(u_{2}))=\mathcal{O}((-\log(u_{2}))^{1-\theta_{0}})$ as $u_{2}\to 1$. Hence, the law of large numbers still holds, but the convergence will be slower due to infinite variance, see also Section 6.2.

6. Simulation experiments

In this section, we perform some simulations that ought to complement the theoretical results that we obtained. First, we focus on Theorem 1 and investigate how well the asymptotic normal distribution approximates the finite-sample distribution of the plug-in estimator for the mutual information and Hellinger distance in case of a Gaussian copula model. Second, we numerically assess how well the estimator (15) performs in the setting of Example 6, and compare the numerical quality of $\widehat{\mathcal{D}}$ given in (19) with the numerical quality of $\widetilde{\mathcal{D}}$ given in (17) within a bivariate Gumbel copula model.

6.1. Asymptotic normality under Gaussian copula model

Recall the asymptotic normality result of Theorem 1 for the estimator $\widehat{\mathcal{D}}^{\mathcal{N}}_{\Phi,n}=\mathcal{D}^{\mathcal{N}}_{\Phi}(\widehat{\mathbf{R}}_{n})$, with $\widehat{\mathbf{R}}_{n}$ the matrix of sample normal scores rank correlation coefficients (8) and $\mathcal{D}^{\mathcal{N}}_{\Phi}$ given in (5) and (6) for $\Phi(t)=t\log(t)$ and $\Phi(t)=(\sqrt{t}-1)^{2}$ respectively. We first turn our attention to Example 4 once more. In Figure 3, we depicted the asymptotic standard deviation $\zeta_{t\log(t)}$ and $\zeta_{(\sqrt{t}-1)^{2}}/2$. If we generate $N$ samples from e.g. a four dimensional multivariate Gaussian distribution with mean zero and covariance matrix $\mathbf{R}$ as in Example 4, we obtain $N$ estimates of $\mathcal{D}_{\Phi}^{\mathcal{N}}$, say $\widehat{\mathcal{D}}^{\mathcal{N},(1)}_{\Phi,n},\dots,\widehat{\mathcal{D}}^{\mathcal{N},(N)}_{\Phi,n}$. An empirical Monte Carlo version of $\sqrt{n}\text{SD}(\widehat{\mathcal{D}}_{\Phi,n}^{\mathcal{N}})$, with SD being the standard deviation, is then given by

[TABLE]

which we can compute for different values of $\rho_{1}$ and $\rho_{2}$, see Figure 5 for some plots in case $n=10\hskip 1.42271pt000$ and $N=1000$. Comparing with Figure 3 of Example 4, this empirically verifies the formula for the asymptotic variance in Theorem 1 in this particular setting. Kernel density estimates for the density of $\widehat{\mathcal{D}}^{\mathcal{N}}_{\Phi,n}$ when $\rho_{1}=\rho_{2}=0.5$ are included as well.

Let now $\widehat{\zeta}_{\Phi,n}$ be the plug-in estimator of the asymptotic standard deviation $\zeta_{\Phi}$ obtained by using $\widehat{\mathbf{R}}_{n}$ instead of the true $\mathbf{R}$. Based on a sample $\mathbf{X}^{(\ell)}$ for $\ell=1,\dots,n$ from a certain multivariate distribution having a Gaussian copula, we are able to compute one realization of the actual sampling distribution of the studentized estimator $\sqrt{n}\left(\widehat{\mathcal{D}}^{\mathcal{N}}_{\Phi,n}-\mathcal{D}^{\mathcal{N}}_{\Phi}(\mathbf{R})\right)/\widehat{\zeta}_{\Phi,n}$, and several replications will give an idea about the entire distribution, which should, according to Theorem 1 approximately be a standard normal one for larger values of $n$. We consider four settings which we can generate samples from:

• •

Setting 1: $k=2,d_{1}=d_{2}=2$, with standard normal marginals and a Gaussian copula having an autoregressive AR(1) correlation matrix with $\rho=0.25$.

• •

Setting 2: as Setting 1, but now with marginals

• $\star$

a $t$ distribution with $3$ degrees of freedom for $X_{11}$

• $\star$

an exponential distribution with mean $1$ for $X_{12}$

• $\star$

a beta distribution with parameters $2$ and $2$ for $X_{21}$

• $\star$

an $F$-distribution with degrees of freedom $2$ and $6$ for $X_{22}$.

• •

Setting 3: similar to Setting 1, but with $\rho=0.8$.

• •

Setting 4: $k=5,d_{1}=4,d_{2}=5,d_{3}=3,d_{4}=1,d_{5}=2$, with standard normal marginals and a Gaussian copula having an equicorrelated correlation matrix with $\rho=0.5$.

Each time, we draw $1000$ samples of sizes $n=50,200,1000,5000$ and make normal Q-Q plots to assess the goodness-of-fit with a standard normal distribution. See Figure 6 for the results of the mutual information, and Figure 7 for the half Hellinger distance.

In each setting, we have a qualitative normal approximation for larger sample sizes. Sampling from a multivariate normal distribution or from a multivariate normal copula with various marginals does not give a significant difference (Setting 1 versus Setting 2). For rather small correlations (Setting 1 and Setting 2), we observe a more pronounced lack-of-fit than for higher correlations (Setting 3). Increasing the total dimension to $q=15$ (Setting 4) results in a large positive bias for small sample size, which is no shock since empirical covariance matrices tend to be more biased when the number of parameters to estimate increases.

6.2. Nested Archimedean copula model

In the context of a general parametric copula family, the estimator (15) relies on an estimator $\widehat{\bm{\theta}}_{n}$ for the copula parameter $\bm{\theta}_{C}$ on the one hand, and on a numerical integral approximation on the other hand. In Example 6, we mathematically illustrated that, in a bivariate Gumbel copula with parameter $\bm{\theta}_{C}=\theta_{0}=3$, the estimator $\widehat{\mathcal{D}}_{(\sqrt{t}-1)^{2},n,M_{n}}$ is doomed to have a slower convergence than the estimator $\widetilde{\mathcal{D}}_{(\sqrt{t}-1)^{2},n,M_{n}}$ relying on a simplified integral approximation as in (17).

Via a small simulation, we now compare the performance of these two estimators when $\theta_{0}=3$ and $n=200$, based on $1000$ samples drawn from this bivariate Gumbel copula. We consider $M_{200}=100,1000,10\hskip 1.42271pt000$ and compare in Table 2 the sample bias, variance and mean squared error of $\widehat{\mathcal{D}}_{(\sqrt{t}-1)^{2},n,M_{n}}/2$ and $\widetilde{\mathcal{D}}_{(\sqrt{t}-1)^{2},n,M_{n}}/2$. The estimator $\widehat{\bm{\theta}}_{n}$ is based on maximizing the pseudo likelihood with non-parametric marginals, i.e. $\widehat{\bm{\theta}}_{C,n}^{\text{NP}}$, using a starting value of $1$. The true value of the dependence coefficient equals $\mathcal{D}_{(\sqrt{t}-1)^{2}}/2=0.20528$, and was computed, not by doing an empirical mean approximation of the two dimensional integral, but using numerical integration. As expected, the performance of $\widehat{\mathcal{D}}_{(\sqrt{t}-1)^{2},n,M_{n}}/2$ is poor due to a large variance that needs very large values of $M_{n}$ to go down.

Next, we investigate how well the estimator $\widetilde{\mathcal{D}}_{(\sqrt{t}-1)^{2},n,M_{n}}/2$ performs in terms of increasing $n$ for the half Hellinger distance between $(X_{1},X_{2})$ and $(X_{3},X_{4})$ having the four dimensional copula given in (18) in Example 6. The estimator $\widehat{\mathcal{D}}_{t\log(t),n,M_{n}}$ for the mutual information is also considered. We look at two settings:

• •

Setting 1: $\bm{\theta}_{C}=(\theta_{0},\theta_{1},\theta_{2})=(1,3,4)$, such that $\mathcal{D}_{(\sqrt{t}-1)^{2}}/2=\mathcal{D}_{t\log(t)}=0$, since $(X_{1},X_{2})$ and $(X_{3},X_{4})$ are independent.

• •

Setting 2: $\bm{\theta}_{C}=(\theta_{0},\theta_{1},\theta_{2})=(3,3,4)$, such that $\mathcal{D}_{(\sqrt{t}-1)^{2}}/2=0.29007$ and $\mathcal{D}_{t\log(t)}=0.99935$, computed via the true $\bm{\theta}_{C}$ and $10\hskip 1.42271pt000\hskip 1.42271pt000$ samples to numerically approximate the integral.

We take $\widehat{\bm{\theta}}_{n}=(\widehat{\theta}_{0,n},\widehat{\theta}_{1,n},\widehat{\theta}_{2,n})=\widehat{\bm{\theta}}_{C,n}^{\text{NP}}$, with $2$ as starting value for $\theta_{0}$ for maximizing the likelihood. The starting values for $\theta_{1}$ and $\theta_{2}$ are taken as the maximizers of the individual pseudo likelihoods (also with non-parametric marginals) corresponding to the marginal samples of $(X_{1},X_{2})$ and $(X_{3},X_{4})$ respectively, with both starting values equal to $2$. In each setting, we take $1000$ Monte Carlo runs and sample sizes $n=50,200,1000,5000$ and fix $M_{n}=10\hskip 1.42271pt000$ for every $n$.

Boxplots are shown in Figure 8. For each dependence coefficient and in each setting, the bias and variance tend to zero when the sample size increases. Table 3 shows empirical Monte Carlo variances $n\widehat{\text{Var}}$ of each estimator. They indicate that cases of stronger dependence (Setting 2) are harder to estimate than cases of weak dependence (Setting 1) where the asymptotic variance is smaller (as we have seen for some Gaussian examples too).

7. A financial application

Quantifying the strength of relationship between variables is fundamental in finance, e.g. in portfolio management. Individual constituents are often (positively) related to one another because of contingency on macro-economic factors, known as systematic risk. For instance, market downturns can have detrimental consequences on the portfolio as association between assets can significantly increase. This phenomenon is known as asymmetric dependence, see e.g. Alcock and Satchell (2018). Closely related is the concept of financial contagion, e.g. Gallegati (2012), Celik (2012), Wang and Hui (2017) and Akhtaruzzaman et al. (2021), among others, evidencing stronger linkages across markets in times of recession.

These markets can be considered within one and the same region, or can be spread over multiple different regions. We might for instance have a random vector $\mathbf{X}_{1}$ describing equity indexes in North America and look at the intra-dependence, or investigate the inter-dependence with European indexes $\mathbf{X}_{2}$, neglecting the within region dependence. Here, we will use the $\Phi$-dependence measures between $\mathbf{X}_{1}$ and $\mathbf{X}_{2}$, intending to illustrate cross-regional financial contagion during the COVID-19 pandemic.

In particular, we analyse historical daily logarithmic returns of stock indexes from North-America (US S&P500, Canadian S&P/TSX Composite Index and Mexican IPC Index), South-America (Brazilian IBOVESPA and Argentina Merval Index), Europe (Euronext 100, German GDAXI, Spanish IBEX 35 and Norwegian OMX Index) and Asia (Japanese Nikkei 225, Chinese SSE Composite Index, Indian S&P BSE 500 and Hong Kong HSI Index), over a time span of Dec 07, 2016 to Dec 06, 2022. The data can freely be accessed and downloaded at https://finance.yahoo.com/. Notice that logarithmic returns are i.i.d. when assuming a random walk market. Each set of index returns in each continent is considered as a random vector $\mathbf{X}_{i}$, and inter-regional financial dependence is assessed as $\mathcal{D}_{\Phi}(\mathbf{X}_{i},\mathbf{X}_{j})$ for $i\neq j$.

We first make use of the method discussed in Section 4, relying on the assumption of a Gaussian copula model. A primary impression of inter-regional financial contagion during the COVID-19 period is obtained by looking at the dependence $\mathcal{D}^{\mathcal{N}}_{\Phi}(\mathbf{X}_{i},\mathbf{X}_{j})$ over time. In total, we have $1099$ non-missing log-returns in the given period that are common for each stock index. We divide this data into $99$ windows of size $101$ with slide step $10$, i.e. $[1,101],[11,111],\dots,[981,1081]$ and a final $100$’th window $[991,1099]$. Next, we compute $\widehat{\mathcal{D}}^{\mathcal{N}}_{\Phi,n}(\mathbf{X}_{i},\mathbf{X}_{j})$ for the returns in each window and assign this to the date corresponding to the left bound of that window. As such, we get an idea of how the dependence between groups of stock indexes across different continents evolved over time, keeping in mind that a certain date reflects the dependence calculated from the $100$ future days available in the dataset. Thanks to Theorem 1, we can also add approximated confidence bounds $\widehat{\mathcal{D}}^{\mathcal{N}}_{\Phi,n}\pm z_{1-\frac{\alpha}{2}}\widehat{\zeta}_{\Phi,n}/\sqrt{n}$ for each window, with $z_{1-\frac{\alpha}{2}}$ the $1-\frac{\alpha}{2}$ lower quantile of a standard normal distribution.

Figure 9 shows the results for the mutual information and half Hellinger distance with $\alpha=0.05$. In all cases, we observe a (for some more pronounced than others) hump during the first year of the pandemic, attracting our attention. In the first half of $2020$, the Corona pandemic gave rise to a stock market crash. Many indexes around the world were recovered by the end of $2020$. See Akhtaruzzaman et al. (2021) and references therein for more detailed information. We define the pre-crisis (period $1$) to be the period Dec 07, 2016 to Nov 29, 2019, the crisis period (period $2$) Dec 02, 2019 to Dec 30, 2020, and the post-crisis period (period $3$) as Jan 04, 2021 to Dec 06, 2022. For these respective periods, there are $n_{1}=548,n_{2}=198$ and $n_{3}=352$ observed log-returns available. Denote $\mathcal{D}^{\mathcal{N}}_{\Phi,m}(\mathbf{X}_{i},\mathbf{X}_{j})$ with $m=1,2,3$ for the Gaussian copula $\Phi$-dependence between $\mathbf{X}_{i}$ and $\mathbf{X}_{j}$ in period $m$, with asymptotic standard deviation $\zeta^{\mathcal{N}}_{\Phi,m}(\mathbf{X}_{i},\mathbf{X}_{j})$, and corresponding sample versions $\widehat{\mathcal{D}}^{\mathcal{N}}_{\Phi,n_{m}}(\mathbf{X}_{i},\mathbf{X}_{j})$ and $\widehat{\zeta}_{\Phi,n_{m}}(\mathbf{X}_{i},\mathbf{X}_{j})$. A test for financial contagion is

[TABLE]

whose rejection provides statistical evidence for stronger linkages across markets during the crisis than before, and

[TABLE]

whose rejection illustrates weaker inter-regional dependence after the crisis than during. Asymptotic approximate $p$-values for these test are respectively given by

[TABLE]

and

[TABLE]

where $z_{12}$ and $z_{23}$ are the corresponding test values.

For two continents neither of which is Asia, we have statistical evidence for inter-regional financial contagion at significance level $1\%$. When Asia is included, we find larger $p$-values, especially for $p_{23}$.

Over the past year, vaccination campaigns in Asia have been losing efficiency, healthcare systems have been struggling, and putting an end to the Chinese Zero-COVID policy has led to a new flare-up.

We now look at the vector dependence between North America, South America and Europa simultaneously, not taking Asia into account. Figure 10 shows the estimated mutual information and Hellinger distance (we multiplied the half Hellinger distance with a factor $2$ in order to see the line more clearly) on the hand hand using the Gaussian copula approach (left) and on the other hand by fitting a nine dimensional nested Clayton copula (right) via the copula pseudo likelihood with non-parametrically estimated marginals using $0.1$ as starting value for the parameter of the root copula, and maximum likelihood estimates (also pseudo likelihood and non-parametric marginals) for the parameters of the child copulas, obtained with all three starting values equal to $2$. We opt for the Clayton family, as possible lower tail dependence might be incorporated as well in this copula family. Using this latter estimation approach, i.e. $\widehat{\mathcal{D}}_{t\log(t),n,M_{n}}$ or $\widetilde{\mathcal{D}}_{(\sqrt{t}-1)^{2},n,M_{n}}$, where we took $M_{n}=10\hskip 1.42271pt000$, the hump of increased dependence during the COVID-19 financial recession is even more expressed.

Finally, notice that returns can also be considered on a different periodicity, such as monthly returns. Longer horizons often allow to better capture the general trend as they are not affected so much by noise that might be present on a daily basis. For instance, as we are considering stock indexes from across the entire globe, trading hours depend on the geographic location and differences in opening and closing hours of the stock exchange might impact the dependence structure.

8. Discussion

In this paper, we started from the general objective of quantifying dependence between a finite, arbitrary amount of random vectors. We postulated desirable properties giving preference to a copula based perspective, and proved conformity for our proposed family of $\Phi$-dependence measures. A notion of perfect dependence is not among the axioms, but a characterization of it is possible for certain choices of $\Phi$, and was illustrated by means of examples.

Assuming a Gaussian copula model, an asymptotic normality result for the suggested plug-in estimator was established and could be interpreted further for specific $\Phi$-functions. Extensions to general parametric copula families were also obtained, focusing on maximum likelihood estimation and again a plug-in approach. Special attention went to nested Archimedean copulas, allowing for different parameters to control the intra- and inter-vector dependence. Simulations investigated finite-sample performances. In a real data analysis, the estimates indicate substantial levels of financial contagion during the COVID-19 recession.

Among the interests for further research are high-dimensional settings ($q$ large), in which regularisation techniques are obliging. Also, semi- or non-parametric modelling of the copula can offer more flexible dependence capturing, and as such dig deeper than associations restricted to correlations (as in Gaussian copulas) and tail dependence (as in Archimedean copulas). These issues will be studied by the authors in future research.

Acknowledgement. The authors gratefully acknowledge support from the Research Fund KU Leuven [C16/20/002 project].