∎

11institutetext: A. Roy 22institutetext: Indian Statistical Institute
Applied Statistics Unit
22email<EMAIL_ADDRESS>33institutetext: A. K. Ghosh
44institutetext: Indian Statistical Institute
Stat Math Unit
44email<EMAIL_ADDRESS>55institutetext: A. Goswami 66institutetext:
Indian Statistical Institute
Stat Math Unit
66email<EMAIL_ADDRESS>77institutetext: C. A. Murthy 88institutetext:
Indian Statistical Institute
Machine Intelligence Unit
88email<EMAIL_ADDRESS>
# Some New Copula Based Distribution-free Tests of Independence among Several
Random Variables

Angshuman Roy  Anil K. Ghosh  Alok Goswami  C. A. Murthy

(Received: date / Accepted: date)

###### Abstract

Over the last couple of decades, several copula based methods have been
proposed in the literature to test for independence among several random
variables. But these existing tests are not invariant under monotone
transformations of the variables, and they often perform poorly if the
dependence among the variables is highly non-monotone in nature. In this
article, we propose a copula based measure of dependency and use it to
construct some distribution-free tests of independence. The proposed measure
and the resulting tests, all are invariant under permutations and strictly
monotone transformations of the variables. Our dependency measure involves a
kernel function with an associated bandwidth parameter. We adopt a multi-scale
approach, where we look at the results obtained for several choices of the
bandwidth and aggregate them judiciously. Large sample properties of the
dependency measure and the resulting tests are derived under appropriate
regularity conditions. Several simulated and real data sets are analyzed to
compare the performance of the proposed tests with some popular tests
available in the literature.

###### Keywords:

Gaussian kernel, Invariance, Large sample consistency, Measure of dependency,
Multi-scale approach

##  1 Introduction

Measuring and testing for dependence among $d~{}(d\geq 2)$ random variables is
a classical problem in statistics. For $d=2$, Pearson’s product-moment
correlation coefficient, being arguably the most popular measure of
dependence, has been used to construct test of independence between two random
variables (see, e.g., Anderson, 2003). But this measure is sensitive against
outliers and extreme values, and it often fails to capture non-linear
dependence between the variables. Other popular measures of association like
Spearman’s rank correlation coefficient $\rho$ (Spearman, 1904), Kendall’s
measure of association $\tau$ (Kendall, 1938), Blomqvist’s quadrant statistic
$\beta$ (Blomqvist, 1950), are robust against outliers. They can also detect
monotone or near monotone relationship between the variables. Tests based on
these rank based statistics have the distribution-free property under the null
hypothesis of independence. Hoeffding (1948) also developed a distribution-
free test of independence using $\phi$ statistic based on bivariate copula.
Reshef et al. (2011) proposed maximal information coefficient as a measure of
dependency, but the tests based on this measure usually have low powers. Tests
of independence between two random vectors include the work of Gieser and
Randles (1997); Taskinen et al. (2003, 2005); Heller et al. (2012, 2013);
Biswas et al. (2016); Sarkar and Ghosh (2018). Székely et al. (2007) developed
a test of independence based on distance correlation, which is known as the
dCov test . Gretton et al. (2007) considered a test based on the Hilbert-
Schmidt norm of the cross-covariance operator, which is popularly known as the
HSIC test.

Um and Randles (2001) generalized Giesser and Randles’ (1997) test for several
variables. Gaißer et al. (2010) developed multivariate extensions of
Hoeffding’s (1948) $\phi$ statistic and the associated test. Pfister et al.
(2017) and Fan et al. (2017) proposed multivariate extensions of the HSIC test
and the dCov test (referred to as the dHSIC test and the mdCov test),
respectively.

Using the ideas of copula, Nelsen (1996) and Úbeda-Flores (2005) proposed
multivariate extensions of Spearman’s $\rho$, Kendall’s $\tau$ and Blomqvist’s
$\beta$ statistics. Tests of independence based on these dependency measures
have the distribution-free property, but they often yield poor performance
when the relationships among the variables are highly non-monotone in nature.
To take care of this problem, in this article, we propose a new copula based
multivariate measure of dependency and use it to test mutual independence
among several random variables.

Our work is motivated by Póczos et al. (2012), where the authors proposed a
dependency measure, which is non-negative and takes the value $0$ if and only
if the random variables $X_{1},X_{2},\ldots,X_{d}$ are jointly independent.
This measure is also invariant under permutation of the coordinates. Our
proposed measure satisfies all these desirable properties. Moreover, while the
dependency measure of Póczos et al. (2012) is only invariant under strictly
increasing transformations of the $X_{i}$’s, our measure is invariant under
strictly monotone transformations. It also satisfies several other desirable
properties. For instance, in the bivariate case, it satisfies all Dependence
Axioms proposed by Schweizer et al. (1981). We propose a data driven estimator
of this measure, which also enjoys similar desirable properties. Unlike Póczos
et al. (2012), here one does not need to generate observations from a uniform
distribution for its construction. Our dependency measure and its estimator
involve a kernel function, and we use the Gaussian kernel for this purpose.

One can use this estimator to a develop distribution-free test. However, the
performance of the test depends on the bandwidth parameter associated with the
Gaussian kernel. While larger bandwidths work well for near monotone
relationships (i.e., when the conditional expectation of one variable is
nearly a monotone function of others) among the variables, smaller bandwidths
are preferred to detect non-monotone relations. So, borrowing idea from multi-
scale classification (see, e.g., Ghosh et al., 2006), here we adopt a multi-
scale approach, where we look at the results for various choices of bandwidth
and then aggregate them judiciously to arrive at the final decision.

The organization of the paper is as follows. In Section 2, we define our
copula based dependency measure and derive some of its desirable properties.
In particular, we prove its invariance over permutations and monotone
transformations of the variables. In Section 3, we propose a nonparametric
estimate of this dependency measure and investigate its theoretical
properties. Some distribution-free tests based on this estimate are
constructed in Section 4, where we also establish the large sample consistency
of these tests. Several simulated and real data sets are analyzed in Section 5
to compare the performance of the proposed tests with some popular tests
available in the literature. Section 6 contains a brief summary of the work
and ends with a discussion on some possible directions for future research.
All proofs and mathematical details are given in Appendix Section.

##  2 The proposed measure of dependency

Our measure of dependency is based on the copula distribution of a
$d$-dimensional random vector. A $d$-dimensional copula is a probability
distribution $C$ on the $d$-dimensional unit hypercube $[0,1]^{d}$ such that
all of its one-dimensional marginals are uniform on $[0,1]$. If $F$ is the
distribution function of a $d$-dimensional random vector
$\mathbf{X}=(X_{1},\ldots,X_{d})$ with continuous one-dimensional marginals
$F_{1},\cdots,F_{d}$, then the copula transformation of $F$ or the copula
distribution of $\mathbf{X}$ is given by

$\displaystyle\mathsf{C}_{\mathbf{X}}(\mathbf{u})=F(F_{1}^{-1}(u_{1}),F_{2}^{-1}(u_{2}),\cdots,F_{d}^{-1}(u_{d})),$
(1)

where $\mathbf{u}=(u_{1},u_{2},\ldots,u_{d})\in[0,1]^{d}$ and
$F_{i}^{-1}(u_{i})\\!=\\!\inf\\{x\\!:\\!F_{i}(x)\\!>\\!u_{i}\\}$ for all
$i=1,2,\ldots,d$ (see, e.g., Nelsen, 2013, for further discussion on copula).
If $\mathsf{C}_{\mathbf{X}}$ is the cumulative distribution function of a
uniform distribution on $[0,1]^{d}$, i.e., $X_{1},\ldots,X_{d}$ are
independent, the associated copula is called the uniform copula and denoted by
$\Pi$. On the other hand, if $X_{1},\ldots,X_{d}$ are comonotonic, i.e. there
exist strictly increasing functions $g_{i}$’s and a random variable $V$ such
that $\mathbf{X}\stackrel{{\scriptstyle
d}}{{=}}(g_{1}(V),g_{2}(V),\ldots,g_{d}(V))$, it is called the maximum copula
and denoted by ${\mathsf{M}}$. So, for every $\mathbf{u}\in[0,1]^{d}$, we have
$\Pi(\mathbf{u})=\prod_{i=1}^{d}u_{i}$ and
${\mathsf{M}}(\mathbf{u})=\min\\{u_{1},\ldots,u_{d}\\}$.

Naturally, larger difference between $\mathsf{C}_{\mathbf{X}}$ and $\Pi$
indicates higher degree of dependence among $X_{1},\ldots,X_{d}$. To measure
the difference between two probability distributions $\mathbb{P}$ and
$\mathbb{Q}$ on $\mathbb{R}^{d}$, we use

$\displaystyle\gamma_{k}(\mathbb{P},\mathbb{Q})=\Bigl{[}\mathrm{E}\left\\{k(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})+k(\mathbf{Z},\mathbf{Z}^{{}^{\prime}})-2k(\mathbf{Y},\mathbf{Z})\right\\}\Bigr{]}^{1/2};$
(2)

where $\mathbf{Y},\mathbf{Y}^{{}^{\prime}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}{\mathbb{P}}$,
$\mathbf{Z},\mathbf{Z}^{{}^{\prime}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}{\mathbb{Q}}$ are independent, and
$k:\mathbb{R}^{d}\times\mathbb{R}^{d}$ $\rightarrow\mathbb{R}$ is a symmetric,
bounded, positive definite kernel. It is known that $\gamma_{k}$ is a pseudo-
metric on the space of all continuous probability distributions on
$\mathbb{R}^{d}$, and it is a metric when $k$ is a characteristic kernel (see
Fukumizu et al., 2008)). Gaussian kernel
$k_{\sigma}(x,y)=\exp\left(-\frac{\|x-y\|^{2}}{2\sigma^{2}}\right)$ with a
bandwidth parameter $\sigma>0$ is a popular choice as a characteristic kernel,
and we shall use it throughout this article.

From the above discussion, it is clear that for any characteristic kernel $k$
on $\mathbb{R}^{d}\times\mathbb{R}^{d}$, one can use
$\gamma_{k}({\mathsf{C}}_{\mathbf{X}},\Pi)$ or
$\gamma^{2}_{k}({\mathsf{C}}_{\mathbf{X}},\Pi)$ as a measure of dependency
among the coordinate variables $X_{1},X_{2},\ldots,X_{d}$. In this article, we
use a scaled version of this measure given by

$\displaystyle
I_{\sigma}(\mathbf{X})={{\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{X}},\Pi)}}/{{\gamma_{k_{\sigma}}(\mathsf{M},\Pi)}}.$
(3)

Note that the denominator $\gamma_{k_{\sigma}}(\mathsf{M},\Pi)$ is strictly
positive. So, $I_{\sigma}(\mathbf{X})$ is well defined. The use of the
Gaussian kernel makes the measure $I_{\sigma}(\mathbf{X})$ invariant under
permutations and strictly monotone transformations of the coordinate
variables. The result is stated below.

###### Lemma 1

$I_{\sigma}(\mathbf{X})$ is invariant under permutations and strictly monotone
transformations of $X_{1},X_{2},\ldots,X_{d}$.

From the definition of $I_{\sigma}(\mathbf{X})$, it is clear that it takes the
value $0$ if and only if the coordinates of $\mathbf{X}$ are independent, and
its value is supposed to increase as the dependency among
$X_{1},X_{2},\ldots,X_{d}$ increases. The following lemma shows that in case
of extreme dependency (i.e., when for each pair of variables, one is a
strictly monotone function of the other), it turns out to be $1$.

###### Lemma 2

For all $i=1,\ldots,d-1$, if $X_{i}$ is a strictly monotone function of
$X_{d}$ with probability one, then $I_{\sigma}(\mathbf{X})$ takes the value
$1$.

This desirable property of $I_{\sigma}(\mathbf{X})$ helps us to properly
assess the degree of dependency among $X_{1},\ldots,X_{d}$. Note that many
well-known dependency measures like the copula based multivariate extensions
of Spearman’s $\rho$, Kendall’s $\tau$, Blomqvist’s $\beta$ and Hoeffding’s
$\phi$ statistic (see, e.g., Úbeda-Flores, 2005; Nelsen, 1996; Gaißer et al.,
2010) do not have this property. But, since $I_{\sigma}(\mathbf{X})$ always
takes value in $[0,1]$, in the case of $d=2$, unlike Spearman’s $\rho$,
Kendall’s $\tau$ statistic or Blomqvist’s $\beta$ statistic, it does not give
us any idea about the direction of dependence between two variables. We know
that the distance correlation measure proposed by Székely et al. (2007) can be
expressed as a weighted squared distance between the characteristic functions
of two distributions. The following theorem shows that
$I_{\sigma}(\mathbf{X})$ also has a similar property.

###### Theorem 2.1

Let $\varphi_{\mathsf{C}_{\mathbf{X}}}$ and $\varphi_{\Pi}$ be the
characteristic functions of $\mathsf{C}_{\mathbf{X}}$ and $\Pi$, respectively.
Define
$C_{\sigma,d}=\kappa\left(\frac{\sigma}{\sqrt{d}}\right)+\kappa^{d}(\sigma)-2\int\displaylimits_{0}^{1}\lambda^{d}(u,\sigma)\,du,$
where
$\kappa(\sigma)=\sqrt{2\pi}\sigma\left[2\Phi\left(\frac{1}{\sigma}\right)-1\right]-2\sigma^{2}\left[1-\exp\left(-\frac{1}{2\sigma^{2}}\right)\right]$,
$\lambda(x,\sigma)=\sqrt{2\pi}\sigma\left[\Phi\left(\frac{x}{\sigma}\right)+\Phi\left(\frac{1-x}{\sigma}\right)-1\right]$
and $\Phi(\cdot)$ is the cumulative distribution function of the N(0,1)
distribution. Then $I^{2}_{\sigma}(\mathbf{X})$ can be expressed as

$I^{2}_{\sigma}(\mathbf{X})=C_{\sigma,d}^{-1}\left(\frac{\sigma}{\sqrt{2\pi}}\right)^{d}\int_{\mathbb{R}^{d}}|\varphi_{\mathsf{C}_{\mathbf{X}}}(\mathbf{w})-\varphi_{\Pi}(\mathbf{w})|^{2}\exp\left(-\frac{\sigma^{2}}{2}\mathbf{w}^{\intercal}\mathbf{w}\right)\,d\mathbf{w}.$

Another interesting property of $I_{\sigma}(\mathbf{X})$ is its continuity.
Note that if $\\{\mathbf{X}_{n};n\geq 1\\}$ is a sequence of random vectors
converging in distribution to $\mathbf{X}$, then $\mathsf{C_{\mathbf{X}_{n}}}$
converges to $\mathsf{C}_{\mathbf{X}}$ weakly. So, using the dominated
convergence theorem, from Theorem 1 it follows that
$I_{\sigma}(\mathbf{X}_{n})$ converges to $I_{\sigma}(\mathbf{X})$ as $n$
increases. This result is stated below.

###### Lemma 3

Let $\\{\mathbf{X}_{n}:n\geq 1\\}$ be a sequence of $d$-dimensional random
vectors with continuous one-dimensional marginals. If $\mathbf{X}_{n}$
converges to $\mathbf{X}$ weakly, we have
$\lim_{n\to\infty}I_{\sigma}(\mathbf{X}_{n})=I_{\sigma}(\mathbf{X}).$

In the case of $d=2$, $I_{\sigma}(\mathbf{X})$ enjoys some additional
properties. For instance, $I^{2}_{\sigma}(\mathbf{X})$ can be viewed as a
product moment correlation coefficient between two random quantities. If
$\mathbf{X}$ follows a bivariate normal distribution with correlation
coefficient $r$, $I_{\sigma}(\mathbf{X})$ turns out to be a strictly
increasing function of $|r|$. These results are shown by the following
theorem.

###### Theorem 2.2

Let $\mathbf{X}=(X_{1},X_{2})$ be a bivariate random vector with continuous
one-dimensional marginals.

  1. $(a)$

Suppose that $\mathbf{T}=(T_{1},T_{2})$ and
$\mathbf{T}^{{}^{\prime}}=(T_{1}^{{}^{\prime}},T_{2}^{{}^{\prime}})$ are
independent, and they follow the distribution $\mathsf{C}_{\mathbf{X}}$, the
copula distribution of $\mathbf{X}$. Define
$V_{i}=k_{\sigma}(T_{i},T_{i}^{{}^{\prime}})-\mathrm{E}\left[k_{\sigma}(T_{i},T_{i}^{{}^{\prime}})\Big{|}T_{i}\right]-\mathrm{E}\left[k_{\sigma}(T_{i},T_{i}^{{}^{\prime}})\Big{|}T_{i}^{{}^{\prime}}\right]+\mathrm{E}\left[k_{\sigma}(T_{i},T_{i}^{{}^{\prime}})\right]$
for $i=1,2$. Then we have
$I^{2}_{\sigma}(\mathbf{X})=\mathrm{Cor}(V_{1},V_{2})$, which takes the value
$1$ if only if $X_{1}$ is a strictly monotone function of $X_{2}$ with
probability one.

  2. $(b)$

If $\mathbf{X}$ follows a bivariate normal distribution with
$\mathrm{Cor}(X_{1},X_{2})=r$, then $I_{\sigma}(\mathbf{X})$ is a strictly
increasing function of $|r|$ with $I_{\sigma}(\mathbf{X})\leq|r|.$

Another interesting property of $I_{\sigma}(\mathbf{X})$ is its
irreducibility. Following Schmid et al. (2010), we call a dependency measure
$I$ irreducible if, for any $d>2$, $I(X_{1},\ldots,X_{d})$ is not a function
of the quantities
$\\{I(X_{i_{1}},\ldots,X_{i_{k}}):\\{i_{1},\cdots,i_{k}\\}\subsetneqq\\{1,\ldots,d\\}\\}$.
Naturally, any reasonable multivariate measure of dependency is expected to be
irreducible. Note that if $I(X_{1},X_{2},X_{3})$ gets completely determined by
$I(X_{1},X_{2})$, $I(X_{2},X_{3})$ and $I(X_{3},X_{1})$, instead of mutual
dependence among $X_{1}$, $X_{2}$ and $X_{3}$, it can only detect pairwise
dependence. The following theorem shows that any copula based multivariate
dependency measure, which takes the value zero only for the uniform copula, is
irreducible.

###### Theorem 2.3

Let $\mathsf{C}_{\mathbf{X}}$ be the copula distribution of $\mathbf{X}$ and
$I(\mathbf{X})=\mathcal{M}(C_{\mathbf{X}})$ be a copula based multivariate
dependency measure. If $\mathcal{M}(\mathsf{C}_{\mathbf{X}})=0$ implies
$\mathsf{C}_{\mathbf{X}}=\Pi$, then $I$ is irreducible.

For any fixed bandwidth parameter $\sigma$, the irreducibility of our proposed
measure $I_{\sigma}(\mathbf{X})$ follows from Theorem 2.3 as a corollary.
However, this property vanishes when $\sigma$ diverges to infinity. In such a
situation, the limiting value of $I_{\sigma}(\mathbf{X})$ turns out to be the
average of squared Spearman’s rank correlations between $\binom{d}{2}$ pairs
of random variables as stated in the following theorem.

###### Lemma 4

As $\sigma$ diverges to infinity, $I^{2}_{\sigma}(\mathbf{X})$ converges to
$\frac{1}{{d\choose 2}}\sum_{1\leq i<j\leq d}\mathrm{Cor}^{2}(S_{i},S_{j})$,
where $\mathbf{S}=(S_{1},S_{2},\cdots,S_{d})\sim\mathsf{C}_{\mathbf{X}}$.

##  3 Estimation of the proposed measure

Let $\mathbf{X}^{(1)},\cdots,\mathbf{X}^{(n)}$ be $n$ independent copies of
the random vector $\mathbf{X}$ taking values in ${\mathbb{R}}^{d}$. For any
fixed $j=1,\ldots,d$, define $R_{j}^{(i)}$ as the rank of $X_{j}^{(i)}$
($i=1,\ldots,n$) in the set $\\{X_{j}^{(1)},\ldots,X_{j}^{(n)}\\}$ to get
$\mathbf{R}^{(i)}=(R_{1}^{(i)},\ldots,R_{d}^{(i)})$, the coordinate-wise rank
of $\mathbf{X}^{(i)}$. We use the normalized rank vectors
$\mathbf{Y}^{(i)}=\mathbf{R}^{(i)}/n~{}$ ($i=1,\ldots,n$) to define the
empirical version of the copula distribution $\mathsf{C}_{\mathbf{X}}$, which
is given by

$\displaystyle\mathsf{C}_{\mathbf{x},n}(u_{1},u_{2},\ldots,u_{d})=\frac{1}{n}\sum_{i=1}^{n}{\mathbb{I}}[Y_{1}^{(i)}\leq
u_{1},Y_{2}^{(i)}\leq u_{2},\ldots,Y_{d}^{(i)}\leq u_{d}],$ (4)

where ${\mathbb{I}}$ is the indicator function. Clearly,
$\mathsf{C}_{\mathbf{x},n}$ is the empirical distribution function based on
$\mathbf{Y}^{(1)},\cdots,\mathbf{Y}^{(n)}$. Indeed, it is the copula transform
of the empirical distribution based on
$\mathbf{X}^{(1)},\ldots,\mathbf{X}^{(n)}$. Similarly, we can define empirical
versions of the maximum copula and the uniform copula as

$\displaystyle\mathsf{M}_{n}(u_{1},u_{2},\ldots,u_{d})=\frac{1}{n}\sum_{i=1}^{n}{\mathbb{I}}[u_{1}\geq
i/n,u_{2}\geq i/n,\ldots,u_{d}\geq i/n]~{}\mbox{and}$ (5)
$\displaystyle\Pi_{n}(u_{1},u_{2},\cdots,u_{d})=\frac{1}{n^{d}}\sum_{1\leq
i_{1},i_{2},\cdots,i_{d}\leq n}{\mathbb{I}}[u_{1}\geq i_{1}/n,u_{2}\geq
i_{2}/n,\cdots,u_{d}\geq i_{d}/n],$

respectively. While $\mathsf{M}_{n}$ puts mass $1/n$ on each of the $n$ points
$\\{({i}/{n},{i}/{n},\cdots,{i}/{n}):$ $~{}1\leq i\leq n\\}$, $\Pi_{n}$
assigns equal mass to $n^{d}$ points of the form
$({i_{1}}/{n},{i_{2}}/{n},\cdots,{i_{d}}/{n})$, for
$i_{1},i_{2},\cdots,i_{d}\in\\{1,2,\cdots,n\\}$. We estimate
$I_{\sigma}(\mathbf{X})$ by its empirical analog

$\widehat{I}_{\sigma,n}(\mathbf{X})={\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\Pi_{n})}/{\gamma_{k_{\sigma}}(\mathsf{M}_{n},\Pi_{n})}.$
(6)

Note that $\widehat{I}_{\sigma,n}(\mathbf{X})$ is well-defined since
$\mathsf{M}_{n}\not=\Pi_{n}$ for every $n>1$. One can also check that
$\widehat{I}_{\sigma,n}(\mathbf{X})$ can be expressed as (replace expectation
by sample mean in equation (2))

$\widehat{I}_{\sigma,n}(\mathbf{X})=\sqrt{\dfrac{s_{1}-2s_{2}+v_{3}}{v_{1}-2v_{2}+v_{3}}},$
(7)

where $s_{1}=\frac{1}{n^{2}}\sum\limits_{1\leq i<j\leq
n}k_{\sigma}(\mathbf{Y}^{(i)},\mathbf{Y}^{(j)})+\frac{1}{n},\;s_{2}=\frac{1}{n^{d+1}}\sum\limits_{i=1}^{n}\prod\limits_{j=1}^{d}\sum\limits_{l=1}^{n}k_{\sigma}(Y_{j}^{(i)},\frac{1}{n}),$
$v_{1}=\frac{2}{n^{2}}\sum\limits_{i=1}^{n-1}(n-i)e^{-\frac{d}{2}\left({i}/{n\sigma}\right)^{2}}+\frac{1}{n},\;v_{2}=\frac{1}{n^{d+1}}\sum\limits_{i=1}^{n}\left[\sum\limits_{j=1}^{n}e^{-\frac{1}{2}\left({(i-j)}/{n\sigma}\right)^{2}}\right]^{d}$
and
$v_{3}=\left[\frac{2}{n^{2}}\sum\limits_{i=1}^{n-1}(n-i)e^{-\frac{1}{2}\left({i}/{n\sigma}\right)^{2}}+\frac{1}{n}\right]^{d}$.

The above formula shows that the computing cost of
$\widehat{I}_{\sigma,n}(\mathbf{X})$ is $O(dn^{2})$. This estimate enjoys some
nice theoretical properties similar to those of $I_{\sigma}(\mathbf{X})$.
These properties are mentioned below.

###### Lemma 5

Suppose that $\mathbf{X}^{(1)},\mathbf{X}^{(2)},\ldots,\mathbf{X}^{(n)}$ are
$n$ independent observations from the distribution of a $d$-dimensional random
vector $\mathbf{X}=(X_{1},X_{2},\ldots,X_{d})$ having continuous one-
dimensional marginals. Then, we have the following results.

  1. (a)

$\widehat{I}_{\sigma,n}(\mathbf{X})$ is invariant under permutation and
strictly monotone transformations of the coordinate variables
$X_{1},X_{2},\ldots,X_{d}$.

  2. (b)

For all $i=1,\ldots,d-1$, if $X_{i}$ is a strictly monotone function of
$X_{d}$ with probability one, then $\widehat{I}_{\sigma,n}(\mathbf{X})$ takes
the value $1$.

Note that other existing copula based dependency measures do not have the
property mentioned in part ($b$) of Lemma 5. For instance, multivariate
extensions of Spearman’s $\rho$, Kendall’s $\tau$, Blomqvist’s $\beta$ and
Hoeffding’s $\phi$ statistics (see, e.g., Nelsen, 1996, 2002; Gaißer et al.,
2010) may not take the value $1$ even when the measurement variables have
monotone relationships among them. To demonstrate this, we considered a simple
example. We generated 10000 observations on
$\mathbf{X}=(X_{1},X_{2},\cdots,X_{d})$, where $X_{i}=V$ or $X_{i}=-V$ for
$V\sim U(0,1)$ and $i=1,2,\ldots,d$. Hence each pair of variables were
monotonically related. We considered three choices of $d$ ($d=3,4,5$), and for
each value of $d$, results are reported in Table 1 for different types of
relationships shown in the orientation column. For example, the
$(\uparrow,\uparrow,\downarrow)$ sign in the orientation column indicates that
$(X_{1},X_{2},X_{3})=(V,V,-V)$. Table 1 clearly shows that all dependency
measures considered here take the value $1$ when the relationships among the
variables are strictly increasing. But, except for
$\widehat{I}_{\sigma,n}(\mathbf{X})$, all other measures fail to have this
property for other monotone relationships among the variables.

Dimension | Orientation | $\widehat{I}_{\sigma,n}(\mathbf{X})$ | $\widehat{\rho}$ | $\widehat{\tau}$ | $\widehat{\beta}$ | $\widehat{\varphi}$
---|---|---|---|---|---|---
3 | $(\uparrow,\uparrow,\uparrow)$ | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
3 | $(\uparrow,\uparrow,\downarrow)$ | 1.000 | -0.333 | -0.333 | -0.333 | 0.517
4 | $(\uparrow,\uparrow,\uparrow,\uparrow)$ | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
4 | $(\uparrow,\uparrow,\uparrow,\downarrow)$ | 1.000 | -0.091 | -0.143 | -0.143 | 0.382
4 | $(\uparrow,\uparrow,\downarrow,\downarrow)$ | 1.000 | -0.212 | -0.143 | -0.143 | 0.327
5 | $(\uparrow,\uparrow,\uparrow,\uparrow,\uparrow)$ | 1.000 | 1.000 | 1.000 | 1.000 | 1.000
5 | $(\uparrow,\uparrow,\uparrow,\uparrow,\downarrow)$ | 1.000 | 0.016 | -0.067 | -0.067 | 0.347
5 | $(\uparrow,\uparrow,\uparrow,\downarrow,\downarrow)$ | 1.000 | -0.108 | -0.067 | -0.067 | 0.236
Table 1: Different measures of dependency when the variables are monotonically
related.

Since $\widehat{I}_{\sigma,n}(\mathbf{X})$ is based on coordinate-wise ranks
of the observations, it is robust against contaminations and outliers
generated from heavy-tailed distributions. Following the results in Póczos et
al. (2012), one can show that addition of a new observation can change its
value by at most ${\bf O}(n^{-1})$. Just like $I_{\sigma}(\mathbf{X})$, its
empirical analog $\widehat{I}_{\sigma,n}(\mathbf{X})$ also enjoys some
additional properties for $d=2$. Theorem 3.1 below shows that result analogous
to Theorem 2.2(a) holds for $\widehat{I}_{\sigma,n}(\mathbf{X})$ as well.

###### Theorem 3.1

Suppose that $\mathbf{X}^{(1)},\cdots,\mathbf{X}^{(n)}$ are independent
observations from a bivariate distribution with continuous one-dimensional
marginals and $\mathbf{Y}^{(1)},\cdots,\mathbf{Y}^{(n)}$ are their normalized
coordinate-wise ranks. For $i,j=1,2,\ldots,n$, define

$\begin{array}[]{ll}V_{1}(i,j)=k_{\sigma}(Y_{1}^{(i)},Y_{1}^{(j)})-\frac{1}{n}\sum_{i=1}^{n}k_{\sigma}(Y_{1}^{(i)},Y_{1}^{(j)})-\frac{1}{n}\sum_{j=1}^{n}k_{\sigma}(Y_{1}^{(i)},Y_{1}^{(j)})\\\
\qquad\qquad+\frac{1}{n^{2}}\sum_{i,j=1}^{n}k_{\sigma}(Y_{1}^{(i)},Y_{1}^{(j)}),~{}\mbox{and}\\\
V_{2}(i,j)=k_{\sigma}(Y_{2}^{(i)},Y_{2}^{(j)})-\frac{1}{n}\sum_{i=1}^{n}k_{\sigma}(Y_{2}^{(i)},Y_{2}^{(j)})-\frac{1}{n}\sum_{j=1}^{n}k_{\sigma}(Y_{2}^{(i)},Y_{2}^{(j)})\\\
\qquad\qquad+\frac{1}{n^{2}}\sum_{i,j=1}^{n}k_{\sigma}(Y_{2}^{(i)},Y_{2}^{(j)}).\end{array}$

Then $\widehat{I}_{\sigma,n}(\mathbf{X})$ can be expressed as
$\widehat{I}_{\sigma,n}(\mathbf{X})=\Big{[}\sum_{i,j}V_{1}(i,j)V_{2}(i,i)\Big{]}/\Big{[}\sum_{i,j}V_{1}^{2}(i,j)\sum_{i,j}V_{2}^{2}(i,j)\Big{]}^{1/2}$.
As a consequence, we have $0\leq\widehat{I}_{\sigma,n}(\mathbf{X})\leq 1$,
where the equality holds if and only if one coordinate variable is a strictly
monotone function of the other.

##  4 Test of independence based on $\widehat{I}_{\sigma,n}(\mathbf{X})$

We have seen that $I_{\sigma}(\mathbf{X})$ serves as a measures of dependence
among the coordinates of $\mathbf{X}$. It is non-negative, and takes the value
$0$ if and only if $X_{1},X_{2},\ldots,X_{d}$ are independent. So, we can use
$\widehat{I}_{\sigma,n}(\mathbf{X})$ as the test statistic and reject $H_{0}$,
the null hypothesis of independence, for large values of
$\widehat{I}_{\sigma,n}(\mathbf{X})$. The large sample distribution of our
test statistic is given by the following theorem.

###### Theorem 4.1

Suppose that $\mathbf{X}$ follows a multivariate distribution with continuous
one-dimensional marginals. Also assume that the associated copula distribution
$\mathsf{C}_{\mathbf{X}}$ has continuous partial derivatives.

  * (a)

$\mbox{If}~{}\mathsf{C}_{\mathbf{X}}=\Pi,\mbox{
then}~{}~{}n\widehat{I}_{\sigma,n}^{2}(\mathbf{X})\overset{\mathcal{L}}{\longrightarrow}\;\sum_{i=1}^{\infty}\lambda_{i}Z_{i}^{2}$,
where the $Z_{i}$’s are i.i.d. $\mathcal{N}(0,1)$ and the $\lambda_{i}$’s are
some positive constants (see the proof of the theorem in the Appendix for
detailed description of the $\lambda_{i}$’s).

  * (b)

If $\mathsf{C}_{\mathbf{X}}\neq\Pi$, then
$\sqrt{n}(\widehat{I}_{\sigma,n}(\mathbf{X})-I_{\sigma}(\mathbf{X}))\overset{\mathcal{L}}{\longrightarrow}\mathcal{N}(0,\delta^{2})$;
where
$\delta^{2}\\!=\\!C_{\sigma,d}^{-2}I^{-2}_{\sigma}(\mathbf{X})\\!\int\limits_{[0,1]^{d}}\\!\int\limits_{[0,1]^{d}}\\!g(u)g(v)\,\mathrm{E}[\,d\mathbb{G}_{\mathsf{C}_{\mathbf{X}}}(u)\,d\mathbb{G}_{\mathsf{C}_{\mathbf{X}}}(v)]$,
$g(u)\\!=\\!\\!\int\limits_{[0,1]^{d}}\\!\\!k_{\sigma}(u,v)d(\mathsf{C}_{\mathbf{X}}\\!-\Pi)(v)$
and $\mathbb{G}_{\mathsf{C}_{\mathbf{X}}}$ is a 0 mean Gaussian process (as
defined in Theorem T1 in the Appendix).

The histograms in Figures 1(a) and 1(b) show the empirical distributions of
$\widehat{I}_{\sigma,n}(\mathbf{X})$ computed based on 5000 independent
samples, each of size 200, generated from bivariate normal distributions with
correlation coefficient $\rho_{0}=0$ and $0.5$, respectively. For
$\rho_{0}=0.5$ (i.e., $\mathsf{C}_{\mathbf{X}}\neq\Pi$) while the empirical
distribution looks like a normal distribution, for $\rho_{0}=0$ (i.e.,
$\mathsf{C}_{\mathbf{X}}=\Pi$), it turns out to be positively skewed. This is
consistent with the result stated in Theorem 4.1.

Figure 1: Empirical distribution of $\widehat{I}_{\sigma,n}(\mathbf{X})$ for
$\sigma=0.2$

The probability convergence of $\widehat{I}_{\sigma,n}(\mathbf{X})$ follows
from Theorem 4.1. But, we also have a stronger result in this context. The
following theorem shows that $\widehat{I}_{\sigma,n}(\mathbf{X})$ converges to
$I_{\sigma}(\mathbf{X})$ almost surely.

###### Theorem 4.2

$\widehat{I}_{\sigma,n}(\mathbf{X})$ converges to $I_{\sigma}(\mathbf{X})$
almost surely as the sample size $n$ tends to infinity.

From Theorem 4.2, it is clear that under the null hypothesis of independence,
$\widehat{I}_{\sigma,n}(\mathbf{X})$ converges to $0$ almost surely, while
under the alternative, it converges to a positive constant. For any fixed
choice of $\sigma$, the large sample consistency of the test follows from it.
However, for practical implementation of the test, one needs to determine the
cut-off. It is difficult to find this cut-off based on the asymptotic null
distribution of the test statistic mentioned in Theorem 4.1 since the
coefficients $\lambda_{i}$’s associated with the chi-square distributions are
all unknown. Here we use the distribution-free property of
$\widehat{I}_{\sigma,n}(\mathbf{X})$ to determine the cut-off. Note that under
$H_{0}$, for each $j=1,2,\ldots,d$, we have
$P(R_{j}^{(1)}=i_{1},R_{j}^{(2)}=i_{2},\ldots,R_{j}^{(n)}=i_{n})=1/n!$ for any
permutation $\\{i_{1},i_{2},\ldots,i_{n}\\}$ of $\\{1,2,\ldots,n\\}$, and for
different values of $j$, they are independent. So, we can easily generate
normalized coordinate-wise ranks to compute the test statistic. We repeat this
procedure 10,000 times to approximate the $(1-\alpha)$-th quantile of the null
distribution of $\widehat{I}_{\sigma,n}(\mathbf{X})$, which is then used as
the cut-off. This whole calculation can be done off-line, and we can prepare a
table of critical values for different choices $n$ and $\sigma$ before
handling actual observations.

Though any fixed choice of the bandwidth $\sigma$ leads to a consistent test
(follows from Theorem 4.2), its power may depend on this choice. The method
commonly used for choosing the bandwidth is based on “median heuristic” (see,
e.g., Fukumizu et al., 2009, Sec 5), where one computes all pairwise distances
among the observations and then the median of those distances is taken as the
bandwidth. Since we are using the kernel on the normalized rank vectors ${\bf
Y}^{(1)},{\bf Y}^{(2)},\ldots,{\bf Y}^{(n)}$ having the null distribution
$\Pi_{n}$, following this idea, we can choose $\sigma$ to be the median of
$\|\mathbf{Z}-\mathbf{Z}^{{}^{\prime}}\|$, where
$\mathbf{Z},\mathbf{Z}^{{}^{\prime}}\stackrel{{\scriptstyle
i.i.d.}}{{\sim}}\Pi_{n}$, and $\|\cdot\|$ denotes the usual Euclidean norm.
Note that the bandwidth chosen in this way is non-random and it is a function
of $n$. We denote it by $\sigma_{n}$. As $n$ increases, since $\Pi_{n}$
converges to $\Pi$, $\sigma_{n}$ converges to the median of
$\|\mathbf{Z}-\mathbf{Z}^{{}^{\prime}}\|$, where
$\mathbf{Z},\mathbf{Z}^{{}^{\prime}}\stackrel{{\scriptstyle
i.i.d}}{{\sim}}\Pi$. Our test remains consistent for such choices of the
bandwidth. This result is stated below.

###### Theorem 4.3

Suppose that $\mathbf{X}^{(1)},\mathbf{X}^{(2)},\ldots,\mathbf{X}^{(n)}$ are
independent copies of $\mathbf{X}$, which follows a multivariate distribution
with continuous univariate marginals and $\mathsf{C}_{\mathbf{X}}\neq\Pi$.
Also consider a sequence of bandwidths $\\{\sigma_{n}:~{}n\geq 1\\}$
converging to some $\sigma_{0}>0$. Then, power of the proposed test based on
${\widehat{I}}_{\sigma_{n},n}(X)$ converges to $1$ as the sample size $n$
diverges to infinity.

In our experiments, we observed that median heuristic performs well when the
relationships among the variables are nearly monotone. But in cases of complex
non-monotone relationships, use of smaller bandwidths often yields better
results. In such cases, instead of median, one can use lower quantiles of
pairwise distances.

Figure 2: Powers of the test for bandwidths based on different quantiles of
pairwise distances.

To demonstrate this, we considered two simple examples involving bivariate
data sets. In one example, observations were generated from the ‘Two
parabolas’-type distribution mentioned in Newton (2009) (see Figure 4(e)) and
in the other example, they were generated from a bivariate normal distribution
with correlation coefficient 0.5. In each case, we generated $25$ observations
and repeated the experiment $10000$ times to estimate the powers of the tests
based on $\widehat{I}_{\sigma,n}(\mathbf{X})$ for different choices of
$\sigma$ based on different quantiles ($0.01,0.02,0.05,0.1,0.2$ and $0.5$) of
pairwise distances. Figure 2 clearly shows that though median of pairwise
distances worked well in the second example, smaller quantiles had better
results in the first. This figure clearly shows that depending on the
underlying distribution of $\mathbf{X}$, sometimes we need to use larger
bandwidth, whereas sometimes smaller bandwidths may perform better. While
larger bandwidths successfully detect global linear or monotone relationships
among the variables, smaller bandwidths are useful for detecting non-monotone
or local patterns. In order to capture both types of dependence, here we adopt
a multi-scale approach, where we look at the results for several choices of
bandwidth and then aggregate them judiciously to come up with the final
decision.

Figures 3(a) and 3(b) show the observed p-values for different choices of the
bandwidth (based on quantiles of pairwise distances) when a sample of size 25
was generated from bivariate normal distributions with correlation coefficient
0 and 0.5, respectively. Clearly, these plots of p-values carry more
information than just the final result. In the first case, higher $p$-values
for all choices of the bandwidth give a visual evidence in favor $H_{0}$,
while smaller $p$-values for a long range of bandwidths in the second case
indicates dependence between the two coordinate variables. Also the pattern of
p-values can reveal the structure of dependence among the variables. For
instance, smaller p-values for larger bandwidths indicate that the
relationship between the two variables is nearly monotone, while those for
smaller bandwidths indicate complex, non-monotone relations.

Figure 3: p-values for bandwidths corresponding to different quantiles of
pairwise distance.

One way of aggregating the results corresponding to $m$ different bandwidths
$\sigma^{(1)},\ldots,\sigma^{(m)}$ is to use
$T_{Sum}=\sum_{i=1}^{m}{\widehat{I}}_{\sigma^{(i)},n}(\mathbf{X})$ or
$T_{Max}=\max_{1\leq i\leq m}{\widehat{I}}_{\sigma^{(i)},n}(\mathbf{X})$ as
the test statistic. Following Sarkar and Ghosh (2018), one can also use
another method based on false discovery rate (FDR). Let $p_{i}$ be the
$p$-value of the test based on $\sigma^{(i)}$ ($i=1,2,\ldots,m$) and
$p_{(1)}\leq p_{(2)}\leq\ldots\leq p_{(m)}$ be the corresponding order
statistics. We reject $H_{0}$ at level $\alpha$ if and only if the set
$\\{i:p_{(i)}<i~{}\alpha/m\\}$ is non-empty. Benjamini and Hochberg (1995)
proposed this method for controlling FDR for a set $m$ independent tests.
Later, Benjamini and Yekutieli (2001) showed that it also controls FDR when
the tests statistics are positively regression dependent. Since we are testing
the same hypothesis for different choices of the bandwidth, this method
controls the level of the test as well (see, e.g., Cuesta-Albertos and
Febrero-Bande, 2010). It is difficult to prove positive regression dependence
among the test statistics corresponding to different choices of bandwidth.
However, all pairwise correlations (computed over 10000 simulations) among
these test statistics were found to be positive in all of our numerical
experiments. This gives an indication of positive regression dependence among
the test statistics and thereby provides an empirical justification for using
the above method. The following theorem shows the large sample consistency of
the multi-scale versions of our tests based on $T_{Sum}$, $T_{Max}$ and $FDR$.

###### Theorem 4.4

Suppose that $\mathbf{X}^{(1)},\mathbf{X}^{(2)},\ldots,\mathbf{X}^{(n)}$ are
independent copies of $\mathbf{X}$ following a multivariate distribution
having continuous univariate marginals and $\mathsf{C}_{\mathbf{X}}\neq\Pi$.
Then, the powers of the multi-scale versions of the proposed test based on
$T_{Sum}$, $T_{Max}$ and FDR converge to $1$ as the sample size tends to
infinity.

##  5 Results from the analysis of simulated and real data sets

We analyzed several simulated and real data sets to compare the performance of
our proposed tests with some popular tests available in the literature. In
particular, we considered the dHSIC test (Pfister et al., 2017), the mdCov
test (Fan et al., 2017) and the tests based on multivariate extensions of
Hoeffding’s $\phi$ (Gaißer et al., 2010) and Spearman’s $\rho$ (Nelsen, 1996)
statistics for comparison. For the implementation of the dHSIC test, we used
the R package “dHSIC” (Pfister and Peters, 2016), where we used the Gaussian
kernel with the default bandwidth chosen based on median heuristic. For the
mdCov test, we used the codes provided by the authors. Following their
suggestion (see Fan et al., 2017, p. 198), we standardized the data set and
used unit bandwidth for all experiments. We used different options for the
weight function available in the codes and reported the best result. For the
tests based on Hoeffding’s $\phi$ and Spearman’s $\rho$ statistics (henceforth
referred to as the Hoeffding test and the Spearmen test, respectively), we
used our own codes. For all these methods, conditional tests based on 10000
random permutations were used. We also considered the tests proposed in Póczos
et al. (2012), where they suggested to compute the cut-offs based on
probability inequalities. But this choice of cut-off makes the resulting tests
very conservative. As a result, they had much lower powers compared to all
other tests considered here. So, we decided not to report those results in
this article. For our proposed tests, we started with the bandwidth based on
median heuristic ($\sigma_{0.5}$, say) and considered other bandwidths of the
form $(0.5)^{i}\times\sigma_{0.5}$ for $i=1,2,\ldots,m$, where
$m=\lceil\log_{2}(\sigma_{0.5}/\sigma_{0.01})\rceil$, for $\sigma_{0.01}$
being the bandwidth based on the first quantile. Results for these bandwidths
were aggregated using the three methods discussed in Section 4. However,
overall performance of the tests based on $T_{Max}$ and FDR was superior than
the test based on $T_{Sum}$. So, here we report the results for the tests
based on $T_{Max}$ and FDR only. Throughout this article, all tests are
considered to have $5\%$ nominal level.

###  5.1 Analysis of simulated data sets

We begin with eight simulated examples involving bivariate observations.
Scatter plots of these data sets are displayed in Figure 4. For each example,
we repeated our experiment 10000 times, and the power of a test was estimated
by proportion of times it rejected $H_{0}$. These estimated powers of
different tests are reported in Figure 5. The first six examples (see Figures
4(a)-4(f)) are taken from Newton (2009), who considered six unusual bivariate
distributions. In all these examples, $X_{1}$ and $X_{2}$ are uncorrelated. In
‘four clouds’ data, they are independent as well. In this example, almost all
tests had powers close to the nomination level of 0.05 (see Figure 5(a)). Only
the test based on FDR had slightly low powers, which is quite expected in view
of the conservative nature of the tests based on FDR.

Figure 4: Different bivariate distributions used in simulation study.

In the next five examples, $X_{1}$ and $X_{2}$ are not independent. In the
example with ‘W’ type data, our proposed test based on FDR had the best
overall performance followed by the test based on $T_{Max}$ and the dHSIC test
(see Figure 5(b)). Powers of all other tests were much lower. Spearman and
Hoeffding tests couldn’t reject $H_{0}$ even on a single occasion. These two
tests had zero power in ‘Circle’-type data as well (see Figure 5(f)). In that
example, the dCov test also had zero power, and the performance of the dHSIC
test was not satisfactory as well. But our proposed test based on $T_{Max}$
and FDR performed well. These two tests outperformed their competitors in ‘Two
parabolas’-type data as well (see Figure 5(e)). In that example, Spearman and
Hoeffding tests again performed poorly, but the performance of mdCov and dHSIC
tests was somewhat better. In the ‘Parabola’-type data, the test based on FDR,
the dHSIC test and the mdCov test had higher powers than all other tests
considered here (see Figure 5(d)). Among the rest, the test based on $T_{Max}$
had better performance. Only in the case of ‘Diamond’-type data, dHSIC and
mdCov tests outperformed our proposed methods (see Figure 5(c)). However, even
in this example, our proposed tests performed well. They had much higher
powers compared to Spearman and Hoeffding tests.

Figure 5: Powers of different tests in simulated bivariate data sets.

Unlike the previous six examples, in our last two examples, $X_{1}$ and
$X_{2}$ are positively correlated. In the example with ‘hyperplane’-type data
(see Figure 3g), we have $X_{1}=U$ and $X_{2}=U+V$. where
$U,V\stackrel{{\scriptstyle i.i.d.}}{{\sim}}\text{U}(-1,1)$. In the example
with ‘normal’ data (see Figure 3h), $(X_{1},X_{2})$ follows a bivariate normal
distribution with correlation coefficient 0.4. In these two examples,
Hoeffding and Spearman tests had the best performance closely followed by the
mdCov test. Our proposed tests and the dHSIC test also had competitive
performance. Among these three tests, the test based on FDR had an edge.

Figure 6: Powers of different tests in eight-dimensional simulated data sets.

Next we carried out our experiments with some eight dimensional data sets,
which can be viewed as multivariate extensions of the bivariate data sets
considered above. For each of the first six examples, we generated two
independent observations from the bivariate distribution, and then four
independent $N(0,1)$ variables were augmented to it to get a vector of
dimension eight. For the ‘hyperplane’-type data, we generated seven i.i.d.
$N(0,1)$ variables $X_{2},X_{3},\cdots,X_{8}$, and then define
$X_{1}=(X_{2}+\ldots+X_{8})+\epsilon$, where $\epsilon\sim N(0,1)$. For the
example with ‘normal’ data, $\mathbf{X}$ was generated from a $8$-dimensional
normal distribution with the mean vector ${\bf 0}$ and the dispersion matrix
$\mathbf{\Sigma}=((a_{i,j}))$, where
$a_{i,j}=0.4^{|i-j|}~{}\forall~{}i,j=1,2,\ldots,8$.

In the example with ‘four clouds’ data, again the test best on FDR had powers
slightly lower than 0.05, but those of all other tests were close to the
nominal level (see Figure 6(a)). Figure 6 clearly shows that except for
‘diamond’-type data, in all other cases, our tests based on $T_{Max}$ and FDR
had best overall performance among the tests considered here. Note that the
dHSIC test needs the sample size to be at least twice the dimension of the
data (i.e., twice the number of coordinate variables) for its implementation.
So, it could not be used in some cases. In such cases, we considered its power
to be zero.

Next, we consider two interesting examples, where none of the lower
dimensional marginals have dependency among the coordinate variables. In
Example-A, we generate four independent $U(-1,1)$ variables
$U_{1},\ldots,U_{4}$, and if their product is positive, we define
$X_{i}=U_{i}$ for $i=1,2,\ldots,4$. In Example-B, we generate
$U_{1},\ldots,U_{4}$ independently from $N(0,1)$ to define
$X_{i}=U_{i}~{}sign(U_{i+1})$ for $i=1,2,3$ and $X_{4}=U_{4}~{}sign(U_{1})$.
In both of these examples, we carried out our experiments 10000 times as
before to compute the powers of different tests. Note that tests based on any
dependency measure, which is not irreducible, will fail to detect the
dependency among the coordinate variables in these examples. Our proposed
methods, particularly the test based on $T_{Max}$ had excellent performance in
these two data sets (see Figure 7). In Example-A, the dHSIC test also had
competitive powers, but performances of all other tests were much inferior.

Figure 7: Powers in different tests in four-dimensional examples.

###  5.2 Analysis of Combined Cycle Power Plant Data

We also analyzed a real data set, namely the Combined cycle power plant data,
for further evaluation of our proposed methods. This data set is available at
the UCI Machine Learning Repository https://archive.ics.uci.edu/ml/datasets/.
It contains 9568 observations from a Combined Cycle Power Plant over a period
of six years (2006-2011), when the plant was set to work with full load. Each
observation consists of hourly average values of ambient temperature, ambient
pressure, relative humidity, exhaust vacuum and electric energy output. The
idea was to predict electric energy output, which is dependent on other
variables. When we used different methods to test for the independence among
these five variables, all tests rejected $H_{0}$ on almost all occasions even
when they were used on random subsets of size $10$ drawn from the data set.
So, next we removed electric energy output from our analysis and carried out
tests for independence among the other four variables.

When we used the whole data set (ignoring electric energy output) for testing,
all tests rejected $H_{0}$. It gives us a clear indication that these four
variables have significant dependence among themselves, and different tests
can be compared based on their powers. But based on a single experiment with
the whole data set, it was not possible to compare among the powers of
different test procedures. So, following the idea of Sarkar and Ghosh (2018),
we carried out our experiments with subsets of different sizes. For each
subset size (i.e., sample size), the experiment was repeated 10000 times to
estimate the powers of different tests by proportion of times they rejected
$H_{0}$. These estimated powers for different tests are shown in Figure 8.
This figure clearly shows that in this example, our proposed test based on FDR
outperformed its all competitors. The test based on $T_{Max}$ also performed
well. Among the rest, only the dHSIC test had satisfactory performance.

Figure 8: Powers of different tests in real datasets.

##  6 Discussion

In this article, we have proposed a copula based multivariate dependency
measure and established some of its theoretical properties. Unlike many other
existing copula based measures, our dependency measure is invariant under
strictly monotone transformations of the coordinate variables. Interestingly,
it takes the value $0$ when the coordinate variables are independent and takes
the value $1$ when for each pair of the coordinate variables, one is a strict
monotone function of the other. A data based estimate of this measure is
proposed and some distribution-free tests of independence are constructed
based on this estimate. Some nice theoretical properties of this estimate have
also been derived and the large sample consistency of the resulting tests has
been proved under appropriate regularity conditions. Unlike the dHSIC test,
our proposed tests can be used even when the sample size is smaller than the
number of variables. However, our proposed methods are not above all
limitations. These rank based methods are mainly applicable when the
coordinate variables are continuous in nature. In the case of discrete data,
one may need to resolve the ties arbitrarily to define the ranks. The choice
of the bandwidth is another issue to be resolved. In this article, we have
adopted a multi-scale approach, where the results for different bandwidths are
aggregated judiciously. The resulting tests worked well in all simulated and
real data sets considered in this article. But instead of taking such a multi-
scale approach, if we can choose a suitable data driven estimate of the
bandwidth, that can further improve the performance of our methods, both in
terms of power and computing time. This can be considered as a problem for
future research.

###### Acknowledgements.

We would like to thank the Associate Editor and two anonymous reviewers for
carefully reading the earlier version of the manuscript and also for providing
us with several helpful comments. We are also thankful to Dr. Pierre Lafaye de
Micheaux for providing the codes of the mdCov test.

## Appendix

######  Proof (Proof of Lemma 1)

For any permutation $\tau$ on $\mathbb{R}^{d}$, we have
$k_{\sigma}(\tau(\mathbf{x}),\tau(\mathbf{y}))$
$=k_{\sigma}(\mathbf{x},\mathbf{y})$ and also, $\mathbf{T}\sim\Pi$ implies
$\tau(\mathbf{T})\sim\Pi$. Using these, one gets

$\displaystyle\mathrm{E}_{(\mathbf{S},\mathbf{S}^{{}^{\prime}})\sim\mathsf{C}_{\tau(\mathbf{X})}\otimes\mathsf{C}_{\tau(\mathbf{X})}}[k_{\sigma}(\mathbf{S},\mathbf{S}^{{}^{\prime}})]=\mathrm{E}_{(\mathbf{S},\mathbf{S}^{{}^{\prime}})\sim\mathsf{C}_{\mathbf{X}}\otimes\mathsf{C}_{\mathbf{X}}}[k_{\sigma}(\tau(\mathbf{S}),\tau(\mathbf{S}^{{}^{\prime}}))]$
$\displaystyle=\mathrm{E}_{(\mathbf{S},\mathbf{S}^{{}^{\prime}})\sim\mathsf{C}_{\mathbf{X}}\otimes\mathsf{C}_{\mathbf{X}}}[k_{\sigma}(\mathbf{S},\mathbf{S}^{{}^{\prime}})]\;\mbox{
and}$
$\displaystyle\mathrm{E}_{(\mathbf{S},\mathbf{T})\sim\mathsf{C}_{\tau(\mathbf{X})}\otimes\Pi}[k_{\sigma}(\mathbf{S},\mathbf{T})]=\mathrm{E}_{(\mathbf{S},\mathbf{T})\sim\mathsf{C}_{\mathbf{X}}\otimes\Pi}[k_{\sigma}(\tau(\mathbf{S}),\tau(\mathbf{T}))]$
$\displaystyle=\mathrm{E}_{(\mathbf{S},\mathbf{T})\sim\mathsf{C}_{\mathbf{X}}\otimes\Pi}[k_{\sigma}(\mathbf{S},\mathbf{T})].$

It follows that
$\gamma_{k_{\sigma}}(\mathsf{C}_{\tau(\mathbf{X})},\Pi)=\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{X}},\Pi)$
and hence $I_{\sigma}(\tau(\mathbf{X}))=I_{\sigma}(\mathbf{X})$.

Next, let $g:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ be a function of the
form $g(x_{1},x_{2}\cdots,x_{d})=(g_{1}(x_{1})$,
$g_{2}(x_{2}),\cdots,g_{d}(x_{d}))$, where $g_{i_{1}},\cdots,g_{i_{s}}$ are
strictly increasing and $g_{j_{1}},\cdots,g_{j_{t}}$ are strictly decreasing
with $s+t=d$. Consider the function
$f:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ given by
$f(x_{1},x_{2},\cdots,x_{d})=(f_{1}(x_{1}),f_{2}(x_{2}),\cdots,f_{d}(x_{d}))$,
where $f_{i_{l}}(x)=x\;\forall\;l=1,2,\cdots,s$ and
$f_{j_{l}}(x)=1-x\;\forall\;l=1,2,\cdots,t$. It can be easily verified that if
$\mathbf{S}\sim\mathsf{C_{g(\mathbf{X})}}$ then
$f(\mathbf{S})\sim\mathsf{C_{\mathbf{X}}}$. Applying this and the fact that
$k_{\sigma}(\mathbf{S},\mathbf{S}^{{}^{\prime}})=k_{\sigma}(f(\mathbf{S}),f(\mathbf{S}^{{}^{\prime}}))$,
we get

$\displaystyle\mathrm{E}_{(\mathbf{S},\mathbf{S}^{{}^{\prime}})\sim\mathsf{C}_{g(\mathbf{X})}\otimes\mathsf{C}_{g(\mathbf{X})}}[k_{\sigma}(\mathbf{S},\mathbf{S}^{{}^{\prime}})]=\mathrm{E}_{(\mathbf{S},\mathbf{S}^{{}^{\prime}})\sim\mathsf{C}_{g(\mathbf{X})}\otimes\mathsf{C}_{g(\mathbf{X}})}[k_{\sigma}(f(\mathbf{S}),f(\mathbf{S}^{{}^{\prime}}))]$
$\displaystyle=\mathrm{E}_{(\mathbf{S},\mathbf{S}^{{}^{\prime}})\sim\mathsf{C}_{\mathbf{X}}\otimes\mathsf{C}_{\mathbf{X}}}[k_{\sigma}(\mathbf{S},\mathbf{S}^{{}^{\prime}})].$

By similar argument and using the fact that $\mathbf{T}\sim\Pi$ implies
$f(\mathbf{T})\sim\Pi$, one gets

$\displaystyle\mathrm{E}_{(\mathbf{S},\mathbf{T})\sim\mathsf{C}_{g(\mathbf{X})\otimes\Pi}}[k_{\sigma}(\mathbf{S},\mathbf{T})]=\mathrm{E}_{(\mathbf{S},\mathbf{T})\sim\mathsf{C}_{g(\mathbf{X})}\otimes\Pi}[k_{\sigma}(f(\mathbf{S}),f(\mathbf{T}))]$
$\displaystyle=\mathrm{E}_{(\mathbf{S},\mathbf{T})\sim\mathsf{C}_{\mathbf{X}}\otimes\Pi}[k_{\sigma}(\mathbf{S},\mathbf{T})].$

Thus
$\gamma_{k_{\sigma}}(\mathsf{C}_{g(\mathbf{X})},\Pi)=\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{X}},\Pi)$,
whence, $I_{\sigma}(g(\mathbf{X}))=I_{\sigma}(\mathbf{X})$, proving the
invariance of $I_{\sigma}(\mathbf{X})$ under strictly monotonic
transformations of $X_{1},X_{2},\ldots,X_{d}$. ∎

######  Proof (Proof of Lemma 2)

Let $\mathbf{X}$ be a random vector with continuous marginals, for which there
is a $j$ such that each $X_{i},\,i\not=j$ is a strictly monotonic function of
$X_{j}$. Then, by Lemma 1, we have
$I_{\sigma}(\mathbf{X})=I_{\sigma}(\mathbf{Y})$ where
$\mathbf{Y}=(X_{j},X_{j},\cdots,X_{j})$. But then $\mathsf{C}_{\mathbf{Y}}$ is
the maximum copula $\mathsf{M},\,$ so that, by definition,
$I_{\sigma}(\mathbf{Y})=1$. ∎

######  Proof (Proof of Theorem 1)

This proof has two steps. At the first step, we prove that
$C_{\sigma,d}=\gamma^{2}_{k_{\sigma}}(\mathsf{M},\Pi)$. At the second step, we
prove that
$\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{X}},\Pi)=\left(\frac{\sigma}{\sqrt{2\pi}}\right)^{d}\int_{\mathbb{R}^{d}}|\varphi_{\mathsf{C}_{\mathbf{X}}}(\mathbf{w})-\varphi_{\Pi}(\mathbf{w})|^{2}\exp\left(-\frac{\sigma^{2}}{2}\mathbf{w}^{\intercal}\mathbf{w}\right)\,d\mathbf{w}$.
Clearly, proving these two steps will complete the proof.

First step: Note that for
$(\mathbf{S},\mathbf{S}^{{}^{\prime}}\\!\\!,\mathbf{T},\mathbf{T}^{{}^{\prime}})\sim\mathsf{M}\otimes\mathsf{M}\otimes\Pi\otimes\Pi$,
we have

$\displaystyle\gamma^{2}_{k_{\sigma}}(\mathsf{M},\Pi)=\mathrm{E}[k_{\sigma}(\mathbf{S},\mathbf{S}^{{}^{\prime}})]-2\mathrm{E}[k_{\sigma}(\mathbf{S},\mathbf{T})]+\mathrm{E}[k_{\sigma}(\mathbf{T},\mathbf{T}^{{}^{\prime}})]$
$\displaystyle=\int^{1}_{0}\\!\\!\\!\int^{1}_{0}\\!\\!e^{-\frac{d(u-v)^{2}}{2\sigma^{2}}}\,du\,dv-2\int^{1}_{0}\\!\\!\left[\int^{1}_{0}\\!\\!e^{-\frac{(u-v)^{2}}{2\sigma^{2}}}\,du\right]^{d}\,dv+\left[\int^{1}_{0}\\!\\!\\!\int^{1}_{0}e^{-\frac{(u-v)^{2}}{2\sigma^{2}}}\,du\,dv\right]^{d}$
$\displaystyle=\kappa\left(\frac{\sigma}{\sqrt{d}}\right)-2\int^{1}_{0}\lambda^{d}(u,\sigma)\,du+\kappa^{d}(\sigma)=C_{\sigma,d}.$

Second step: We use the well-known formula for Fourier transform of the
$d$-dimensional Gaussian density:

$\exp\left(-\frac{1}{2\sigma^{2}}\mathbf{x}^{\intercal}\mathbf{x}\right)=\int_{\mathbb{R}^{d}}e^{-\sqrt{-1}\mathbf{x}^{\intercal}\mathbf{w}}\cdot\left(\frac{\sigma}{\sqrt{2\pi}}\right)^{d}\exp\left(-\frac{\sigma^{2}}{2}\mathbf{w}^{\intercal}\mathbf{w}\right)\,d\mathbf{w},\;\mathbf{x}\in\mathbb{R}^{d}.$

This gives us

$k_{\sigma}(\mathbf{x},\mathbf{y})=\left(\frac{\sigma}{\sqrt{2\pi}}\right)^{d}\int_{\mathbb{R}^{d}}e^{-\sqrt{-1}\mathbf{x}^{\intercal}\mathbf{w}}\cdot
e^{\sqrt{-1}\mathbf{y}^{\intercal}\mathbf{w}}\exp\left(-\frac{\sigma^{2}}{2}\mathbf{w}^{\intercal}\mathbf{w}\right)\,d\mathbf{w},\;\mathbf{x},\mathbf{y}\in\mathbb{R}^{d}.$

Using the representation of $\gamma^{2}_{k}$ from equation (2) and Fubini’s
theorem, one gets

$\displaystyle\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{X}},\Pi)$
$\displaystyle=\\!\left(\frac{\sigma}{\sqrt{2\pi}}\right)^{\\!d}\\!\\!\int_{\mathbb{R}^{d}}\Big{[}\,\varphi_{\mathsf{C}_{\mathbf{X}}}(\mathbf{w})\overline{\varphi_{\mathsf{C}_{\mathbf{X}}}}(\mathbf{w})+\varphi_{\Pi}(\mathbf{w})\overline{\varphi_{\Pi}}(\mathbf{w})$
$\displaystyle\qquad\qquad\qquad-2\varphi_{\mathsf{C}_{\mathbf{X}}}(\mathbf{w})\overline{\varphi_{\Pi}}(\mathbf{w})\,\Big{]}\,\exp\left(-\frac{\sigma^{2}}{2}\mathbf{w}^{\intercal}\mathbf{w}\right)\,d\mathbf{w},$

from which the second part follows. ∎

###### Lemma L1

Let $(\mathbf{X},\mathbf{Y})$ and
$(\mathbf{X}^{{}^{\prime}},\mathbf{Y}^{{}^{\prime}})$ be independent and
identically distributed random vectors taking values in
$\mathcal{\mathbf{X}}\times\mathcal{\mathbf{Y}}$. Given symmetric measurable
functions $k:\mathcal{X}\times\mathcal{X}\to\mathbb{R}$ and
$\overline{k}:\mathcal{Y}\times\mathcal{Y}\to\mathbb{R}$, define

$\displaystyle
V=k(\mathbf{X},\mathbf{X}^{{}^{\prime}})-\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\Big{|}\mathbf{X}\right]-\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\Big{|}\mathbf{X}^{{}^{\prime}}\right]+\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\right]$
$\displaystyle
W=\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})-\mathrm{E}\left[\,\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\Big{|}\mathbf{Y}\right]-\mathrm{E}\left[\,\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\Big{|}\mathbf{Y}^{{}^{\prime}}\right]+\mathrm{E}\left[\,\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\right].$

Then, we have

$\displaystyle\mathrm{E}\left[\,VW\,\right]=\,\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\,\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\right]+\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\right]\mathrm{E}\left[\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\right]$
$\displaystyle-2\,\mathrm{E}\left[\,\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\Big{|}\mathbf{X}\right]\mathrm{E}\left[\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\Big{|}\mathbf{Y}\right]\,\right].$

###### Proof

The proof is based on expanding the product $VW$ and then taking term-by-term
expectations. One and only one term gives
$\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\,\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\right]$.
The seven terms, where at least one of
$\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\right]$ or
$\mathrm{E}\left[\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\right]$
appear as a factor, and the two terms
$\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\Big{|}\mathbf{X}\right]\cdot\mathrm{E}\left[\,\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\Big{|}\mathbf{Y}^{{}^{\prime}}\right]$
and
$\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\Big{|}\mathbf{X}^{{}^{\prime}}\right]\cdot\mathrm{E}\left[\,\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\Big{|}\mathbf{Y}\right]$,
will all give the same expectation
$\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\right]\mathrm{E}\left[\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\right]$
(the last two because of independence of $(\mathbf{X},\mathbf{Y})$ and
$(\mathbf{X}^{{}^{\prime}},\mathbf{Y}^{{}^{\prime}})$). Taking into account
the signs of these nine terms with the same expectation, we would be left with
just one with a positive sign. Next, the remaining six terms will all have the
same expectation, namely,
$\mathrm{E}\left[\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\Big{|}\mathbf{X}\right]\mathrm{E}\left[\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\Big{|}\mathbf{Y}\right]\right]$.
For two of the terms, this is straightforward. But the other four terms need
judicious use of properties of conditional expectation. For example, by
independence of $(\mathbf{X},\mathbf{Y})$ and
$(\mathbf{X}^{{}^{\prime}},\mathbf{Y}^{{}^{\prime}})$, we have
$\mathrm{E}\left[\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\Big{|}\mathbf{Y}\right]=\mathrm{E}\left[\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\Big{|}(\mathbf{X},\mathbf{Y})\right]$
and similarly
$\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\Big{|}(\mathbf{X},\mathbf{Y})\right]=\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\Big{|}\mathbf{X}\right]$.
Using these, we get

$\displaystyle\mathrm{E}\left[\,k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\mathrm{E}\left[\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\Big{|}\mathbf{Y}\right]\,\right]$
$\displaystyle=\mathrm{E}\left[\,k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\mathrm{E}\left[\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\Big{|}(\mathbf{X},\mathbf{Y})\right]\,\right]$
$\displaystyle=\mathrm{E}\left[\,\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\big{|}(\mathbf{X},\mathbf{Y})\right]\mathrm{E}\left[\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\Big{|}(\mathbf{X},\mathbf{Y})\right]\,\right]$
$\displaystyle=\mathrm{E}\left[\,\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\big{|}\mathbf{X}\right]\mathrm{E}\left[\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\Big{|}\mathbf{Y}\right]\,\right].$

One can similarly handle other three terms. Considering the signs of these six
terms with the same expectation, one is left with
$-2\mathrm{E}\left[\mathrm{E}\left[k(\mathbf{X},\mathbf{X}^{{}^{\prime}})\big{|}\mathbf{X}\right]\mathrm{E}\left[\overline{k}(\mathbf{Y},\mathbf{Y}^{{}^{\prime}})\Big{|}\mathbf{Y}\right]\right]$.
This completes the proof. ∎

######  Proof (Proof of Theorem 2)

(a) By definition, $V_{1}$ and $V_{2}$ have zero means. Using Lemma L1 with
kernels $k=\overline{k}=k_{\sigma}$ on
$\mathcal{X}=\mathcal{Y}=\mathbb{R}^{d}$, we get

$\displaystyle\mathrm{Cov}[V_{1},V_{2}]=$
$\displaystyle\mathrm{E}[V_{1}V_{2}]=\mathrm{E}\left[k_{\sigma}(T_{1},T_{1}^{{}^{\prime}})k_{\sigma}(T_{2},T_{2}^{{}^{\prime}})\right]+\mathrm{E}\left[k_{\sigma}(T_{1},T_{1}^{{}^{\prime}})\right]\mathrm{E}\left[k_{\sigma}(T_{2},T_{2}^{{}^{\prime}})\right]$
$\displaystyle-2\,\mathrm{E}\left[\mathrm{E}\left[k_{\sigma}(T_{1},T_{1}^{{}^{\prime}})\Big{|}T_{1}\right]\mathrm{E}\left[k_{\sigma}(T_{2},T_{2}^{{}^{\prime}})\Big{|}T_{2}\right]\right]=\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{(X,Y)},\Pi).$

One can similarly show that
$\mathrm{Var}[V_{1}]=\mathrm{Var}[V_{2}]=\gamma^{2}_{k_{\sigma}}(\mathsf{M},\Pi)$
and hence $I^{2}_{\sigma}(X,Y)=\mathrm{Cor}[V_{1},V_{2}]$. The inequality
$I_{\sigma}(X_{1},X_{2})\leq 1$ follows from it. Further, from the condition
for equality in Cauchy-Schwartz inequality and the fact that $V_{1}$ and
$V_{2}$ are identically distributed, it follows that
$I_{\sigma}(X_{1},X_{2})=1$ if and only if $V_{1}=V_{2}$ almost surely.

Since $T_{1}$ and $T_{1}^{{}^{\prime}}$ are independent and uniformly
distributed random variables on $[0,1]$ and so also are $T_{2}$ and
$T_{2}^{{}^{\prime}}$, it follows that $V_{1}=V_{2}$ almost surely if and only
if $g(T_{1},T_{1}^{{}^{\prime}})=g(T_{2},T_{2}^{{}^{\prime}})$ almost surely,
where $g(x,y)=k_{\sigma}(x,y)-\lambda(x,\sigma)-\lambda(y,\sigma)$ with
$\lambda(\cdot,\sigma)$ as defined in Theorem 1.

Now, using the facts that $(T_{1},T_{2})$ and
$(T_{1}^{{}^{\prime}},T_{2}^{{}^{\prime}})$ are independent and identically
distributed with values in $[0,1]^{2}$ and that the function $g$ is uniformly
continuous on the compact set $[0,1]^{2}$, one can easily deduce that
$g(T_{1},T_{1}^{{}^{\prime}})=g(T_{2},T_{2}^{{}^{\prime}})~{}\text{a.s.}$
implies $g(T_{1},T_{1})=g(T_{2},T_{2})$ a.s. But, this, in turn, implies that
$\Phi\left(\frac{T_{1}}{\sigma}\right)+\Phi\left(\frac{1-T_{1}}{\sigma}\right)=\Phi\left(\frac{T_{2}}{\sigma}\right)+\Phi\left(\frac{1-T_{2}}{\sigma}\right)~{}\text{a.s.}$.
From this, we may conclude that $\mathrm{Pr}\left[T_{2}=T_{1}~{}\text{or
}T_{2}=1-T_{1}\right]=1$ and also
$\mathrm{Pr}\left[\lambda(T_{1},\sigma)=\lambda(T_{2},\sigma)\right]=1$. Of
course, the same would be true of the pair
$(T_{1}^{{}^{\prime}},T_{2}^{{}^{\prime}})$, which is moreover independent of
the pair $(T_{1},T_{2})$.

Using these in the equality
$g(T_{1},T_{1}^{{}^{\prime}})=g(T_{2},T_{2}^{{}^{\prime}})~{}\text{a.s.}$, one
obtains
$k_{\sigma}(T_{1},T_{1}^{{}^{\prime}})=k_{\sigma}(T_{2},T_{2}^{{}^{\prime}})~{}\text{a.s.}$,
which implies that
$|T_{1}-T_{1}^{{}^{\prime}}|=|T_{2}-T_{2}^{{}^{\prime}}|~{}\text{a.s.}$ We
conclude that either $T_{2}=T_{1}~{}\text{a.s.}$ or
$T_{2}=1-T_{1}~{}\text{a.s.}$ Thus the copula distribution of $(X_{1},X_{2})$
is either the distribution of $(T_{1},T_{1})$ or that of $(T_{1},1-T_{1})$
where $T_{1}$ is uniformly distributed on $[0,1]$. So, $X_{1}$ and $X_{2}$ are
almost surely strictly monotone functions of each other.

(b) For $|r|<1$, let $\phi_{r}$ denote the density of the standard bivariate
normal distribution with correlation coefficient $r$. Also, let $\Phi$ and
$\phi$ denote respectively the cumulative distribution function and the
density function of the standard univariate normal distribution. It is well-
known that the copula distribution of any bivariate normal distribution with
correlation coefficient $r$ is the same as that of the standard bivariate
normal distribution with the same correlation coefficient. Using the well-
known Mehler’s representation (see Kibble (1945), Page 1) of standard
bivariate normal density with correlation $r$, one then gets that, for
$|r|<1$, the copula distribution $\mathcal{C}(r)$ of any bivariate normal
distribution with correlation coefficient $r$ has density given by

$\frac{\phi_{r}(\Phi^{-1}(u),\Phi^{-1}(v))}{\phi(\Phi^{-1}(u))\phi(\Phi^{-1}(v))}=\sum_{i=0}^{\infty}\frac{r^{i}}{i!}H_{i}((\Phi^{-1}(u))H_{i}(\Phi^{-1}(v)),\quad(u,v)\in[0,1]^{2},\vspace{-5pt}$

where $\\{H_{i}(x),i\geq 0\\}$ are the well-known Hermite polynomials. Using
this, we get that if $(S,T)\sim\mathcal{C}(r_{1})\otimes\mathcal{C}(r_{2})$
with $|r_{1}|<1,|r_{2}|<1$, then

$\displaystyle\mathrm{E}[k_{\sigma}(S,T)]\\!=$
$\displaystyle\\!\\!\\!\int\limits_{[0,1]^{4}}\\!e^{-\frac{(s_{1}-t_{1})^{2}+(s_{2}-t_{2})^{2}}{2\sigma^{2}}}\sum_{i=0}^{\infty}\frac{r_{1}^{i}}{i!}H_{i}(\Phi^{-1}(s_{1}))H_{i}(\Phi^{-1}(s_{2}))$
$\displaystyle\times\sum_{j=0}^{\infty}\\!\frac{r_{2}^{j}}{j!}H_{j}(\Phi^{-1}(t_{1}))H_{j}(\Phi^{-1}(t_{2}))\,ds_{1}\,ds_{2}\,dt_{1}\,dt_{2}\vspace{-10pt}$
(E1)

We now claim that in the above expression, the double summation and
integration can be interchanged. To justify this, we recall that the Hermite
polynomials $\\{H_{i}(\cdot),i\geq 0\\}$ form a complete orthonormal basis for
$L_{2}(\mathbb{R},\phi(x)dx)$ and, in particular, for any $i\geq 0$,
$\int_{[0,1]}\left|H_{i}(\Phi^{-1}(s))\right|\,ds=\int_{\mathbb{R}}\left|H_{i}(x)\right|\phi(x)\,dx\leq\left[\int_{\mathbb{R}}H_{i}^{2}(x)\phi(x)\,dx\right]^{\frac{1}{2}}$
$=1$. As a consequence,

$\displaystyle\sum_{i=0}^{\infty}\sum_{j=0}^{\infty}\int_{[0,1]^{4}}\left|e^{-\frac{(s_{1}-t_{1})^{2}+(s_{2}-t_{2})^{2}}{2\sigma^{2}}}\right|\,\left|\frac{r_{1}^{i}r_{2}^{j}}{i!j!}\right|\,\left|H_{i}(\Phi^{-1}(s_{1}))\right|\,\left|H_{i}(\Phi^{-1}(s_{2}))\right|$
$\displaystyle\qquad\times\left|H_{j}(\Phi^{-1}(t_{1}))\right|\left|H_{j}(\Phi^{-1}(t_{2}))\right|\,ds_{1}\,ds_{2}\,dt_{1}\,dt_{2}$
$\displaystyle\leq\sum_{i=0}^{\infty}\sum_{j=0}^{\infty}\frac{1}{i!j!}\left[\int\limits_{[0,1]}\\!\\!\left|H_{i}(\Phi^{-1}(s))\right|ds\right]^{2}\\!\left[\int\limits_{[0,1]}\\!\\!\left|H_{j}(\Phi^{-1}(t))\right|dt\right]^{2}\\!\\!\leq\\!\sum_{i=0}^{\infty}\sum_{j=0}^{\infty}\frac{1}{i!j!}<\infty.$

We can, therefore, interchange the double summation and integration in the
right-hand-side of the equation (E1) above to obtain that, for any
$r_{1},r_{2}$ with $|r_{1}|<1,|r_{2}|<1$,

$\displaystyle\mathrm{E}_{(S,T)\sim\mathcal{C}(r_{1})\otimes\mathcal{C}(r_{2})}[k_{\sigma}(S,T)]=\sum_{i=0}^{\infty}\sum_{j=0}^{\infty}a_{i,j}r_{1}^{i}r_{2}^{j}\,,$
$\displaystyle{\mbox{w}here}~{}~{}a_{i,j}$
$\displaystyle:=\frac{1}{i!j!}\left[\int_{[0,1]^{2}}e^{-\frac{(u-v)^{2}}{2\sigma^{2}}}H_{i}(\Phi^{-1}(s))H_{j}(\Phi^{-1}(t))\,ds\,dt\right]^{\\!2}$
$\displaystyle=\frac{1}{i!j!}\left[\int_{\mathbb{R}^{2}}e^{-\frac{(\Phi(x)-\Phi(y))^{2}}{2\sigma^{2}}}H_{i}(x)H_{j}(y)\phi(x)\phi(y)\,dx\,dy\right]^{\\!2}.$
(E2)

Observe that $a_{i,j}\geq 0,a_{i,j}=a_{j,i}$ and also,
$(i!j!)a_{i,j}\\!\leq\\!\left[\int\limits_{\mathbb{R}}\left|H_{i}(x)\right|\phi(x)\,dx\right]^{2}\\!\leq\\!1$.
Note that for any bivariate normal random vector $(X_{1},X_{2})$ with
correlation coefficient $r$ (where $|r|<1$), we have
$\gamma^{2}_{k_{\sigma}}(\mathcal{C}_{(X_{1},X_{2})},\Pi)=\gamma^{2}_{k_{\sigma}}(\mathcal{C}(r),\mathcal{C}(0))$,
which equals

$\displaystyle\mathrm{E}_{{}_{(S,S^{{}^{\prime}})\sim\mathcal{C}(r)\otimes\mathcal{C}(r)}}\\![k_{\sigma}(S,\\!S^{{}^{\prime}})]\\!-\\!2\mathrm{E}_{{}_{(S,T)\sim\mathcal{C}(r)\otimes\mathcal{C}(0)}}\\![k_{\sigma}(S,\\!T)]\\!+\\!\mathrm{E}_{{}_{(T,T^{{}^{\prime}})\sim\mathcal{C}(0)\otimes\mathcal{C}(0)}}\\![k_{\sigma}(T,\\!T^{{}^{\prime}})]$
$\displaystyle=\\!\sum_{k=0}^{\infty}\\!\sum_{\begin{subarray}{c}i,j\geq 0\\\
i+j=k\end{subarray}}\\!a_{i,j}\\!r^{k}\\!-\\!2\,\sum_{k=0}^{\infty}a_{k,0}\,r^{k}\\!+\\!a_{0,0}=\\!\sum_{k=1}^{\infty}\\!\sum_{\begin{subarray}{c}i,j\geq
1\\\
i+j=k\end{subarray}}\\!a_{i,j}\;r^{k}\underset{(a)}{=}\,\sum_{k=1}^{\infty}\\!\sum_{\begin{subarray}{c}i,j\geq
1\\\ i+j=2k\end{subarray}}\\!a_{i,j}\;r^{2k}.$

Equality (a) is due to the fact that the $i$th Hermite polynomial $H_{i}$ is
an even or an odd function according as $i$ is even or odd, so that if exactly
one of $i$ and $j$ is odd, then $a_{i,j}=0$, as can easily be seen from
equation (Proof (Proof of Theorem 2)). Therefore, we have

$\displaystyle I^{2}_{\sigma}(X_{1},X_{2})$
$\displaystyle=\gamma^{-2}_{k_{\sigma}}(\mathsf{M},\Pi)\gamma^{2}_{k_{\sigma}}(\mathcal{C}_{(X_{1},X_{2})},\Pi)=\gamma^{-2}_{k_{\sigma}}(\mathsf{M},\Pi)\sum_{k=1}^{\infty}\sum_{\begin{subarray}{c}i,j\geq
1,i+j=2k\end{subarray}}\\!\\!a_{i,j}\;r^{2k}.\vspace{-10pt}$

So, $I^{2}_{\sigma}(X_{1},X_{2})=r^{2}g(r)$, where
$g(r)=\gamma^{-2}_{k_{\sigma}}(\mathsf{M},\Pi)\sum\limits_{k=1}^{\infty}\;\sum\limits_{i,j\geq
1:\,i+j=2k}\\!\\!a_{i,j}\;r^{2(k-1)}$ is a power series in $r^{2}$ with
positive coefficients and hence increasing in $|r|$. So,
$I^{2}_{\sigma}(X_{1},X_{2})=r^{2}\cdot g(r)$ is an increasing function of
$|r|$. ∎

######  Proof (Proof of Theorem 3)

It is enough to show that for every dimension $d~{}(\geq 3)$, there exist two
$d$ dimensional copulas $C_{1}$ and $C_{2}$ with
$\mathcal{M}(C_{1})\neq\mathcal{M}(C_{2})$, such that for any choice of co-
ordinates $\\{i_{1},\cdots,i_{k}\\}\subsetneqq\\{1,\ldots,d\\}$, if
$C_{1}^{{}^{\prime}}$ and $C_{2}^{{}^{\prime}}$ are the associated marginal
copulas arising out of of $C_{1}$ and $C_{2}$, then
$\mathcal{M}(C_{1}^{{}^{\prime}})=\mathcal{M}(C_{2}^{{}^{\prime}})$.

Take $C_{1}$ to be the $d$-dimensional uniform copula $\Pi$. Then
$\mathcal{M}(C_{1})=0$, and also for any lower dimensional marginal copula
$C_{1}^{{}^{\prime}}$ of $C_{1}$, $\mathcal{M}(C_{1}^{{}^{\prime}})=0$. We now
exhibit a $d$-dimensional copula $C_{2}\neq\Pi$ such that any lower
dimensional marginal copula $C_{2}^{{}^{\prime}}$ of $C_{2}$ is uniform
copula. We would then have $\mathcal{M}(C_{1})=0\neq\mathcal{M}(C_{2})$ but
$\mathcal{M}(C_{1}^{{}^{\prime}})=\mathcal{M}(C_{2}^{{}^{\prime}})=0$, which
will complete the proof.

We take $C_{2}$ to be the copula given by the copula density $\mathcal{C}_{2}$
defined as

$\mathcal{C}_{2}(u_{1},u_{2},,\cdots,u_{d})=2\,\mathbb{I}\left[\left(u_{1}-\frac{1}{2}\right)\left(u_{2}-\frac{1}{2}\right)\cdots\left(u_{d}-\frac{1}{2}\right)\geq
0\right],$

where $\mathbb{I}$ denotes the indicator function. To show that all lower
dimensional marginal copulas of $C_{2}$ are uniform, it is enough to show that
the marginal copula $C_{2}^{{}^{\prime}}$ that we get from $C_{2}$ discarding
the $d^{\text{th}}$ co-ordinate, is uniform. Now, note that the density of
$C_{2}^{{}^{\prime}}$ is given by

$\displaystyle\mathcal{C}^{{}^{\prime}}_{2}(u_{1},u_{2},,\cdots,u_{d-1})$
$\displaystyle=\int_{0}^{1}2\mathbb{I}\left[\left(u_{1}-\frac{1}{2}\right)\left(u_{2}-\frac{1}{2}\right)\cdots\left(u_{d}-\frac{1}{2}\right)\geq
0\right]\,du_{d}$
$\displaystyle=\int_{0}^{\frac{1}{2}}2\mathbb{I}\left[\left(u_{1}-\frac{1}{2}\right)\left(u_{2}-\frac{1}{2}\right)\cdots\left(u_{d-1}-\frac{1}{2}\right)\leq
0\right]\,du_{d}$
$\displaystyle\quad+\int_{\frac{1}{2}}^{1}2\mathbb{I}\left[\left(u_{1}-\frac{1}{2}\right)\left(u_{2}-\frac{1}{2}\right)\cdots\left(u_{d-1}-\frac{1}{2}\right)\geq
0\right]\,du_{d}$
$\displaystyle=\mathbb{I}\left[\left(u_{1}-\frac{1}{2}\right)\left(u_{2}-\frac{1}{2}\right)\cdots\left(u_{d-1}-\frac{1}{2}\right)\leq
0\right]$
$\displaystyle\quad+\mathbb{I}\left[\left(u_{1}-\frac{1}{2}\right)\left(u_{2}-\frac{1}{2}\right)\cdots\left(u_{d-1}-\frac{1}{2}\right)\geq
0\right]=1.~{}~{}~{}~{}\hfill\qed$

######  Proof (Proof of Lemma 4)

We shall prove that as $\sigma\to\infty$,
$\sigma^{4}\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{X}},\Pi)\to\sum_{1\leq
i<j\leq d}\mathrm{Cov}^{2}(S_{i},S_{j})$. It will imply that
$\sigma^{4}\gamma^{2}_{k_{\sigma}}(\mathsf{M},\Pi)\to{d\choose
2}\mathrm{Var}^{2}(S_{1})$ as $\sigma\to\infty$, which in turn will imply that
$I^{2}_{\sigma}(\mathbf{X})=\frac{\sigma^{4}\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{X}},\Pi)}{\sigma^{4}\gamma^{2}_{k_{\sigma}}(\mathsf{M},\Pi)}\to\frac{1}{{d\choose
2}}\sum_{1\leq i<j\leq
d}\frac{\mathrm{Cov}^{2}(S_{i},S_{j})}{\mathrm{Var}^{2}(S_{1})}=\frac{1}{{d\choose
2}}\sum_{1\leq i<j\leq
d}\frac{\mathrm{Cov}^{2}(S_{i},S_{j})}{\mathrm{Var}(S_{i})\mathrm{Var}(S_{j})}=\frac{1}{{d\choose
2}}\sum_{1\leq i<j\leq d}\mathrm{Cor}^{2}(S_{i},S_{j})$ as $\sigma\to\infty$;
which is our desired result.

Observe that
$\mathrm{E}k_{\sigma}(\mathbf{S},\mathbf{T})=1-\frac{1}{2\sigma^{2}}\mathrm{E}\|\mathbf{S}-\mathbf{T}\|^{2}_{2}+\frac{1}{8\sigma^{4}}\mathrm{E}\|\mathbf{S}-\mathbf{T}\|^{4}_{2}+\mathcal{O}\left(\frac{1}{\sigma^{6}}\right).$
Assume that $\mathbf{S},\mathbf{S}^{{}^{\prime}},\mathbf{T}$ and
$\mathbf{T}^{{}^{\prime}}$ are four random vectors such that
$(\mathbf{S},\mathbf{S}^{{}^{\prime}},\mathbf{T},\mathbf{T}^{{}^{\prime}})\sim\mathsf{C}_{\mathbf{X}}\otimes\mathsf{C}_{\mathbf{X}}\otimes\Pi\otimes\Pi$.
Then,

$\displaystyle\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{X}},\Pi)$
$\displaystyle=\mathrm{E}k_{\sigma}(\mathbf{S},\mathbf{S}^{{}^{\prime}})-2\mathrm{E}k_{\sigma}(\mathbf{S},\mathbf{T})+\mathrm{E}k_{\sigma}(\mathbf{T},\mathbf{T}^{{}^{\prime}})$
$\displaystyle=-\frac{1}{2\sigma^{2}}\mathrm{E}\left[\|\mathbf{S}-\mathbf{S}^{{}^{\prime}}\|_{2}^{2}+\|\mathbf{T}-\mathbf{T}^{{}^{\prime}}\|_{2}^{2}-2\|\mathbf{S}-\mathbf{T}\|_{2}^{2}\right]$
$\displaystyle\quad+\frac{1}{8\sigma^{4}}\mathrm{E}\left[\|\mathbf{S}-\mathbf{S}^{{}^{\prime}}\|_{2}^{4}+\|\mathbf{T}-\mathbf{T}^{{}^{\prime}}\|_{2}^{4}-2\|\mathbf{S}-\mathbf{T}\|_{2}^{4}\right]+\mathcal{O}\left(\frac{1}{\sigma^{6}}\right).$
Now,
$\displaystyle\mathrm{E}\left[\|\mathbf{S}-\mathbf{S}^{{}^{\prime}}\|_{2}^{2}+\|\mathbf{T}-\mathbf{T}^{{}^{\prime}}\|_{2}^{2}-2\|\mathbf{S}-\mathbf{T}\|_{2}^{2}\right]$
$\displaystyle\quad\quad\quad\quad\quad\quad\quad=\mathrm{E}\sum_{i=1}^{d}\left[(S_{i}-S_{i}^{{}^{\prime}})^{2}+(T_{i}-T_{i}^{{}^{\prime}})^{2}-2(S_{i}-T_{i})^{2}\right]=0,\mbox{
and }$
$\displaystyle\mathrm{E}\left[\|\mathbf{S}-\mathbf{S}^{{}^{\prime}}\|_{2}^{4}+\|\mathbf{T}-\mathbf{T}^{{}^{\prime}}\|_{2}^{4}-2\|\mathbf{S}-\mathbf{T}\|_{2}^{4}\right]$
$\displaystyle=\mathrm{E}\sum_{i=1}^{d}\left[(S_{i}-S_{i}^{{}^{\prime}})^{4}+(T_{i}-T_{i}^{{}^{\prime}})^{4}-2(S_{i}-T_{i})^{4}\right]$
$\displaystyle+2\mathrm{E}\sum_{1\leq i<j\leq
d}\Big{[}(S_{i}-S_{i}^{{}^{\prime}})^{2}(S_{j}-S_{j}^{{}^{\prime}})^{2}+(T_{i}-T_{i}^{{}^{\prime}})^{2}(T_{j}-T_{j}^{{}^{\prime}})^{2}-2(S_{i}-T_{i})^{2}(S_{j}-T_{j})^{2}\Big{]}$
$\displaystyle=2\mathrm{E}\sum_{1\leq i<j\leq
d}\Big{[}(S_{i}-S_{i}^{{}^{\prime}})^{2}(S_{j}-S_{j}^{{}^{\prime}})^{2}+(T_{i}-T_{i}^{{}^{\prime}})^{2}(T_{j}-T_{j}^{{}^{\prime}})^{2}-2(S_{i}-T_{i})^{2}(S_{j}-T_{j})^{2}\Big{]}.$

Hence
$\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{X}},\Pi)=\frac{1}{4\sigma^{4}}\mathrm{E}\sum_{1\leq
i<j\leq
d}\Big{[}(S_{i}-S_{i}^{{}^{\prime}})^{2}(S_{j}-S_{j}^{{}^{\prime}})^{2}+(T_{i}-T_{i}^{{}^{\prime}})^{2}(T_{j}-T_{j}^{{}^{\prime}})^{2}-2(S_{i}-T_{i})^{2}(S_{j}-T_{j})^{2}\Big{]}+\mathcal{O}\left(\frac{1}{\sigma^{6}}\right).$
Now, by some straight-forward but tedious calculations, it can be shown that
$\mathrm{E}\Big{[}(S_{i}-S_{i}^{{}^{\prime}})^{2}(S_{j}-S_{j}^{{}^{\prime}})^{2}+(T_{i}-T_{i}^{{}^{\prime}})^{2}(T_{j}-T_{j}^{{}^{\prime}})^{2}-2(S_{i}-T_{i})^{2}(S_{j}-T_{j})^{2}\Big{]}=4\mathrm{Cov}^{2}(S_{i},S_{j})$.
This implies that as $\sigma\to\infty$,
$\sigma^{4}\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{X}},\Pi)\to\sum_{1\leq
i<j\leq d}\mathrm{Cov}^{2}(S_{i},S_{j}).$ ∎

######  Proof (Proof of Lemma 5)

(a) Clearly, applying a permutation to the coordinates of the observation
vectors $\mathbf{X}^{(i)},i=1,2,\cdots,n,$ changes the coordinates of the
$\mathbf{Y}^{(i)}$’s by the same permutation. Since $s_{1}$ and $s_{2}$ from
equation (7) are both invariant under permutation of coordinates of the
$\mathbf{Y}^{(i)}$’s, the proof is complete.

For proving invariance under monotonic transformation, it is enough to
consider the case when only one of the coordinates in the observation vectors
is changed by a strictly monotonic non-identity transformation. Assume,
therefore, that only the $s^{\text{th}}$ coordinate of the
$\mathbf{X}^{(i)}$’s is changed by a strictly monotonic transformation, while
the other coordinates are kept the same. This will affect only the
$s^{\text{th}}$ coordinate of the $\mathbf{Y}^{(i)}$’s. Denoting the changed
$\mathbf{Y}^{(i)}$’s as $\mathbf{Y}^{*(i)}$’s, it is clear that $Y^{*(i)}_{s}$
will equal $Y^{(i)}_{s}$ or $1+\frac{1}{n}-Y^{(i)}_{s}$ for all
$i=1,2,\cdots,n$, according as the transformation is strictly increasing or
strictly decreasing. In either case,
$k_{\sigma}(\mathbf{Y}^{*(i)},\mathbf{Y}^{*(j)})=k_{\sigma}(\mathbf{Y}^{(i)},\mathbf{Y}^{(j)})$,
so that $s_{1}$ in equation (7) remains unchanged. One can easily see that
$s_{2}$ also remains unchanged as well.

(b) Without loss of generality, we may assume that the first coordinates of
the $\mathbf{X}^{(i)}$’s are in ascending order. Now, suppose that every other
coordinate of the $\mathbf{X}^{(i)}$’s is in a strictly monotonic relation
with the first coordinate; then, for $j=2,\cdots,d$, the $j^{\text{th}}$
coordinates of the $\mathbf{X}^{(i)}$’s will be in either ascending or
descending order. By monotonic transformation invariance property, we may
assume, without loss of generality, that all the coordinates of the
$\mathbf{X}^{(i)}$’s are in ascending order. But then, the
$\mathbf{Y}^{(i)}$’s are clearly given by $Y_{j}^{(i)}=\frac{i}{n}$, for all
$j$ and one can then see that $s_{1}=v_{1}$ and $s_{2}=v_{2}$, whence it
follows that $\widehat{I}_{\sigma,n}(\mathbf{X})=1$. ∎

######  Proof (Proof of Theorem 4)

For two independent random vectors $(T_{1},T_{2})$ and
$(T_{1}^{{}^{\prime}},T_{2}^{{}^{\prime}})$ with both having distribution
$\mathsf{C}_{\mathbf{x},n}$, one has

$V_{1}=k_{\sigma}(T_{1},T_{1}^{{}^{\prime}})-\mathrm{E}\left[k_{\sigma}(T_{1},T_{1}^{{}^{\prime}})\Big{|}T_{1}\right]-\mathrm{E}\left[k_{\sigma}(T_{1},T_{1}^{{}^{\prime}})\Big{|}T_{1}^{{}^{\prime}}\right]+\mathrm{E}\left[k_{\sigma}(T_{1},T_{1}^{{}^{\prime}})\right]\mbox{
and}$
$V_{2}=k_{\sigma}(T_{2},T_{2}^{{}^{\prime}})-\mathrm{E}\left[k_{\sigma}(T_{2},T_{2}^{{}^{\prime}})\Big{|}T_{2}\right]-\mathrm{E}\left[k_{\sigma}(T_{2},T_{2}^{{}^{\prime}})\Big{|}T_{2}^{{}^{\prime}}\right]+\mathrm{E}\left[k_{\sigma}(T_{2},T_{2}^{{}^{\prime}})\right]\mbox{
so that}\vspace{-1.5ex}$
$\displaystyle\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\Pi_{n})$
$\displaystyle=\mathrm{E}\left[k_{\sigma}(T_{1},T_{1}^{{}^{\prime}})k_{\sigma}(T_{2},T_{2}^{{}^{\prime}})\right]-2\mathrm{E}\left[\mathrm{E}\left[k_{\sigma}(T_{1},T_{1}^{{}^{\prime}})\Big{|}T_{1}\right]\mathrm{E}\left[k_{\sigma}(T_{2},T_{2}^{{}^{\prime}})\Big{|}T_{2}\right]\right]$
$\displaystyle+\mathrm{E}\left[k_{\sigma}(T_{1},T_{1}^{{}^{\prime}})\right]\mathrm{E}\left[k_{\sigma}(T_{2},T_{2}^{{}^{\prime}})\right]$
$\displaystyle=\mathrm{E}\left[\left\\{k_{\sigma}(T_{1},T_{1}^{{}^{\prime}})-\mathrm{E}\left[k_{\sigma}(T_{1},T_{1}^{{}^{\prime}})\Big{|}T_{1}\right]-\mathrm{E}\left[k_{\sigma}(T_{1},T_{1}^{{}^{\prime}})\Big{|}T_{1}^{{}^{\prime}}\right]+\mathrm{E}\left[k_{\sigma}(T_{1},T_{1}^{{}^{\prime}})\right]\right\\}\right.$
$\displaystyle\quad\quad\quad\left.\left\\{k_{\sigma}(T_{2},T_{2}^{{}^{\prime}})-\mathrm{E}\left[k_{\sigma}(T_{2},T_{2}^{{}^{\prime}})\Big{|}T_{2}\right]-\mathrm{E}\left[k_{\sigma}(T_{2},T_{2}^{{}^{\prime}})\Big{|}T_{2}^{{}^{\prime}}\right]+\mathrm{E}\left[k_{\sigma}(T_{2},T_{2}^{{}^{\prime}})\right]\right\\}\right]$
$\displaystyle=\frac{1}{n^{2}}\sum_{1\leq i,j\leq n}V_{1}(i,j)V_{2}(i,j)\,.$

Here the second last equality follows from Lemma L1. Similarly, one can show
that

$\gamma^{2}_{k_{\sigma}}(\mathsf{M}_{n},\Pi_{n})=\frac{1}{n^{2}}\sum_{1\leq
i,j\leq n}V_{1}^{2}(i,j)=\frac{1}{n^{2}}\sum_{1\leq i,j\leq
n}V_{2}^{2}(i,j)\,.$

Cauchy-Schwartz (CS) inequality immediately gives
$\widehat{I}_{\sigma,n}(\mathbf{X})^{2}\leq 1$. Further, by the necessary and
sufficient condition for equality in the CS inequality and using the fact that
$\sum_{1\leq i,j\leq n}V^{2}_{1}(i,j)=\sum_{1\leq i,j\leq n}V^{2}_{2}(i,j)$,
one gets that $\widehat{I}_{\sigma,n}(\mathbf{X})=1$ if and only if
$V_{1}(i,j)=V_{2}(i,j)~{}\forall~{}i,j$.

Now, if one coordinate of the observation vectors is a monotone function of
the other coordinate, then either $Y_{2}^{(i)}=Y_{1}^{(i)}~{}\forall~{}i$ or
$Y_{2}^{(i)}=\frac{n+1}{n}-Y_{1}^{(i)}~{}\forall~{}i$. In either case,
$|Y_{1}^{(i)}-Y_{1}^{(j)}|=|Y_{2}^{(i)}-Y_{2}^{(j)}|~{}\forall~{}i,j$, which
will clearly imply that $V_{1}(i,j)=V_{2}(i,j)~{}\forall~{}i,j$.

To prove the converse, first observe that for any $i$,

$\displaystyle\sum_{l=1}^{n}k_{\sigma}(Y_{1}^{(i)},Y_{1}^{(l)})=\sum_{l=1}^{n}k_{\sigma}(Y_{1}^{(l)},Y_{1}^{(i)})=\sum_{l=1}^{n}k_{\sigma}(Y_{1}^{(i)},l/n)\text{~{}and~{}}$
$\displaystyle\sum_{l=1}^{n}k_{\sigma}(Y_{2}^{(i)},Y_{2}^{(l)})=\sum_{l=1}^{n}k_{\sigma}(Y_{2}^{(l)},Y_{2}^{(i)})=\sum_{l=1}^{n}k_{\sigma}(Y_{2}^{(i)},l/n).$

Now suppose that $V_{1}(i,j)=V_{2}(i,j)~{}\forall~{}i,j$. Then, taking $i=j$,
one deduces that

$\sum_{l=1}^{n}k_{\sigma}(Y^{(i)}_{1},l/n)=\sum_{l=1}^{n}k_{\sigma}(Y^{(i)}_{2},l/n)\quad\forall~{}i\in\\{1,2,\cdots,n\\}.$
(E4)

Using this now in $V_{1}(i,j)=V_{2}(i,j)$, one gets

$k_{\sigma}(Y_{1}^{(i)},Y_{1}^{(j)})=k_{\sigma}(Y_{2}^{(i)},Y_{2}^{(j)}),\text{
i.e.,
}|Y_{1}^{(i)}-Y_{1}^{(j)}|=|Y_{2}^{(i)}-Y_{2}^{(j)}|~{}~{}\forall~{}i,j.$ (E5)

We now claim that for $i,i^{{}^{\prime}}\in\\{1,2,\cdots,n\\}$,
$\sum_{l=1}^{n}k_{\sigma}(i/n,l/n)=\sum_{l=1}^{n}k_{\sigma}(i^{{}^{\prime}}/n,l/n)$
if and only if either $i^{{}^{\prime}}=i$ or $i^{{}^{\prime}}=n+1-i$. The ‘if’
part is easy to see; if $i^{{}^{\prime}}=n+1-i$, the equality is obtained by
observing that
$k_{\sigma}(i^{{}^{\prime}}/n,j/n)=k_{\sigma}(i/n,(n+1-j)/n)~{}\forall~{}j$
and then making a change of variable ($j\mapsto n+1-j$) in the summation. The
’only if’ part can now be completed by observing that whenever $i<n+1-i$,
$\sum_{l=1}^{n}k_{\sigma}((i+1)/n,l/n)-\sum_{l=1}^{n}k_{\sigma}(i/n,l/n)=e^{-\frac{i^{2}}{2n^{2}\sigma^{2}}}-e^{-\frac{(n-i)^{2}}{2n^{2}\sigma^{2}}}>0$,
implying that $\sum_{l=1}^{n}k_{\sigma}(i/n,l/n)$ is strictly increasing in
$i$ whenever $i<n+1-i$.

Using this, (E4) implies that, for each $i$, we have either
$Y^{(i)}_{2}=Y^{(i)}_{1}$ or $Y^{(i)}_{2}=1+\frac{1}{n}-Y^{(i)}_{1}$. Next,
let $i$ be such that $Y_{1}^{(i)}=1/n$. We know that either
$Y_{2}^{(i)}=Y_{1}^{(i)}$ or $Y_{2}^{(i)}=1+1/n-Y_{1}^{(i)}$. Suppose first
that $Y_{2}^{(i)}=Y_{1}^{(i)}$. Now, take any $j\not=i$. We know $Y_{2}^{(j)}$
equals either $Y_{1}^{(j)}$ or $1+1/n-Y_{1}^{(j)}$. But then (E5) rules out
the possibility that $Y_{2}^{(j)}=1+1/n-Y_{1}^{(j)}$. Thus we have
$Y_{2}^{(j)}=Y_{1}^{(j)}$ for all $j$. Similarly, if
$Y_{2}^{(i)}=1+1/n-Y_{1}^{(i)}$, one can show that
$Y_{2}^{(j)}=1+1/n-Y_{1}^{(j)}$ for all $j$. Thus we conclude that either
$Y_{2}^{(j)}=Y_{1}^{(j)}~{}\forall~{}j$ or
$Y_{2}^{(j)}=1+1/n-Y_{1}^{(j)}~{}\forall~{}j$. But this means that one
coordinate of the observation vectors is either an increasing or a decreasing
function of the other coordinate. ∎

The following well-known result, which can be found in Tsukahara (2005), is
crucial in our derivation of the limiting distributions of
${\widehat{I}}_{\sigma,n}(\mathbf{X})$ in both under the null and the
alternative hypotheses.

######  Theorem T1 (Weak convergence of copula process)

Let $\mathbf{X^{(1)}},\cdots,\mathbf{X^{(n)}}$ be independent observations on
the random vector $\mathbf{X}$ with copula distribution
$\mathsf{C}_{\mathbf{X}}$ and let $\mathsf{C}_{\mathbf{x},n}$ be the empirical
copula based on $\mathbf{X^{(1)}},\cdots,\mathbf{X^{(n)}}$. If, for all
$i=1,2,\cdots,d$, the $i^{\text{th}}$ partial derivatives
$D_{i}\mathsf{C}_{\mathbf{X}}(u)$ of $\mathsf{C}_{\mathbf{X}}$ exist and are
continuous, then the process
$\sqrt{n}(\mathsf{C}_{\mathbf{x},n}-\mathsf{C}_{\mathbf{X}})$ converges weakly
in $l^{\infty}([0,1]^{d})$ to the process
$\mathbb{G}_{\mathsf{C}_{\mathbf{X}}}$ given by

$\mathbb{G}_{\mathsf{C}_{\mathbf{X}}}(u)=\mathbb{B}_{\mathsf{C}_{\mathbf{X}}}(u)-\sum_{i=1}^{d}D_{i}\mathsf{C}_{\mathbf{X}}(u)\mathbb{B}_{\mathsf{C}_{\mathbf{X}}}(u^{(i)}),\vspace{-2ex}$

where $\mathbb{B}_{\mathsf{C}_{\mathbf{X}}}$ is a $d$-dimensional Brownian
bridge on $[0,1]^{d}$ with covariance function
$\mathrm{E}[\mathbb{B}_{\mathsf{C}_{\mathbf{X}}}(u)\mathbb{B}_{\mathsf{C}_{\mathbf{X}}}(v)]=\mathsf{C}_{\mathbf{X}}(u)\wedge\mathsf{C}_{\mathbf{X}}(v)-\mathsf{C}_{\mathbf{X}}(u)\mathsf{C}_{\mathbf{X}}(v)$,
and for each $i$, $u^{(i)}$ denotes the vector obtained from $u$ by replacing
its all coordinates, except the $i^{\text{th}}$ th one, by $1$.

Let the distribution function of $\mathbf{X}$ be denoted by $\mathbf{F}$ and
its marginals by $F_{1},F_{2},\cdots,F_{d}$. With
$\mathbf{X}^{(i)}\\!=\\!(X_{1}^{(i)},X_{2}^{(i)},\cdots
X_{d}^{(i)}),\;1\\!\leq\\!i\\!\leq\\!n$, denoting i.i.d. observations from
$\mathbf{X}$, define vectors
$\mathbf{Z}^{(i)}\\!=\\!(Z_{1}^{(i)},Z_{2}^{(i)},\allowbreak\cdots,Z_{d}^{(i)}),\;1\\!\leq\\!i\\!\leq\\!n$,
where $Z_{j}^{(i)}=F_{j}(X_{j}^{(i)})$. We will denote the empirical
distribution based on
$\mathbf{Z}^{(1)},\mathbf{Z}^{(2)},\cdots,\mathbf{Z}^{(n)}$ by
$\mathsf{C}_{\mathbf{z},n}$, and the empirical distribution function based on
the $\mathbf{X}^{(i)}$’s by $\hat{\mathbf{F}}$.

###### Lemma L2

Assume that $\\{P_{n}\\}_{n\geq 1}$ is a sequence of distributions over
$[0,1]^{d}$. Then

  1. 1.

$\left|\gamma^{2}_{k_{\sigma}}(\Pi_{n},P_{n})-\gamma^{2}_{k_{\sigma}}(\Pi,P_{n})\right|=\mathcal{O}(n^{-2})$

  2. 2.

$\left|\gamma^{2}_{k_{\sigma}}(\mathsf{M}_{n},P_{n})-\gamma^{2}_{k_{\sigma}}(\mathsf{M},P_{n})\right|=\mathcal{O}(n^{-2}).$

###### Proof

We prove the first part only. The proof of the second part is similar.

$\displaystyle\left|\gamma^{2}_{k_{\sigma}}(\Pi_{n},P_{n})\\!-\\!\gamma^{2}_{k_{\sigma}}(\Pi,P_{n})\right|$
$\displaystyle\leq\left|\mathrm{E}_{{}_{(S,S^{{}^{\prime}})\sim\Pi_{n}\otimes\Pi_{n}}}\\![k_{\sigma}(S,S^{{}^{\prime}})]-\mathrm{E}_{{}_{(S,S^{{}^{\prime}})\sim\Pi\otimes\Pi}}\\![k_{\sigma}(S,S^{{}^{\prime}})]\right|$
$\displaystyle\quad+2\mathrm{E}_{{}_{T\sim
P_{n}}}\left|\mathrm{E}_{{}_{S\sim\Pi_{n}}}\\![k_{\sigma}(S,T)]-\mathrm{E}_{{}_{S\sim\Pi}}\\![k_{\sigma}(S,T)]\right|.$

The first term on the right hand side of the above inequality is easily seen
to be bounded above by

$\frac{1}{n^{2d}}\\!\\!\sum_{\begin{subarray}{c}\mu=(i_{1}/n,\cdots,i_{d}/n)\\\
1\leq i_{1},\cdots,i_{d}\leq
n\end{subarray}}\\!\,\sum_{\begin{subarray}{c}\nu=(j_{1}/n,\cdots,j_{d}/n)\\\
1\leq j_{1},\cdots,j_{d}\leq
n\end{subarray}}\;\int\limits_{[\mu-1/n,\mu]}\int\limits_{[\nu-1/n,\nu]}\\!\\!\\!\\!\\!\left|k_{\sigma}(\mu,\nu)\\!-\\!k_{\sigma}(\zeta,\eta)\right|d\zeta\,d\eta,\vspace{-1ex}$

where for any $u=(u_{1},u_{2},\cdots,u_{d})\in[0,1]^{d}$ and $\delta>0$,
$[u-\delta,u]$ denotes the rectangle
$[u_{1}-\delta,u_{1}]\times[u_{2}-\delta,u_{2}]\times\cdots\times[u_{d}-\delta,u_{d}]$.
The last expression is clearly bounded above by

$\max_{\begin{subarray}{c}\mu=(i_{1}/n,\cdots,i_{d}/n)\\\ 1\leq
i_{1},\cdots,i_{d}\leq
n\end{subarray}}\;\max_{\begin{subarray}{c}\nu=(j_{1}/n,\cdots,j_{d}/n)\\\
1\leq j_{1},\cdots,j_{d}\leq
n\end{subarray}}\;\sup_{\zeta\in[\mu-1/n,\mu]}\;\sup_{\eta\in[\nu-1/n,\nu]}\left|k_{\sigma}(\mu,\nu)-k_{\sigma}(\zeta,\eta)\right|.\vspace{-1ex}$

Using Lemma 6 of Póczos et al. (2012), one can further deduce that the last
expression is bounded above by $dn^{-2}$.
Similar technique can be used for the second term to get the upper bound

$2\,\mathrm{E}_{T\sim
P_{n}}\;\max_{\begin{subarray}{c}\mu=(i_{1}/n,i_{2}/n,\cdots,i_{d}/n)\\\ 1\leq
i_{1},i_{2},\cdots,i_{d}\leq
n\end{subarray}}\;\sup_{\eta\in[\mu-1/n,\mu]}\left|k_{\sigma}(\mu,T)-k_{\sigma}(\eta,T)\right|\leq
2dn^{-2}.\vspace{-1ex}$

Combining these two bounds, we get
$\left|\gamma^{2}_{k_{\sigma}}(\Pi_{n},P_{n})-\gamma^{2}_{k_{\sigma}}(\Pi,P_{n})\right|=\mathcal{O}(n^{-2})$.
∎

###### Lemma L3

$\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\mathsf{C}_{\mathbf{z},n})\rightarrow
0$ almost surely as $n\rightarrow\infty$.

######  Proof (Sketch of the proof)

Since the essential idea of the proof is contained in Póczos et al. (2012)
[Appendix E], we only describe the two main steps.

First, we use the definition of $\mathsf{C}_{\mathbf{z},n}$ and Lemma 6 of
Póczos et al. (2012) to get the inequality
$\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\mathsf{C}_{\mathbf{z},n})\leq
2\sqrt{d}L\,\max\limits_{1\leq j\leq
d}\;\sup\limits_{x\in\mathbb{R}}|\hat{F}_{j}(x)-F_{j}(x)|$, where
$\hat{F}_{1},\ldots,\hat{F}_{d}$ are the marginals of the empirical
distribution $\hat{\mathbf{F}}$ based on
$\mathbf{X}^{(1)},\ldots,\mathbf{X}^{(n)}$.

Then using the above inequality and the Kiefer-Dvoretzky-Wolfowitz Theorem
(see Massart (1990), Page 1269), for any $\epsilon>0$, we get
$\mathrm{Pr}\left[\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\mathsf{C}_{\mathbf{z},n})>\epsilon\right]\leq
2d\exp\left(-\frac{n\epsilon^{2}}{2dL^{2}}\right)$. The result now follows
from the Borel-Cantelli Lemma. ∎

The next lemma and its proof are based on the ideas in Gretton et al. (2012)
[Appendix A2].

###### Lemma L4

$\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{z},n},\mathsf{C}_{\mathbf{X}})\rightarrow
0$ almost surely as $n\rightarrow\infty$.

###### Proof

It is enough to prove
$\mathrm{E}\left[\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{z},n},\mathsf{C}_{\mathbf{X}})\right]\leq\frac{2}{\sqrt{n}}\;\mbox{
and }\;$
$\mathrm{Pr}\left[\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{z},n},\mathsf{C}_{\mathbf{X}})-\mathrm{E}[\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{z},n},\mathsf{C}_{\mathbf{X}})]>\epsilon\right]\leq\exp\left(\\!-\frac{n\epsilon^{2}}{2}\right).$

Denoting $\mathcal{F}$ to be the unit ball in the RKHS associated to the
kernel $k_{\sigma}$ on $\mathbb{R}^{d}$, one gets
$\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{z},n},\mathsf{C}_{\mathbf{X}})=\sup_{f\in\mathcal{F}}\left|\frac{1}{n}\sum_{i=1}^{n}f(\mathbf{Z}^{(i)})-\mathrm{E}_{\mathbf{Z}\sim\mathsf{C}_{\mathbf{X}}}f(\mathbf{Z})\right|$
(see Sriperumbudur et al. (2010)).

Letting $\mathbf{Z}^{*(i)},\;1\leq i\leq n$ to be i.i.d. with the same
distribution as and independent of $\mathbf{Z}^{(i)},\;1\leq i\leq n$ and
$\delta_{i},\;1\leq i\leq n$ to be i.i.d. random variables taking values $\pm
1$ with equal probabilities, independent of the
$\mathbf{Z}^{(i)},\mathbf{Z}^{*(i)},\;1\leq i\leq n$, it is easy to see that

$\displaystyle\mathrm{E}[\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{z},n},\mathsf{C}_{\mathbf{X}})]$
$\displaystyle=\mathrm{E}\left[\sup_{f\in\mathcal{F}}\left|\frac{1}{n}\sum_{i=1}^{n}f(\mathbf{Z}^{(i)})-\mathrm{E}_{\mathbf{Z}\sim\mathsf{C}_{\mathbf{X}}}f(\mathbf{Z})\right|\right]$
$\displaystyle\leq\mathrm{E}\left[\sup_{f\in\mathcal{F}}\left|\frac{1}{n}\sum_{i=1}^{n}f(\mathbf{Z}^{(i)})-\frac{1}{n}\sum_{i=1}^{n}f(\mathbf{Z}^{*(i)})\right|\right]$
$\displaystyle=\mathrm{E}\left[\sup_{f\in\mathcal{F}}\left|\frac{1}{n}\sum_{i=1}^{n}\delta_{i}\left(f(\mathbf{Z}^{(i)})-f(\mathbf{Z}^{*(i)})\right)\right|\right]\underset{(a)}{\leq}\frac{2}{\sqrt{n}}.\vspace{-10pt}$

For the last inequality (a), we used a well-known result referred to as “Bound
on Rademacher Complexity” (see Bartlett and Mendelson (2003), Page 478) .

We next calculate the upper bound of change in magnitude due to change in a
particular coordinate. Consider
$\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{z},n},\mathsf{C}_{\mathbf{X}})$ as a
function of $\mathbf{Z}^{(i)}$. It is easy to verify that changing any
coordinate $\mathbf{Z}^{(i)}$, the change in
$\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{z},n},\mathsf{C}_{\mathbf{X}})$ will
be at most $2n^{-1}$. We use now the well-known McDiarmid’s inequality (see
McDiarmid (1989), Page 149) to get

$\mathrm{Pr}\left[\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{z},n},\mathsf{C}_{\mathbf{X}})-\mathrm{E}[\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{z},n},\mathsf{C}_{\mathbf{X}})]>\epsilon\right]\leq\exp\left(-\frac{2\epsilon^{2}}{n.(2/n)^{2}}\right)=\exp\left(-\frac{n\epsilon^{2}}{2}\right).\qed$

###### Lemma L5

Suppose that the assumptions of Theorem T1 hold. Then, we have the following
results.

  * If $\mathsf{C}_{\mathbf{X}}\neq\Pi$, $\sqrt{n}(\gamma_{k_{\sigma}}^{2}(\mathsf{C}_{\mathbf{x},n},\Pi)-\gamma_{k_{\sigma}}^{2}(\mathsf{C}_{\mathbf{X}},\Pi))\overset{\mathcal{L}}{\rightarrow}\mathcal{N}(0,\delta_{0}^{2}),$ where
$\delta^{2}_{0}=$
$4\displaystyle\int_{\mathbb{R}^{d}}\displaystyle\int_{\mathbb{R}^{d}}{g(u)}g(v)\,\mathrm{E}[\,d\mathbb{G}_{\mathsf{C}_{\mathbf{X}}}(u)\,d\mathbb{G}_{\mathsf{C}_{\mathbf{X}}}(v)]$
and
$g(u)=\displaystyle\int_{[0,1]^{d}}{k_{\sigma}(u,v)\,d(\mathsf{C}_{\mathbf{X}}-\Pi)(v)}$.

  * If $\mathsf{C}_{\mathbf{X}}=\Pi$, $n\gamma_{k_{\sigma}}^{2}(\mathsf{C}_{\mathbf{x},n},\Pi)\overset{\mathcal{L}}{\rightarrow}\displaystyle\int_{\mathbb{R}^{d}}\displaystyle\int_{\mathbb{R}^{d}}k_{\sigma}(u,v)\,d\mathbb{G}_{\Pi}(u)\,d\mathbb{G}_{\Pi}(v).$

###### Proof

When $\mathsf{C}_{\mathbf{X}}\neq\Pi$:  Denoting $\mathcal{D}([0,1]^{d})$ to
be the space of right continuous real valued uniformly bounded functions on
$[0,1]^{d}$ with left limits, equipped with max-sup norm, one can easily
verify that the function
$\psi(\mathsf{D})=\gamma^{2}_{k_{\sigma}}(\mathsf{D},\Pi)$ on
$\mathcal{D}([0,1]^{d})$ is Hadamard-differentiable and the derivative at
$\mathsf{C}_{\mathbf{X}}$ is given by

$\psi^{{}^{\prime}}_{\mathsf{C}_{\mathbf{X}}}(\mathsf{D})=2\int_{[0,1]^{d}}\int_{[0,1]^{d}}k_{\sigma}(u,w)\,d(\mathsf{C}_{\mathbf{X}}-\Pi)(w)\,d\mathsf{D}(u).$

To prove this, consider a real sequence $\\{t_{n}\\}$ converging to 0 and a
$\mathcal{D}([0,1]^{d})$-valued sequence $\\{D_{n}\\}$ converging to
$D\in\mathcal{D}([0,1]^{d})$ such that
$\mathsf{C}_{\mathbf{X}}+t_{n}D_{n}\in\mathcal{D}([0,1]^{d})$. For any
$D\in\mathcal{D}([0,1]^{d})$, define
$\mu_{D}(w)=\int_{[0,1]^{d}}k_{\frac{\sigma}{\sqrt{2}}}(u,w)\,dD(u)~{}\forall
w\in\mathbb{R}^{d}$. Then

$\displaystyle\frac{\varphi(\mathsf{C}_{\mathbf{X}}+t_{n}D_{n})-\varphi(\mathsf{C}_{\mathbf{X}})}{t_{n}}$
$\displaystyle=\left(\frac{1}{\sigma}\sqrt{\frac{2}{\pi}}\right)^{d}\frac{1}{t_{n}}\int_{\mathbb{R}^{d}}\left(\mu_{\mathsf{C}_{\mathbf{X}}}(w)+t_{n}\mu_{D_{n}}(w)-\mu_{\Pi}(w)\right)^{2}\,dw$
$\displaystyle\qquad-\left(\frac{1}{\sigma}\sqrt{\frac{2}{\pi}}\right)^{d}\frac{1}{t_{n}}\int_{\mathbb{R}^{d}}\left(\mu_{\mathsf{C}_{\mathbf{X}}}(w)-\mu_{\Pi}(w)\right)^{2}\,dw$
$\displaystyle=\left(\frac{1}{\sigma}\sqrt{\frac{2}{\pi}}\right)^{d}\frac{1}{t_{n}}\int_{\mathbb{R}^{d}}t_{n}\mu_{D_{n}}(w)\left(2\mu_{\mathsf{C}_{\mathbf{X}}}(w)+t_{n}\mu_{D_{n}}(w)-2\mu_{\Pi}(w)\right)\,dw$
$\displaystyle=\left(\frac{1}{\sigma}\sqrt{\frac{2}{\pi}}\right)^{d}\left[2\int_{\mathbb{R}^{d}}\mu_{D_{n}}(w)\left(\mu_{\mathsf{C}_{\mathbf{X}}}(w)-\mu_{\Pi}(w)\right)\,dw+t_{n}\int_{\mathbb{R}^{d}}\mu^{2}_{D_{n}}(w)\,dw\right]$
(8)

Now, using the fact
$\int_{\mathbb{R}^{d}}\left(\frac{1}{\sigma}\sqrt{\frac{2}{\pi}}\right)^{d}k_{\frac{\sigma}{\sqrt{2}}}(u,w)k_{\frac{\sigma}{\sqrt{2}}}(v,w)\,dw=k_{\sigma}(u,v)$,
it is quite straight-forward to check that
$\left(\frac{1}{\sigma}\sqrt{\frac{2}{\pi}}\right)^{d}\int_{\mathbb{R}^{d}}\mu_{D_{n}}(w)\left(\mu_{\mathsf{C}_{\mathbf{X}}}(w)-\mu_{\Pi}(w)\right)\,dw=$
$\int_{[0,1]^{d}}\int_{[0,1]^{d}}k_{\sigma}(u,v)\,dD_{n}(u)\,d(\mathsf{C}_{\mathbf{X}}-\Pi)(v)$.
From this identity and Equation (8), we get

$\displaystyle\psi^{{}^{\prime}}_{\mathsf{C}_{\mathbf{X}}}(\mathsf{D})$
$\displaystyle=\lim_{n\to\infty}\frac{\varphi(\mathsf{C}_{\mathbf{X}}+t_{n}D_{n})-\varphi(\mathsf{C}_{\mathbf{X}})}{t_{n}}$
$\displaystyle=2\int_{[0,1]^{d}}\int_{[0,1]^{d}}k_{\sigma}(u,v)\,dD(u)\,d(\mathsf{C}_{\mathbf{X}}-\Pi)(v).$

This Lemma then follows easily from Theorem T1 and the functional delta
method. The only thing that one needs to verify is that
$\psi^{{}^{\prime}}_{\mathsf{C}_{\mathbf{X}}}(\mathbb{G}_{\mathsf{C}_{\mathbf{X}}})$
is a normally distributed with 0 mean and variance
$\delta^{2}_{0}=4\int_{[0,1]^{d}}\int_{[0,1]^{d}}g(u)g(v)\,\mathrm{E}[\,d\mathbb{G}_{\mathsf{C}_{\mathbf{X}}(}u)\,d\mathbb{G}_{\mathsf{C}_{\mathbf{X}}}(v)]$
where
$g(u)=\int_{[0,1]^{d}}k_{\sigma}(u,w)\,d(\mathsf{C}_{\mathbf{X}}-\Pi)(w)$. But
this is straightforward from the formula for the derivative
$\psi^{{}^{\prime}}_{\mathsf{C}_{\mathbf{X}}}$.

When $\mathsf{C}_{\mathbf{X}}=\Pi$:  Clearly the map
$\mathsf{D}\rightarrow\left(\frac{1}{\sigma}\sqrt{\frac{2}{\pi}}\right)^{d}\\!\int_{\mathbb{R}^{d}}\left(\int_{[0,1]^{d}}k_{\frac{\sigma}{\sqrt{2}}}(u,v)\,d\mathsf{D}(u)\right)^{2}\\!dv$
from $\mathcal{D}([0,1]^{d})$ to $\mathbb{R}$ is continuous. So, the fact that
$\sqrt{n}(\mathsf{C}_{\mathbf{x},n}-\Pi)\rightarrow\mathbb{G}_{\Pi}$ and the
continuous mapping theorem gives

$\displaystyle n\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\Pi)$
$\displaystyle\overset{\mathcal{L}}{\longrightarrow}\left(\frac{1}{\sigma}\sqrt{\frac{2}{\pi}}\right)^{d}.\int_{\mathbb{R}^{d}}\left(\int_{[0,1]^{d}}k_{\frac{\sigma}{\sqrt{2}}}(u,v)\,d\mathbb{G}_{\Pi}(u)\right)^{2}\,dv$
$\displaystyle=\int_{[0,1]^{d}}\int_{[0,1]^{d}}k_{\sigma}(u,v)\,d\mathbb{G}_{\Pi}(u)\,d\mathbb{G}_{\Pi}(v).\quad\quad\quad\quad\quad\quad\quad\quad~{}~{}~{}~{}~{}~{}~{}~{}~{}\qed$

######  Proof (Proof of Theorem 5)

When $\mathsf{C}_{\mathbf{X}}\neq\Pi$:  Write
$\sqrt{n}(\widehat{I}_{\sigma,n}^{2}(\mathbf{X})\\!-\\!I^{2}_{\sigma})\\!=\\!\sqrt{n}\left(\frac{\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\Pi_{n})}{\gamma^{2}_{k_{\sigma}}(\mathsf{M}_{n},\Pi_{n})}\\!-\\!\frac{\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{X}},\Pi)}{\gamma^{2}_{k_{\sigma}}(\mathsf{M},\Pi)}\right)$
as
$A_{1,n}+A_{2,n}+A_{3,n}$, where

$\displaystyle
A_{1,n}=\sqrt{n}\left(\frac{\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\Pi_{n})}{\gamma^{2}_{k_{\sigma}}(\mathsf{M}_{n},\Pi_{n})}-\frac{\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\Pi_{n})}{\gamma^{2}_{k_{\sigma}}(\mathsf{M},\Pi)}\right),$
$\displaystyle
A_{2,n}=\sqrt{n}\left(\frac{\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\Pi_{n})}{\gamma^{2}_{k_{\sigma}}(\mathsf{M},\Pi)}-\frac{\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\Pi)}{\gamma^{2}_{k_{\sigma}}(\mathsf{M},\Pi)}\right)\\!$
$\displaystyle\mbox{ and
}A_{3,n}=\sqrt{n}\left(\frac{\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\Pi)}{\gamma^{2}_{k_{\sigma}}(\mathsf{M},\Pi)}-\frac{\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{X}},\Pi)}{\gamma^{2}_{k_{\sigma}}(\mathsf{M},\Pi)}\right).\vspace{-10pt}$

Clearly $A_{2,n}\rightarrow 0$ almost surely by Lemma L2. The same is true of
$A_{1,n}$ as well, once again by Lemma L2, because it is bounded above by

$\underbrace{\frac{\gamma^{2}_{k_{\sigma}}\\!(\mathsf{C}_{\mathbf{x},n},\\!\Pi_{n})}{\gamma^{2}_{k_{\sigma}}\\!(\mathsf{M}_{n},\\!\Pi_{n})\gamma^{2}_{k_{\sigma}}\\!(\mathsf{M},\\!\Pi)}}_{\text{bounded
sequence}}\left(\underbrace{\sqrt{n}\left|\gamma^{2}_{k_{\sigma}}\\!(\mathsf{M}_{n},\\!\Pi_{n})\\!-\\!\gamma^{2}_{k_{\sigma}}\\!(\mathsf{M},\\!\Pi_{n})\right|}_{\text{goes
to
0}}+\underbrace{\sqrt{n}\left|\gamma^{2}_{k_{\sigma}}\\!(\mathsf{M},\\!\Pi_{n})\\!-\\!\gamma^{2}_{k_{\sigma}}\\!(\mathsf{M},\\!\Pi)\right|}_{\text{goes
to 0}}\right).\vspace{-10pt}$

Therefore, using Lemma L5, we can conclude that
$\sqrt{n}(\widehat{I}_{\sigma,n}^{2}(\mathbf{X})-I^{2}_{\sigma}(\mathbf{X}))\overset{\mathcal{L}}{\rightarrow}\mathcal{N}(0,C_{\sigma,d}^{-2}.\delta^{2}_{0})$.
Now, applying the delta method, one gets
$\sqrt{n}(\widehat{I}_{\sigma,n}(\mathbf{X})-I_{\sigma}(\mathbf{X}))\overset{\mathcal{L}}{\rightarrow}\mathcal{N}(0,\delta^{2})$,
where
$\delta^{2}=C_{\sigma,d}^{-2}I^{-2}_{\sigma}(\mathbf{X})\int_{[0,1]^{d}}\int_{[0,1]^{d}}\\!g(u)g(v)\,\mathrm{E}[\,d\mathbb{G}_{\mathsf{C}_{\mathbf{X}}}(u)\,d\mathbb{G}_{\mathsf{C}_{\mathbf{X}}}(v)].$

When $\mathsf{C}_{\mathbf{X}}=\Pi$:  As a consequence of the Lemma L5 and
Lemma L2, under null hypothesis and assumptions of Theorem T1, we have

$\displaystyle n\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\Pi_{n})$
$\displaystyle=n\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\Pi)+\left(n\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\Pi_{n})-n\gamma^{2}_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\Pi)\right)$
$\displaystyle\overset{\mathcal{L}}{\longrightarrow}\;\int_{[0,1]^{d}}\int_{[0,1]^{d}}k_{\sigma}(u,v)\,d\mathbb{G}_{\Pi}(u)\,d\mathbb{G}_{\Pi}(v).\vspace{-10pt}$

It is enough to show that
$\int_{[0,1]^{d}}\int_{[0,1]^{d}}k_{\sigma}(u,v)\,d\mathbb{G}_{\Pi}(u)\,d\mathbb{G}_{\Pi}(v)\overset{\mathcal{L}}{=}\sum_{i=1}^{\infty}\alpha_{i}Z_{i}^{2},$
for some $\alpha_{i}>0$ and
$Z_{i}\overset{\text{i.i.d.}}{\sim}\mathcal{N}(0,1)$. The actual result will
follow putting $\lambda_{i}=\alpha_{i}C_{\sigma,d}$.

To this end, we define
$X(w):=\int_{[0,1]^{d}}k_{\frac{\sigma}{\sqrt{2}}}(u,w)\,d\mathbb{G}_{\Pi}(u),~{}\forall
w\in\mathbb{R}^{d}$. So, $\\{X(w):w\in\mathbb{R}^{d}\\}$ is then a zero-mean
continuous path Gaussian process. This implies that

$\int_{\mathbb{R}^{d}}\left(\int_{[0,1]^{d}}k_{\frac{\sigma}{\sqrt{2}}}(u,w)\,d\mathbb{G}_{\Pi}(u)\right)^{2}\\!dw\;=\int_{\mathbb{R}^{d}}(X(w))^{2}\,dw\;\overset{\mathcal{L}}{=}\;\sum_{i=1}^{\infty}\beta_{i}Z_{i}^{2}.\vspace{-1.5ex}$

Now, using the fact that
$\left(\frac{1}{\sigma}\sqrt{\frac{2}{\pi}}\right)^{d}\int_{\mathbb{R}^{d}}k_{\frac{\sigma}{\sqrt{2}}}(u,w)k_{\frac{\sigma}{\sqrt{2}}}(v,w)dw=k_{\sigma}(u,v)$,
one can easily see that the last equality yields the desired result

$\int_{[0,1]^{d}}\int_{[0,1]^{d}}k_{\sigma}(u,v)\,d\mathbb{G}_{\Pi}(u)\,d\mathbb{G}_{\Pi}(v)\;\overset{\mathcal{L}}{=}\;\sum_{i=1}^{\infty}\alpha_{i}Z_{i}^{2},\;\;\mbox{
with
}\alpha_{i}=\left(\frac{1}{\sigma}\sqrt{\frac{2}{\pi}}\right)^{d}\\!\\!\beta_{i}.~{}~{}~{}\qed\vspace{-1.5ex}$

######  Proof (Proof of Theorem 6)

Triangle inequality and $|a-b|^{2}\leq|a^{2}-b^{2}|$ for $a,b\geq 0$ gives

$\left|\gamma_{k_{\sigma}}\\!(\mathsf{M}_{n},\\!\Pi_{n})\\!-\\!\gamma_{k_{\sigma}}\\!(\mathsf{M},\\!\Pi)\right|\\!\leq\\!\left|\gamma^{2}_{k_{\sigma}}\\!(\mathsf{M}_{n},\\!\Pi_{n})\\!-\\!\gamma^{2}_{k_{\sigma}}\\!(\mathsf{M}_{n},\\!\Pi)\right|^{\\!\frac{1}{2}}+\left|\gamma^{2}_{k_{\sigma}}\\!(\mathsf{M}_{n},\\!\Pi)\\!-\\!\gamma^{2}_{k_{\sigma}}\\!(\mathsf{M},\\!\Pi)\right|^{\\!\frac{1}{2}}.$

Using Lemma L2, we get
$\lim_{n\rightarrow\infty}\left|\gamma_{k_{\sigma}}(\mathsf{M}_{n},\Pi_{n})-\gamma_{k_{\sigma}}(\mathsf{M},\Pi)\right|=0$
a.s.. Using again the same inequalities and the fact that
$\gamma_{k_{\sigma}}$ is a metric, one gets

$\displaystyle\left|\gamma_{k_{\sigma}}\\!(\mathsf{C}_{\mathbf{x},n},\\!\Pi_{n})\\!-\\!\gamma_{k_{\sigma}}\\!(\mathsf{C}_{\mathbf{X}},\\!\Pi)\right|\\!$
$\displaystyle\leq\\!\left|\gamma^{2}_{k_{\sigma}}\\!(\mathsf{C}_{\mathbf{x},n},\\!\Pi_{n})\\!-\\!\gamma^{2}_{k_{\sigma}}\\!(\mathsf{C}_{\mathbf{x},n},\\!\Pi)\right|^{\\!\frac{1}{2}}$
$\displaystyle\quad+\gamma_{k_{\sigma}}\\!(\mathsf{C}_{\mathbf{x},n},\\!\mathsf{C}_{\mathbf{z},n})+\gamma_{k_{\sigma}}\\!(\mathsf{C}_{\mathbf{z},n},\\!\mathsf{C}_{\mathbf{X}}).$

Again, using Lemmas L2, L3 and L4, we get
$\left|\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\Pi_{n})-\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{X}},\Pi)\right|\stackrel{{\scriptstyle
a.s.}}{{\rightarrow}}0$ as $n\rightarrow 0$, and as a consequence, we conclude
that as $n\rightarrow\infty$,

$\widehat{I}_{\sigma,n}(\mathbf{X})=\frac{\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{x},n},\Pi_{n})}{\gamma_{k_{\sigma}}(\mathsf{M}_{n},\Pi_{n})}\;\rightarrow\;\frac{\gamma_{k_{\sigma}}(\mathsf{C}_{\mathbf{X}},\Pi)}{\gamma_{k_{\sigma}}(\mathsf{M},\Pi)}=I_{\sigma}(\mathbf{X})\mbox{
almost surely}\quad\quad\quad\quad~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\qed.$

###### Lemma L6

Let $P_{n}$ and $Q_{n}$ be sequence of probability distribution over
$[0,1]^{d}$. Let $\sigma_{n}$ be a sequence of positive real numbers that
converges to $\sigma_{0}>0$. Then as $n\to\infty$,
$|\gamma_{k_{\sigma_{n}}}^{2}(P_{n},Q_{n})-\gamma_{k_{\sigma_{0}}}^{2}(P_{n},Q_{n})|\to
0.$

###### Proof

First we observe that

$\displaystyle|\gamma_{k_{\sigma_{n}}}^{2}(P_{n},Q_{n})-\gamma_{k_{\sigma_{0}}}^{2}(P_{n},Q_{n})|$
$\displaystyle\leq\mathrm{E}_{(\mathbf{S},\mathbf{S}^{{}^{\prime}})\sim
P_{n}\otimes
P_{n}}|k_{\sigma_{n}}(\mathbf{S},\mathbf{S}^{{}^{\prime}})-k_{\sigma_{0}}(\mathbf{S},\mathbf{S}^{{}^{\prime}})|$
$\displaystyle\quad+2\mathrm{E}_{(\mathbf{S},\mathbf{T})\sim P_{n}\otimes
Q_{n}}|k_{\sigma_{n}}(\mathbf{S},\mathbf{T})-k_{\sigma_{0}}(\mathbf{S},\mathbf{T})|$
$\displaystyle\quad+\mathrm{E}_{(\mathbf{T},\mathbf{T}^{{}^{\prime}})\sim
Q_{n}\otimes
Q_{n}}|k_{\sigma_{n}}(\mathbf{S},\mathbf{T}^{{}^{\prime}})-k_{\sigma_{0}}(\mathbf{T},\mathbf{T}^{{}^{\prime}})|.$

Applying Lemma 6 of Póczos et al. (2012), one can get an upper bound of
$|k_{\sigma_{n}}(\mathbf{S},\mathbf{T})$
$-k_{\sigma_{0}}(\mathbf{S},\mathbf{T})|$ in the following way

$\displaystyle|k_{\sigma_{n}}(\mathbf{S},\mathbf{T})-k_{\sigma_{0}}(\mathbf{S},\mathbf{T})|$
$\displaystyle\leq
L\left\|\frac{\mathbf{S}}{\sigma_{n}}-\frac{\mathbf{S}}{\sigma_{0}}\right\|+L\left\|\frac{\mathbf{T}}{\sigma_{n}}-\frac{\mathbf{T}}{\sigma_{0}}\right\|\leq
2L\sqrt{d}\left|\frac{1}{\sigma_{n}}-\frac{1}{\sigma_{0}}\right|,$

where $L$ is a constant. Thus we can conclude that
$|\gamma_{k_{\sigma_{n}}}^{2}(P_{n},Q_{n})-\gamma_{k_{\sigma_{0}}}^{2}(P_{n},Q_{n})|\leq
8L\sqrt{d}\left|\frac{1}{\sigma_{n}}-\frac{1}{\sigma_{0}}\right|.$ This
completes the proof. ∎

######  Proof (Proof of Theorem 7)

Note that

$\displaystyle|\gamma^{2}_{\sigma_{n}}(\mathsf{C}_{\mathbf{x},n},\Pi_{n})-\gamma^{2}_{\sigma_{0}}(\mathsf{C}_{\mathbf{X}},\Pi)|$
$\displaystyle\leq|\gamma^{2}_{\sigma_{n}}(\mathsf{C}_{\mathbf{x},n},\Pi_{n})-\gamma^{2}_{\sigma_{0}}(\mathsf{C}_{\mathbf{x},n},\Pi_{n})|$
$\displaystyle\quad+|\gamma^{2}_{\sigma_{0}}(\mathsf{C}_{\mathbf{x},n},\Pi_{n})-\gamma^{2}_{\sigma_{0}}(\mathsf{C}_{\mathbf{X}},\Pi)|.$

As $n\to\infty$, the first term in the right hand side goes to 0 due to Lemma
L6 and the second term almost surely converges to 0 due to Theorem 6.
Similarly, one can show that
$|\gamma^{2}_{\sigma_{n}}(\mathsf{M}_{n},\Pi_{n})-\gamma^{2}_{\sigma_{0}}(\mathsf{M},\Pi)|\to
0$ as $n\to\infty$. This implies $\widehat{I}_{\sigma_{n},n}(\mathbf{X})\to
I_{\sigma_{0}}(\mathbf{X})$ almost surely. Because of the fact that
$I_{\sigma_{0}}(\mathbf{X})=0$ if and only if the coordinates of $\mathbf{X}$
are independent, test of independence based on the statistic
$\widehat{I}_{\sigma_{n},n}(\mathbf{X})$ is consistent. ∎

######  Proof (Proof of Theorem 8)

From Theorem (7), it follows that $\widehat{I}_{\sigma^{(i)},n}(\mathbf{X})$’s
are consistent test statistics for all $i=1,2,\cdots,m$. Since $m$ is finite,
by the virtue of the definition of $T_{max}$ and $T_{sum}$, they converge to 0
almost surely if and only if the coordinates of $\mathbf{X}$ are independent.
Otherwise, they converge to positive quantities. This property makes the
resulting tests consistent. Again, under the alternative hypothesis, for any
$i$, the p-value $p_{i}$ corresponding to the test statistic
$\widehat{I}_{\sigma^{(i)},n}(\mathbf{X})$ converges to zero almost surely.
So, for sufficiently large $n$, almost surely, there would exist at least one
$i$ such that $p_{i}$ is less than $\alpha/m$, which makes the set
$\\{i:p_{(i)}<\alpha/m\\}$ non-empty. Thus the power of the test based on FDR
tends to be $1$ as sample size tends to infinity. ∎

## References

  * Anderson (2003) Anderson TW (2003) An Introduction to Multivariate Statistical Analysis. Wiley, New York
  * Bartlett and Mendelson (2003) Bartlett PL, Mendelson S (2003) Rademacher and Gaussian complexities: risk bounds and structural results. The Journal of Machine Learning Research 3:463–482
  * Benjamini and Hochberg (1995) Benjamini Y, Hochberg Y (1995) Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society Series B 57:289–300
  * Benjamini and Yekutieli (2001) Benjamini Y, Yekutieli D (2001) The control of the false discovery rate in multiple testing under dependency. Annals of Statistics 29(4):1165–1188
  * Biswas et al. (2016) Biswas M, Sarkar S, Ghosh AK (2016) On some exact distribution-free tests of independence between two random vectors of arbitrary dimensions. Journal of Statistical Planning and Inference 175:78–86
  * Blomqvist (1950) Blomqvist N (1950) On a measure of dependence between two random variables. The Annals of Mathematical Statistics 21(4):593–600
  * Cuesta-Albertos and Febrero-Bande (2010) Cuesta-Albertos J, Febrero-Bande M (2010) A simple multiway ANOVA for functional data. Test 19(3):537–557
  * Fan et al. (2017) Fan Y, de Micheaux PL, Penev S, Salopek D (2017) Multivariate nonparametric test of independence. Journal of Multivariate Analysis 153:189 – 210
  * Fukumizu et al. (2008) Fukumizu K, Gretton A, Sun X, Schölkopf B (2008) Kernel measures of conditional dependence. In: Advances in Neural Information Processing Systems, pp 489–496
  * Fukumizu et al. (2009) Fukumizu K, Gretton A, Lanckriet GR, Schölkopf B, Sriperumbudur BK (2009) Kernel choice and classifiability for rkhs embeddings of probability distributions. In: Advances in Neural Information Processing Systems, pp 1750–1758
  * Gaißer et al. (2010) Gaißer S, Ruppert M, Schmid F (2010) A multivariate version of Hoeffding’s phi-square. Journal of Multivariate Analysis 101(10):2571–2586
  * Ghosh et al. (2006) Ghosh AK, Chaudhuri P, Sengupta D (2006) Classification using kernel density estimates: multiscale analysis and visualization. Technometrics 48(1):120–132
  * Gieser and Randles (1997) Gieser PW, Randles RH (1997) A nonparametric test of independence between two vectors. Journal of the American Statistical Association 92(438):561–567
  * Gretton et al. (2007) Gretton A, Fukumizu K, Teo CH, Song L, Schölkopf B, Smola A (2007) A kernel statistical test of independence. In: Advances in Neural Information Processing Systems, pp 585–592
  * Gretton et al. (2012) Gretton A, Borgwardt KM, Rasch MJ, Schölkopf B, Smola A (2012) A kernel two-sample test. Journal of Machine Learning Research 13(1):723–773
  * Heller et al. (2012) Heller R, Gorfine M, Heller Y (2012) A class of multivariate distribution-free tests of independence based on graphs. Journal of Statistical Planning and Inference 142(12):3097–3106
  * Heller et al. (2013) Heller R, Heller Y, Gorfine M (2013) A consistent multivariate test of association based on ranks of distances. Biometrika 100(2):503–510
  * Hoeffding (1948) Hoeffding W (1948) A non-parametric test of independence. The Annals of Mathematical Statistics 19(4):546–557
  * Kendall (1938) Kendall MG (1938) A new measure of rank correlation. Biometrika 30(1–2):81–93
  * Kibble (1945) Kibble WF (1945) An extension of a theorem of Mehler’s on Hermite polynomials. In: Mathematical Proceedings of the Cambridge Philosophical Society, Cambridge Univ Press, vol 41, pp 12–15
  * Massart (1990) Massart P (1990) The tight constant in the Dvoretzky-Kiefer-Wolfowitz inequality. The Annals of Probability 18(3):1269–1283
  * McDiarmid (1989) McDiarmid C (1989) On the method of bounded differences. Surveys in Combinatorics 141(1):148–188
  * Nelsen (1996) Nelsen RB (1996) Nonparametric measures of multivariate association. In: Rüschendorf L, Schweizer B, Taylor MD (eds) Distributions with Fixed Marginals and Related Topics, Lecture Notes–Monograph Series, vol 28, Institute of Mathematical Statistics, Hayward, C A, pp 223–232
  * Nelsen (2002) Nelsen RB (2002) Concordance and copulas: a survey. In: Cuadras CM, Fortiana J, Rodriguez-Lallena JA (eds) Distributions with Given Marginals and Statistical Modelling, Springer Netherlands, Dordrecht, pp 169–177
  * Nelsen (2013) Nelsen RB (2013) An Introduction to Copulas. Springer
  * Newton (2009) Newton MA (2009) Introducing the discussion paper by Székely and Rizzo. The Annals of Applied Statistics 3(4):1233–1235
  * Pfister and Peters (2016) Pfister N, Peters J (2016) dHSIC: Independence testing via Hilbert Schmidt independence criterion. URL https://CRAN.R-project.org/package=dHSIC, r package version 1.0
  * Pfister et al. (2017) Pfister N, Bühlmann P, Schölkopf B, Peters J (2017) Kernel-based tests for joint independence. Journal of the Royal Statistical Society: Series B 80(1):5–31
  * Póczos et al. (2012) Póczos B, Ghahramani Z, Schneider J (2012) Copula-based kernel dependency measures. In: Langford J, Pineau J (eds) Proceedings of the 29th International Conference on Machine Learning, New York, pp 775–782
  * Reshef et al. (2011) Reshef DN, Reshef YA, Finucane HK, Grossman SR, McVean G, Turnbaugh PJ, Lander ES, Mitzenmacher M, Sabeti PC (2011) Detecting novel associations in large data sets. Science 334(6062):1518–1524
  * Sarkar and Ghosh (2018) Sarkar S, Ghosh AK (2018) Some multivariate tests of independence based on ranks of nearest neighbors. Technometrics 60(1):101–111
  * Schmid et al. (2010) Schmid F, Schmidt R, Blumentritt T, Gaißer S, Ruppert M (2010) Copula-based measures of multivariate association. In: Jaworski P, Durante F, Hardle WK, Rychlik T (eds) Copula Theory and Its Applications, Springer, pp 209–236
  * Schweizer et al. (1981) Schweizer B, Wolff EF, et al. (1981) On nonparametric measures of dependence for random variables. The Annals of Statistics 9(4):879–885
  * Spearman (1904) Spearman C (1904) The proof and measurement of association between two things. The American Journal of Psychology 15(1):72–101
  * Sriperumbudur et al. (2010) Sriperumbudur BK, Gretton A, Fukumizu K, Schölkopf B, Lanckriet GRG (2010) Hilbert space embeddings and metrics on probability measures. Journal of Machine Learning Research 11:1517–1561
  * Székely et al. (2007) Székely GJ, Rizzo ML, Bakirov NK (2007) Measuring and testing dependence by correlation of distances. The Annals of Statistics 35(6):2769–2794
  * Taskinen et al. (2003) Taskinen S, Kankainen A, Oja H (2003) Sign test of independence between two random vectors. Statistics and Probability Letters 62(1):9–21
  * Taskinen et al. (2005) Taskinen S, Oja H, Randles RH (2005) Multivariate nonparametric tests of independence. Journal of the American Statistical Association 100(471):916–925
  * Tsukahara (2005) Tsukahara H (2005) Semiparametric estimation in copula models. Canadian Journal of Statistics 33(3):357–375
  * Úbeda-Flores (2005) Úbeda-Flores M (2005) Multivariate versions of Blomqvist’s beta and Spearman’s footrule. Annals of the Institute of Statistical Mathematics 57:781–788
  * Um and Randles (2001) Um Y, Randles RH (2001) A multivariate nonparametric test of independence among many vectors. Journal of Nonparametric Statistics 13(5):699–708