arXiv:1706.05940v2 [math.ST] 22 Nov 2017

Detection of Block-Exchangeable Structure
in Large-Scale Correlation Matrices
Samuel Perreault ∗
Thierry Duchesne
Département de mathématiques et de statistique, Université Laval
Johanna G. Nešlehová
Department of Mathematics and Statistics, McGill University
November 23, 2017
Correspondence:
Samuel Perreault [samuel.perreault.3@ulaval.ca]
Département de mathématiques et de statistique
1045, avenue de la Médecine
Université Laval
Québec, QC, G1V 0A6
CANADA

∗

This work has been funded by individual operating grants from the Natural sciences and engineering
research council of Canada to TD (RGPIN-2016-05883) and JN (RGPIN-2015-06801), team grants from
the Fonds de recherche du Québec – Nature et technologies to TD and JN (2015–PR–183236) and from the
Canadian Statistical Sciences Institutes to JGN and a graduate scholarship from the Fonds de recherche
du Québec – Nature et technologie to SP.

1

Abstract
Correlation matrices are omnipresent in multivariate data analysis. When the
number d of variables is large, the sample estimates of correlation matrices are typically noisy and conceal underlying dependence patterns. We consider the case when
the variables can be grouped into K clusters with exchangeable dependence; an assumption often made in applications in finance and econometrics. Under this partial
exchangeability condition, the corresponding correlation matrix has a block structure and the number of unknown parameters is reduced from d(d − 1)/2 to at most
K(K + 1)/2. We propose a robust algorithm based on Kendall’s rank correlation
to identify the clusters without assuming the knowledge of K a priori or anything
about the margins except continuity. The corresponding block-structured estimator
performs considerably better than the sample Kendall rank correlation matrix when
K < d. Even in the unstructured case K = d, though there is no gain asymptotically,
the new estimator can be much more efficient in finite samples. When the data are
elliptical, the results extend to linear correlation matrices and their inverses. The
procedure is illustrated on financial stock returns.

Keywords: Agglomerative clustering; Constrained maximum likelihood; Copula; Kendall’s
tau; Parameter clustering; Shrinkage

2

1

INTRODUCTION

Relationships between the components of a random vector X = (X1 , . . . , Xd ), say, are
of prime interest in many fields where statistical methods are used. Traditionally, this
dependence is summarized through a correlation matrix. When X is multivariate Normal,
the classical choice is the linear correlation matrix, which is omnipresent in inference. When
multivariate Normality fails, however, as it frequently does in applications in econometrics,
finance, insurance, and risk management, linear correlation can be grossly misleading and
may not even exist (Embrechts et al., 2002). For this reason, it is safer to use a rank
correlation matrix, such as the matrix of pair-wise Kendall’s taus or Spearman’s rhos.
In high dimensions, empirical correlation matrices typically conceal underlying dependence patterns. This is due to their sheer size, but also to the inherent imprecision of
the estimates, especially when the sample size is small compared to dimension d. As an
example, which we detail in Section 7, consider the log-returns of 106 stocks included in the
NASDAQ100 index from the beginning of 2017 up to 2017/09/16, giving 177 observations.
The empirical Kendall rank correlation matrix based on residuals from a fitted stochastic
volatility model is displayed in the left panel of Figure 1. The darker the color, the larger
the correlation. Hardly any dependence pattern is visible from this matrix.
Noisiness of sample correlation matrices is well documented. Several strategies have
been proposed to remedy for it, most notably shrinkage (Ledoit and Wolf, 2003a,b, 2004;
Schäfer and Strimmer, 2005). Alternative procedures developed in the context of graphical
models consist of decomposing a noisy inverse covariance matrix into a low-rank matrix
and a sparse matrix (Chandrasekaran et al., 2010; Ma et al., 2013; Agarwal et al., 2012).
We follow a different path in this article. Motivated by the above NASDAQ example
and many others, we focus on applications in which it makes sense to assume that the

3

Figure 1: The matrix of pair-wise empirical Kendall’s tau corresponding to the original
labeling of the d = 106 stocks included in the NASDAQ100 index (left panel) and after relabeling (middle panel). The right panel shows the improved estimate obtained by applying
Algorithm 1 to the stocks data and selecting a cluster structure with α = .95.
correlation matrix has a block structure. This arises when the variables can be grouped
into K disjoint clusters in such a way that for any two clusters A and B, and any Xi ∈ A and
Xj ∈ B, the correlation between Xi and Xj satisfies ρ(Xi , Xj ) = ρAB . In other words, all
variables within each cluster are equicorrelated and the between-cluster correlation depends
only on the clusters but not the representatives. Assuming a block structure is a way
to reduce dimensionality, for the number of unknown pairwise correlations decreases from
d(d−1)/2 to at most K(K +1)/2. Prime examples are portfolio and credit risk modelling or
more generally econometrics and finance, where variables can be grouped, e.g., by industry
classifications or risk types; a model that makes explicit use of correlation matrices with a
block structure is the block DECO model (Engle and Kelly, 2012). Correlation matrices
with a block structure also arise in the modelling of categorical, clustered, or gene expression
data, and in genomics studies. In the NASDAQ100 example, a block structure does indeed
emerge upon relabeling of the variables, as shown in the middle panel of Figure 1. The
4

picture is still very noisy, however.
This article describes a technique for learning the cluster structure from data and shows
how to use the latter to devise a more efficient estimator of the correlation matrix. No
prior knowledge of the clusters, their number or composition is assumed. We only require
that the underlying cluster structure satisfies a partial exchangeability condition, i.e., that
the dependence within each cluster is exchangeable. To identify the clusters from the
sample correlation matrix, we then devise an iterative algorithm that is similar to, but
different from, agglomerative clustering. The knowledge of K is not assumed a priori.
Though at first one might find this problem similar to the one tackled by model-based
clustering, they are different. The latter attempts to cluster together observations that
come from the same subpopulations of a multivariate mixture distribution, while the current
proposal aims at identifying the elements of a correlation matrix that are equal. The partial
exchangeability assumption imposes a particular set of constraints on the rows and columns
of the correlation matrix, thus the need for a novel approach. The algorithm also outputs
an improved estimate of the correlation matrix which has a block structure, and an estimate
of its asymptotic covariance matrix. In the above example of stock returns, the relabeling
in the middle panel was done using the clusters identified through the proposed algorithm;
the improved estimate of the correlation matrix is displayed in the right panel. As we
prove asymptotically and illustrate via simulations, the improvement of the estimator can
be substantial, in particular for K  d. Even in the unstructured case when K = d and
there is no gain asymptotically, the new estimator can perform substantially better in finite
samples due to a bias-variance tradeoff, particularly when n is small compared to d.
All procedures developed in this paper are based on the matrix T of pair-wise Kendall
rank correlations, which turned out to be slightly more convenient than Spearman’s rank
correlation matrix, for reasons stated in Section 2. A clear advantage of this approach
5

over the linear correlation matrix is that Kendall’s tau is margin-free, well-defined and
well-behaved irrespective of the distribution of X. This makes our methodology entirely
nonparameteric and margin-free. No Normality assumption is required; the only assumption made throughout the paper is that the marginal distributions of X1 , . . . , Xd are continuous. As we explain in Section 2, the partial exchangeability assumption only concerns
the dependence in X and not its univariate marginals. In particular, it is not assumed that
the variables in the same cluster are equally distributed. When X is multivariate Normal,
Student t, generalized hyperbolic or more generally elliptical, there is a one-to-one relationship between T and the linear correlation matrix, provided the latter exists (Lindskog
et al., 2002; Hult and Lindskog, 2002). This means that the linear correlation matrix or
its inverse, often termed the precision matrix, may be obtained from an estimator of T;
this idea has been exploited, e.g., in nonparanormal graphical models (Liu et al., 2012; Xue
and Zou, 2012) and copula regression (Cai and Zhang, 2017; Zhao and Genest, 2017). The
improved estimator of T developed in this paper may thus be used to obtain more efficient
estimators of the linear correlation matrix as well as the precision matrix.
Beyond the estimation of correlation itself, our procedure can be used as a first step in
more complex dependence model constructions. When d is large, a model for the distribution of X needs to be both flexible and parsimonious. Within the Normal or more generally
elliptical model, this means that the number of free parameters in the correlation matrix
needs to be reduced, and the block structure identified through our algorithm can serve
precisely this purpose. Outside the Normal model, dependence in X can be conveniently
described through copulas. Due to the result of Sklar (1959), the joint distribution of X
can be rewritten, for all x1 , . . . , xd ∈ R, as
P [X1 ≤ x1 , . . . , Xd ≤ xd ] = C (F1 (x1 ), . . . , Fd (xd )) ,

(1)

where F1 , . . . , Fd are the univariate marginal distributions of X and C is a copula, i.e., a
6

joint distribution function with standard uniform marginals. By combining an arbitrary
copula with arbitrary univariate marginals, Sklar’s result can be used to construct a wide
variety of distributions (Genest and Nešlehová, 2012; Joe, 2015). To achieve flexibility and
parsimony when d is large, dimensionality reduction needs to take place through a clever
construction of C; examples are vines (Kurowicka and Joe, 2011), factor models (Krupskii
and Joe, 2013; Hua and Joe, 2017), or hierarchical constructions (Mai and Scherer, 2012;
Brechmann, 2014). The cluster algorithm proposed in this paper is particularly well suited
for such approaches: Equicorrelated clusters can first be identified through it and modeled
by exchangeable lower-dimensional copulas. Dependence between clusters can then be
achieved subsequently through vines or factors.
The article is organized as follows. Section 2 specifies the partial exchangeability assumption and discusses its implications. In Section 3, we construct an improved estimator
of T when the cluster structure is known a priori, derive its asymptotic distribution, and
show that it is superior to the empirical Kendall rank correlation matrix when K < d
in terms of asymptotic variance. The way this estimator is constructed is then used in
Section 4 to derive an algorithm through which K and the cluster structure can be learned
from data. The algorithm also returns an improved estimate of T and an estimate of its
finite-sample covariance matrix. In Section 5 we discuss the special case when X is Normal
or elliptical, and explain how the clustering algorithm and the improved estimate of T can
be used to estimate the linear correlation matrix. The performance of the algorithm and of
the improved estimate of T is studied through simulations in Section 6 and illustrated on
NASDAQ100 stock returns in Section 7. Section 8 concludes the paper. Technical results
and proofs are relegated to Appendices A and B. Additional simulation results and details
about the clusters obtained for the application on stock returns are given in Appendices C
and D, respectively. All appendices are available in the Online Supplement, as well as the
7

R-code to implement the clustering algorithm.

2

PARTIAL EXCHANGEABILITY ASSUMPTION

Throughout, let X = (X1 , . . . , Xd ) be a random vector with continuous univariate marginals,
denoted F1 , . . . , Fd . In this case, the copula C in Sklar’s decomposition (1) is unique; in
fact, it is the joint distribution of the vector (F1 (X1 ), . . . , Fd (Xd )). The following partial
exchangeability assumption plays a central role in this paper.
Partial Exchangeability Assumption (PEA). For j = 1, . . . , d, let Uj = Fj (Xj ). A
partition G := {G1 , . . . , GK } of {1, . . . , d} satisfies the Partial Exchangeability Assumption
(PEA) if for any u1 , . . . , ud ∈ [0, 1] and any permutation π of 1, . . . , d such that for all
j ∈ {1, . . . , d} and all k ∈ {1, . . . , K}, j ∈ Gk if and only if π(j) ∈ Gk , one has
C(u1 , . . . , ud ) = C(uπ(1) , . . . , uπ(d) )
L

L

or equivalently, (U1 , . . . , Ud ) = (Uπ(1) , . . . , Uπ(d) ), where = denotes equality in distribution.
To understand the PEA, note first that the partition G = {{1, . . . , d}} with K = 1
satisfies the PEA only if C is fully exchangeable, meaning that for all u1 , . . . , ud ∈ [0, 1]
and any permutation π of 1, . . . , d, C(u1 , . . . , ud ) = C(uπ(1) , . . . , uπ(d) ); examples of fully
exchangeable copulas are Gaussian or Student t with an equicorrelation matrix, and all
Archimedean copulas. When K > 1, the PEA is a weaker version of full exchangeability.
A partition G for which PEA holds divides X1 , . . . , Xd into clusters such that the copula
C is invariant under within-cluster permutations. In particular, for any k ∈ {1, . . . , K},
the copula of (Xj , j ∈ Gk ) is fully exchangeable. In contrast to many standard clustering
contexts, the PEA does not imply that variables in the same cluster are equally distributed,
8

because it only concerns the copula C of X and not the marginals F1 , . . . , Fd . Also note
that the partition for which the PEA holds may not be unique; we discuss this point in
greater detail in Section 4. Yet, the coarsest partition G that satisfies the PEA is unique,
viz.
G=

arg min
G∗

(|G ∗ |).

(2)

satisfies the PEA

We should also emphasize that the PEA is not restrictive in any way. In the completely
unstructured case, G in (2) is {{1}, . . . , {d}}. However, as we demonstrate in Section 6,
there may still be a substantial advantage in considering a coarser partition for inference
purposes in finite samples. When the partition G in (2) is such that |G| = K < d, a dimension reduction from d to K takes place. Several commonly used models make an implicit
or explicit use of this kind of dimension reduction; examples include latent variable models
(e.g., frailty or random effects models), Markov random fields or graphical models, hierarchical copula models (Brechmann, 2014; Mai and Scherer, 2012), factor copulas (Gregory
and Laurent, 2004; Krupskii and Joe, 2015), or nested Archimedean models (Joe, 2015).
Definition 1. For a partition G that satisfies the PEA, we write Xi ∼ Xj whenever i, j ∈ Gk
for some k ∈ {1, . . . , K}. Furthermore, the cluster membership matrix ∆ is a d × d matrix
whose (i, j)th element is given, for all i, j ∈ {1, . . . , d}, by ∆ij = 1(Xi ∼ Xj ).
Next, let T be the d × d matrix of pairwise Kendall correlation coefficients. Specifically,
for all i, j ∈ {1, . . . , d}, the (i, j)th element of T is Tij = τ (Xi , Xj ), where
τ (Xi , Xj ) = Pr((Xi − Xi∗ )(Xj − Xj∗ ) > 0) − Pr((Xi − Xi∗ )(Xj − Xj∗ ) < 0)
is the population version of Kendall’s tau between Xi and Xj , i.e. the difference between
the probabilities of concordance and discordance of (Xi , Xj ) and its independent copy
(Xi∗ , Xj∗ ). Because τ (Xi , Xj ) depends only on the copula Cij of (Xi , Xj ), viz.
Z
τ (Xi , Xj ) = −1 + 4 Cij (ui , uj ) dCij (ui , uj ),
9

(3)

see, e.g., Nelsen (2006), it is not surprising that under the PEA, several entries in T are
identical. This is specified in the next result, which follows directly from the PEA and (3).
Proposition 1. Suppose that the partition G of {1, . . . , d} satisfies the PEA and that
Xi1 ∼ Xi2 and Xj1 ∼ Xj2 where i1 6= j1 and i2 6= j2 . Then the copulas Ci1 j1 and Ci2 j2 are
identical and, consequently, Ti1 j1 = Ti2 j2 .
Remark 1. It follows from Proposition 1 that whenever Xi1 ∼ Xi2 and Xj1 ∼ Xj2 where
i1 6= j1 and i2 6= j2 , any copula-based measure of association κ will satisfy κ(Xi1 , Xj1 ) =
κ(Xi2 , Xj2 ); examples of such measures include Spearman’s rho or Gini’s gamma (Nelsen,
2006). The reason why we focus on Kendall’s tau in this paper is that the asymptotic
and finite-sample variance of the empirical Kendall’s tau matrix have a tractable form, for
example much more so than that of Spearman’s rho (Borkowf, 2002). Having said that, the
procedures proposed here could in principle be extended to other measures of association.
Suppose now that a partition G with |G| > 1 satisfies the PEA. Proposition 1 then
implies that when Xi ∼ Xj for some i 6= j, the ith and jth rows and columns in T are
identical, once the diagonal entries are aligned. Consequently, if the variables are relabeled
so that the clusters are contiguous, then the cluster membership matrix ∆ is block-diagonal
and T becomes a block matrix. As an example, consider the following simple scenario,
which is used throughout the paper to illustrate the main concepts.
Example 1. Consider a random vector X∗ of dimension d = 10 such that the partition
G ∗ = {G1∗ , G2∗ , G3∗ } of size K = 3 given by G1∗ = {1, 3, 6, 9}, G2∗ = {5, 7, 8} and G3∗ = {2, 4, 10}
satisfies the PEA. The corresponding cluster membership matrix ∆∗ and the matrix T∗ of
pairwise Kendall correlations are shown in Figure 2 (the exact distribution of X∗ does not
matter at this point). In this case, the clusters are not contiguous, i.e., ∆∗ is not block
diagonal, and the block structure of T∗ is not easily seen. Once the variables are relabeled as
10

1.0

0.5

0.0

∆∗

T∗

∆

T

Figure 2: The Kendall correlation and cluster membership matrices corresponding to the
vector X∗ in Example 1 before (T∗ and ∆∗ ) and after (T and ∆) relabeling of the variables.
∗
X = (X1 , . . . , X10 ) = (X1∗ , X3∗ , X6∗ , X9∗ , X5∗ , X7∗ , X8∗ , X2∗ , X4∗ , X10
), the partition that satisfies

the PEA becomes G := {G1 , G2 , G3 }, where
G1 := {1, 2, 3, 4},

G2 := {5, 6, 7},

G3 := {8, 9, 10}.

(4)

The clusters are now contiguous, ∆ is block-diagonal and T has an apparent block structure,
viz. Figure 2. Every time we revisit this example, we work with the relabeled vector X.
Although we use examples in which the matrix ∆ is block-diagonal for illustrative
purposes, contiguity of the clusters is not required. In fact, given that G is unknown, the
variables are unlikely to be labeled so that T has an apparent block structure. To describe
the latter, we first need additional notation. To this end, let R be an arbitrary symmetric
d × d matrix. The entries above the main diagonal can be stacked in a vector, say ρ, of
length p = d(d − 1)/2. Note that the diagonal elements of R play no role at this point. The
particular way the vectorization is done is irrelevant, as long as it is the same throughout.
For example, one may use the lexicographical ordering viz.
ρ = (R12 , . . . , R1d , R23 , . . . , R2d , . . . , R(d−1)d )> .
11

(5)

For arbitrary r ∈ {1, . . . , p}, (ir , jr ) refers to the pair of indices ir < jr such that ρr = Rir jr .
Now any partition G = {G1 , . . . , GK } of {1, . . . , d} induces a partition of the elements of ρ,
or, equivalently, of {1, . . . , p}. For any k1 ≤ k2 ∈ {1, . . . , K}, let
Bk1 k2 = {r ∈ {1, . . . , p} : (ir , jr ) ∈ (Gk1 × Gk2 ) ∪ (Gk2 × Gk1 )}.

(6)

Note that the total number of nonempty blocks Bk1 k2 is
L = K(K − 1)/2 +

K
X

1(|Gk | > 1)

(7)

i=1

because when k1 = k2 , Bk1 k2 is nonempty only if |Gk1 | > 1. Refering to the sets Bk1 k2 using
a single index, the partition of {1, . . . , p} is then given by
BG := {B` : 1 ≤ ` ≤ L}.

(8)

In analogy to the cluster membership matrix ∆, we define a p×L block membership matrix
B; for all r ∈ {1, . . . , p} and ` ∈ {1, . . . , L}, its (r, `)th element is given by
Br` = 1(r ∈ B` ).

(9)

Finally, define the set TG of all symmetric matrices with a block structure given by BG , viz.
TG := {R ∈ Rd×d : R symmetric and ∀` ∈ {1, . . . , L} r, s ∈ B` ⇒ Rir jr = Ris js }.

(10)

Note that only the elements of R that are above the main diagonal enter the definition of
TG in (10), so that the diagonal elements of R play no role.
Now suppose that G is a partition of {1, . . . , d} such that the PEA holds and that
the elements above the main diagonal of T are stacked in τ . By Proposition 1, for any
` ∈ {1, . . . , L} and r, s ∈ B` , τ r = τ s , or, equivalently, Tir jr = Tis js . This means that
12

T ∈ TG ; when no confusion can arise, we will also write τ ∈ TG . Consequently, there are
only L distinct elements in τ . Storing these in a vector τ ∗ ∈ [−1, 1]L , we thus have
τ = Bτ ∗ .

(11)

This means that when PEA holds, the number of free parameters in T is reduced from
d(d − 1)/2 to L given in (7). We revisit Example 1 to illustrate.
Example 2. Consider the matrix T corresponding to X in Example 1, and stack it in a
vector τ of length p = 45 constructed as in (5). Given that there are K = 3 clusters given
in (4), L = 6. Consequently, τ ∗ is of length 6, i.e., the cluster structure G reduces the
number of free parameters in T from 45 to 6. The 6 distinct blocks can be seen in Figure 2
or more clearly in the left panel of Figure 10 in the Online Supplement.

3

IMPROVED ESTIMATION OF T

Suppose that Xi = (Xi1 , . . . , Xid )> , i = 1, . . . , n, is a random sample from X. The classical
nonparametric estimator of T is T̂; for i, j ∈ {1, . . . , d}, its (i, j)th element is given by
T̂ij := −1 +

X
4
1(Xri ≤ Xsi )1(Xrj ≤ Xsj ).
n(n − 1) r6=s

(12)

As we explained in Section 2, if the PEA holds for some partition G with |G| < d, T has a
block structure and the number of free parameters reduces from d(d − 1)/2 to L. In this
section, we show that an a priori knowledge of G leads to a more efficient estimator of T.
Recall first that for all i 6= j ∈ {1, . . . , d}, T̂ij is a U -statistic and thus unbiased
and asymptotically Normal (Hoeffding, 1947, 1948). The behavior of T̂ was studied by
El Maache and Lepage (2003) and Genest et al. (2011); results pertaining to the closely related coefficient of agreement appear in Ehrenberg (1952). If τ and τ̂ denote the vectorized
13

versions of T and T̂ respectively, one has, as n → ∞,
√
where

n(τ̂ − τ )

N (0p , Σ∞ ),

(13)

denotes convergence in distribution and 0p is the p-dimensional vector of zeros.

Genest et al. (2011) provide expressions for the asymptotic variance Σ∞ as well as the
finite sample variance Σ of τ̂ ; the latter is also given in Part A of the Online Supplement.
The asymptotic Normality of τ̂ specified in (13) suggests using the following loss function ` : [−1, 1]d → [0, ∞) as a basis for inference:
`(t|τ̂ , Σ) := (τ̂ − t)> Σ−1 (τ̂ − t).

(14)

This loss function is a (Mahalanobis) distance between t and τ̂ , accounting for the heterogeneous variability of the entries of τ̂ . The fact that the finite sample variance Σ is
unknown is irrelevant for now; it will only become a concern in Section 4.
Considering an arbitrary t ∈ [−1, 1]d , it is obvious that ` attains its minimum at τ̂
since `(t|τ̂ , Σ) ≥ 0 = `(τ̂ |τ̂ , Σ). Now suppose that G is a partition of {1, . . . , d} such that
the PEA holds. Unless |G| = d, it is extremely unlikely that T̂ has the block structure
implied by G, i.e., T̂ ∈ TG . By transforming the loss function ` in (14) using (11), we can
now introduce the structural constraints implied by G into the estimation procedure.
Theorem 1. Suppose that G is a partition of {1, . . . , d} such that the PEA holds. Then
for `(t|τ̂ , Σ) as in (14) and the block membership matrix B defined in (9),
τ̃ (τ̂ |G) := arg min `(t|τ̂ , Σ) = Γτ̂ ,

(15)

t∈TG

where Γ := BB+ and

+

denotes the Moore–Penrose generalized inverse. Furthermore, for

any ` ∈ {1, . . . , L} and r ∈ B` ,
τ̃ (τ̂ |G)r =

1 X
τ̂ s .
|B` | s∈B
`

14

Proof. Observe that any t ∈ TG can be expressed as Bt∗ for some t∗ ∈ [−1, 1]L . Solving
∂
`(Bt∗ |τ̂ , Σ)
∂t∗

= 0 for t∗ gives as the unique solution

−1 > −1
τ̃ ∗ = B> Σ−1 B
B Σ τ̂ .

From Proposition 6 and Lemma 4 in the Online Supplement, we have that B> Σ−1 =
B> Σ−1 BB+ . Consequently, τ̃ ∗ = B+ τ̂ and τ̃ = Γτ̂ given that τ̃ = Bτ̃ ∗ . The expression
for τ̃ (τ̂ |G)r follows immediately from (33) in the proof of Lemma 4.
Remark 2. The block-wise averages that appear in τ̃ are akin to the pair-wise linear correlation averaging of Elton and Gruber (1973) and Ledoit and Wolf (2003a), and particularly
the block DECO of Engle and Kelly (2012) which uses block-wise averages of linear correlations. Although τ̃ seems very natural, the PEA plays a crucial role in the proof of Theorem
1. It is needed to invoke Proposition 6 and Lemma 4 and get a Σ matrix with a structure
that will lead to the required simplifications. In particular, it is not enough to assume that
τ ∈ TG for some partition G.
When it introduces no confusion, we refer to τ̃ (τ̂ |G) as τ̃ and to its matrix version
T̃(T̂|G) thereof as T̃. What is crucial in Theorem 1 is that τ̃ consists of the cluster averages
of the elements of τ̂ and as such does not involve the unknown finite sample variance Σ
of τ̂ , so an estimator of Σ is not needed to compute τ̃ . The information contained in G is
introduced by projecting T̂ onto TG . The resulting estimator T̃ is expected to be closer to
the original matrix T because the entries that estimate a same value are averaged over, thus
reducing the estimation variance. In fact, for any r ∈ {1, . . . , p}, the asymptotic variance
of τ̃ r is less than or equal to that of τ̂ r as a result of the following theorem.
Theorem 2. Let G be a partition of {1, . . . , d} such that the PEA holds. For B given by
(9), let Γ = BB+ and τ̃ = Γτ̂ be as in Theorem 1. Then the following statements hold:
15

T

T̂

T̃

Figure 3: The matrices T, T̂ and T̃ in Example 3.
(i) As n → ∞,

√

n(τ̃ − τ )

N (0, ΓΣ∞ ).

(ii) The matrix Σ∞ − ΓΣ∞ is nonnegative definite.
Proof. Because τ̃ − τ = Γ(τ̂ − τ ), (i) follows from (13) and the fact that ΓΣ∞ Γ = ΓΣ∞ by
Lemma 5 in the Online Supplement. From Lemma 5 therein, Σ∞ − ΓΣ∞ = (Ip − Γ)Σ∞ =
(Ip − Γ)> Σ∞ (Ip − Γ), where Ip denotes the p × p identity matrix. Consequently, (ii) follows
from the fact that Σ∞ is nonnegative definite.
To conclude this section, we illustrate τ̃ using the setup in Example 1.
Example 3. Consider a random sample of size n = 70 from the vector X in Example 1; we
used X to be Normally distributed with Kendall correlation matrix T displayed in Figure 2.
Figure 3 displays T, T̂ and T̃. For this one simulated sample, it is clear that T̃ is a better
estimate of T than T̂.

16

LEARNING THE STRUCTURE G

4

Because the cluster structure G is typically unknown, the improved estimator τ̃ derived in
the previous section cannot be directly used. In this section, we propose a way to learn G
from data and to obtain an improved estimator of τ as a by-product. The task of learning
G is split into two subtasks: we first identify d candidate structures in Section 4.1 and then
choose one among them in Section 4.2.

4.1

Creating a set of candidate structures

The first observation to be made is that the cluster structure G for which the PEA holds
may not be unique. This is because if the PEA holds for some cluster structure G, it holds
for any refinement G ∗ thereof defined as follows.
Definition 2. Let G = {G1 , . . . , GK } be a partition of {1, . . . , d}. A refinement of G is a
partition G ∗ = {G1 , . . . , GK ∗ } of {1, . . . , d} such that K ∗ > K and
∀k ∗ ∈ {1, . . . , K ∗ } ∃k ∈ {1, . . . , K}

such that

Gk∗ ⊆ Gk .

The block structure implied by a refinement of G is consistent with the the block structure implied by G. This is formalized in the next proposition, which follows easily from (10).
Proposition 2. For two partitions G and G ∗ of {1, . . . , d}, G ∗ is a refinement of G if and
only if TG ⊆ TG ∗ .
Proposition 2 implies in particular that if T ∈ TG for some partition G, then for any
refinement G ∗ thereof, T ∈ TG ∗ . This is illustrated in Example 4 below.
Example 4. Consider the partition G given in (4) in Example 1. The partition G ∗ =
{G1∗ , . . . , G5∗ } with G1∗ = {1, 2}, G2∗ = {3, 4}, G3∗ = {5, 6, 7}, G4∗ = {8, 9, 10} is a refinement
17

∆†

T̃(T̂|G ∗ )

∆

T̃(T̂|G)

Figure 4: The estimates T̃(T̂|G ∗ ) and T̃(T̂|G) from Example 4, along with the cluster
membership matrices ∆ and ∆∗ of the partitions G and G ∗ , respectively.
of G since G1∗ , G2∗ ⊆ G1 , G3∗ ⊆ G2 and G4∗ ⊆ G3 . Consequently, G ∗ satisfies the PEA as well.
Figure 4 shows the cluster membership matrices ∆ and ∆∗ corresponding to G and G ∗ ,
respectively. Also displayed are the estimates T̃(T̂|G ∗ ) and T̃(T̂|G); one can see that the
block structure of the former is embedded in the latter but not conversely.
While the partition for which the PEA holds may not be unique, the coarsest partition
G in (2) is. The fact that any refinement of G in (2) also satisfies the PEA motivates
us to start with the finest possible partition G (d) := {{1}, . . . , {d}} for which the PEA
always holds, and to merge the clusters one at a time in a way that resembles hierarchical
agglomerative clustering. Specifically, we will create a path P = {G (d) , . . . , G (1) } through
the set of all possible partitions of {1, . . . , d} with |G (i) | = i for i = 1, . . . , d, with the aim
that G given in (2) is an element of P. The construction of P is motivated by the following
observation.
Proposition 3. Let G be an arbitrary partition of {1, . . . , d}, B the associated block membership matrix (9) and Γ = BB+ , where B+ is the generalized Moore–Penrose inverse of
−1

B. Let also τ̃ = Γτ̂ . If Σ∞ is positive definite, and Σ̂ is an estimator of Σ such that Σ̂
exists and, as n → ∞, nΣ̂ → Σ∞ element-wise in probability, the following holds.
18

−1

(i) If G fulfils the PEA, (τ̂ − τ̃ )> (Σ̂ /n)(τ̂ − τ̃ ) → 0 in probability.
−1

(ii) If G does not fulfill the PEA, (τ̂ − τ̃ )> (Σ̂ /n)(τ̂ − τ̃ ) → τ > (Ip − Γ)Σ−1
∞ (Ip − Γ)τ
in probability; if Γτ 6= τ , the limit is strictly positive.
−1
Proof. First, note that as n → ∞, Σ−1 /n → Σ−1
is a
∞ in probability given that A 7→ A

continuous map for nonsingular matrices (Stewart, 1969). Now write
−1

−1

(τ̂ − τ̃ )> (Σ̂ /n)(τ̂ − τ̃ ) = τ̂ > (Ip − Γ)(Σ̂ /n)(Ip − Γ)τ̂ .
−1

Because as n → ∞, τ̂ → τ in probability by (13), τ̂ > (Ip − Γ)(Σ̂ /n)(Ip − Γ)τ̂ converges
to τ > (Ip − Γ)Σ−1
∞ (Ip − Γ)τ in probability, proving the statement (ii). When G fulfills the
PEA, Γτ = τ and hence (Ip − Γ)τ = 0p , proving (i).
Remark 3. Note that in Proposition 3 (ii), Γτ = τ could indeed occur even if G does not
satisfy the PEA. This is because Kendall’s τ of two distinct copulas may be equal. However
the PEA must be met in order for Theorems 1 and 2 to hold.
Motivated by Proposition 3, the construction of P relies on slowly introducing information through constraints under which the loss function ` in (14) is minimized. Before
proceeding with this idea however, the estimation of the unknown finite-sample variance Σ
of τ̂ needs to be considered. While Σ does not appear in the estimator τ̃ in Theorem 1, it
is relevant for the construction of P. In Appendix A in the Online Supplement, we explain
how to obtain an estimator of Σ for a given partition G. Because Σ inherits a certain block
structure if the PEA holds for G, this estimator is a function of an empirical estimator Σ̂ of
Σ, the structure G and a shrinkage parameter w; hence we denote it by Σ̃(Σ̂ |τ̂ , G, w). In
Section A.3 of the Online Supplement it is shown that nΣ̃(Σ̂ |τ̂ , G, w) → Σ∞ element-wise
in probability if w shrinks to zero as n → ∞.
19

(i)

Now suppose that the ith partition G (i) has been selected; let Σ̃w := Σ̃(Σ̂ |τ̂ , G (i) , w)
denote the corresponding estimate of Σ. To select the (i − 1)-st cluster structure G (i−1) ,
merge two clusters at a time and choose the optimal merger, in the sense that
G (i−1) :=

arg min
G ∗ :TG ∗ ⊂TG (i) , |G ∗ |=i−1



(i)
` τ̃ (τ̂ |G ∗ ) |τ̂ , Σ̃w .

(16)

The minimization in (16) is done by simply going through all i(i − 1)/2 possible mergers;
TG ∗ ⊂ TG (i) indicates that G (i) must be a refinement of G ∗ , so that the previously introduced
equality constraints are carried. We then update the estimate of τ to τ̃ (τ̂ |G (i−1) ) as in
(i−1)

Theorem 1, the estimate of Σ to Σ̃w

= Σ̃(Σ̂ |τ̂ , G (i−1) , w) and iterate the above steps

until i = 1. The whole procedure is formalized in Algorithm 1 below.
Algorithm 1.
(0) Initialization. Fix w ∈ [0, 1] and set
G (d) := {{1}, . . . , {d}} and

(d)

Σ̃w := Σ̃(Σ̂|τ̂ , G (d) , w).

(1) Iteration. For i ∈ {d − 1, . . . , 1},
(1a) Structure selection. Set G (i) =

arg min
G ∗ :TG ∗ ⊂TG (i+1) , |G ∗ |=i



(i+1)
` τ̃ (τ̂ |G ∗ )|τ̂ , Σ̃w
;

(i)

(1b) Update. Set Σ̃w = Σ̃(Σ̂|τ̂ , G (i) , w);
(2) Output. Return P := {G (d) , . . . , G (1) }.
Algorithm 1 returns a sequence P = {G (d) , . . . , G (1) } of decreasingly complex structures;
note that G (i) is a refinement of G (i−1) for all i = 2, . . . , d. From each G (i) in P, we can then
compute the cluster membership matrix ∆(i) , the improved estimate τ̃ (i) = τ̃ (τ̂ |G (i) ) of τ
(i)

defined as in Theorem 1, its matrix version T̃ , and an estimate Σ̃
20

(i)

of Σ (upon setting

w = 0). Furthermore, we may also construct the corresponding Γ matrix, say Γ(i) , and
(i)

obtain a consistent estimator of the covariance matrix of τ̃ (i) , Γ(i) Σ̃ Γ(i) , which simplifies
to Γ(i) Σ̃

(i)

by Lemma 6.

If present on the path, the coarsest possible structure G in (2) is always G (K) , where
K = |G|. Note that T̃

(d)

(1)

= T̂ and T̃ , the equicorrelation matrix corresponding to

G (1) = {{1, . . . , d}}, are inevitable outputs of the algorithm. This is why we refer to P as
a path: in the space of all d × d matrices, it is the path we took to go from T̃

(d)

(1)

to T̃ .

Example 5. Figure 5 depicts the application of Algorithm 1 on T̂ constructed from the
random sample in Example 3 of X from Example 1. Here, the true cluster structure G in
(2) is given by (4). It indeed lies on the path; it corresponds to ∆(3) .
Remark 4. Consider a partition G for which the PEA holds, and let G ∗ be a refinement
thereof. Let B and B∗ be the block membership matrices given by (9) corresponding to
G and G ∗ , respectively, and set Γ = BB+ and Γ∗ = B∗ (B∗ )+ . Then because Σ−1 ∈
SG ⊂ SG ∗ , where SG and SG ∗ are as defined in Appendix A.2, Lemma 6 therein applies and
(Ip − Γ∗ )> ΣΓ∗ = 0. Furthermore, for τ̃ ∗ = Γ∗ τ̂ and τ̃ = Γτ̂ , Γ∗ τ̃ = τ̃ . Hence,
(τ̂ − τ̃ )> Σ−1 (τ̂ − τ̃ ) = (τ̂ − τ̃ ∗ )> Σ−1 (τ̂ − τ̃ ∗ ) + (τ̃ ∗ − τ̃ )> Σ−1 (τ̃ ∗ − τ̃ ).

(17)

Now set K = |G|. If G (i) , i = K, . . . , d is a sequence of partitions such that G (K) = G,
G (d) = {{1}, . . . {d}}, and for each i = K, . . . , d − 1, G (i+1) is a refinement of G (i) . For
all i = K, . . . , d, let B(i) be as in (9), Γ(i) = B(i) (B(i) )+ and τ̃ (i) = Γ(i) τ̂ . A successive
application of (17) then gives
(τ̂ − τ̃

(K) >

−1

) Σ (τ̂ − τ̃

(K)

)=

d
X

(τ̃ (i) − τ̃ (i−1) )> Σ−1 (τ̃ (i) − τ̃ (i−1) )

(18)

i=K+1

In particular, for any i = K, . . . , d − 1, `(τ̃ (i) |τ̂ , Σ) ≥ `(τ̃ (i+1) |τ̂ , Σ). This motivates that
in the iteration Step (1a) of Algorithm 1, only two clusters are merged at a time.
21



(10)

∆



(5)

, T̃

∆ , T̃

(10)

(5)









(9)

(9)

∆ , T̃

(4)

(4)

∆ , T̃









Figure 5: The pairs of matrices



(8)

∆ , T̃

(3)

∆ , T̃

(i)

(i)

∆ , T̃



(8)

(3)







(7)
(7)
∆ , T̃



(2)
(2)
∆ , T̃





(6)

∆ , T̃

(1)

∆ , T̃

(6)

(1)





, i = 1, . . . , d corresponding to the path P

obtained by applying Algorithm 1 in Example 5.

4.2

Structure selection

We need to identify a final estimate of G among the d structures in P. Proposition 3
suggests that the loss will increase sharply when the clustering has become too coarse. The

22

following result, whose proof can be found in Section B.2 of the Online Supplement, offers
a way to determine when this sharp increase may have occurred.
Proposition 4. Let G = {G1 , . . . , GK } be a partition of {1, . . . , d} satisfying the PEA, let
B be the corresponding block membership matrix (9), Γ = BB+ and τ̃ = Γτ̂ . If Σ∞
is positive definite, and Σ̂ is any estimator of Σ such that Σ̂

−1

exists and, as n → ∞,

nΣ̂ → Σ∞ element-wise in probability, then, as n → ∞,
−1

`(τ̃ |τ̂ , Σ̂) = (τ̂ − τ̃ )> Σ̂ (τ̂ − τ̃ )

χ2p−L ,

where L is the number of distinct blocks given in (7).
(i)

At each iteration of Algorithm 1, Σ is estimated by Σ̃w . Proposition 4 and Section
(i)

A.3 of the Online Supplement suggests using `(τ̃ (i) |τ̂ , Σ̃w ) to get a rough idea of when too
much clustering has been applied through
h
i
(i)
α(i) := P χ2p−Li > `(τ̃ (i) |τ̂ , Σ̃w ) ,

(19)

where Li is the number of blocks given by (7) corresponding to the ith partition G (i) . For
n large enough, we expect that a sharp decrease in α(i) will occur at the first i such that
TG 6⊆ TG (i) , that is, when the Γ matrix corresponding to G (i) becomes inadmissible. We
do not use the criterion (19) as a formal p-value, but rather as a tool that can help with
structure selection. A naive automated selection procedure based on (19) is used in the
simulations of Section 6, where it leads to good estimators of τ .
Example 6. Consider again the random sample of size n = 70 from X in Example 3.
We computed α(i) , i = 1, . . . , 10, given by (19) for the path obtained with Algorithm 1 in
Example 5. As can be seen in Figure 6, the gap between α(3) and α(2) strongly suggests that
the best structure is G (3) , which is indeed the true structure in this case.
23

1.0
0.8
0.6
0.4
0.2
0.0
10

7

4

1

Figure 6: The pairs (i, α(i) ), i = 10 . . . , 1, computed in Example 6.

5

ESTIMATION OF LINEAR CORRELATION

In this section, we show how the PEA can be used to obtain improved estimates of the
classical linear correlation matrix P with entries Pij = Cor(Xi , Xj ), i, j = 1, . . . , d when
X is elliptical. Recall that an absolutely continuous random vector X follows an elliptical
distribution with mean vector µ ∈ Rd , positive definite d × d dispersion matrix D and
density generator g : [0, ∞) → [0, ∞), in notation X ∼ E(µ, D, g), if its density f satisfies,
for all x ∈ Rd ,
−1/2

f (x) = |D|


g

(x − µ)0 D−1 (x − µ)
2


,

(20)

where |D| denotes the determinant of D. The equation (20) means that the level curves
of f are concentric ellipses centred at µ. For a book treatment of elliptical distributions,
see, e.g., Fang et al. (1990) or Fang and Zhang (1990). Well-known examples of elliptical
distributions are the multivariate Normal, Student t or generalized hyperbolic distributions;
for their use in finance and risk modelling, see McNeil et al. (2015).
Note that when X ∼ E(µ, D, g), µ and D are not necessarily the mean and covariance
matrix of X, respectively; the former may not even exist. However, if E(Xi2 ) < ∞ for
all i ∈ {1, . . . , d}, E(X) = µ and there exists a constant c > 0 such that Cov(X) = cD

24

(Fang and Zhang, 1990, Theorem 2.6.4). Consequently, if the linear correlation matrix
p
P of X exists, one has, for all i, j ∈ {1, . . . , d}, Pij = Dij / Dii Djj . Surprisingly, for
all i 6= j ∈ {1, . . . , d}, the correlation coefficient Pij is in one-to-one relationship with
Kendall’s correlation τ (Xi , Xj ), viz.

Pij = sin

πTij
2


,

(21)

as shown by Fang et al. (2002) and Lindskog et al. (2002). Because the map in (21) is a
bijection, (21) can be used to construct an estimator of P, given, for all i 6= j ∈ {1, . . . , d},
by P̂ij = sin(π T̂ij /2). As illustrated by Lindskog et al. (2002), the resulting estimator
P̂ can be considerably more efficient than the sample correlation matrix, especially when
the margins of X are heavy-tailed. Recently, P̂ has been employed, e.g, in the context of
nonparanormal graphical models (Liu et al., 2012; Xue and Zou, 2012) or Gaussian copula
regression (Cai and Zhang, 2017).
Now suppose that G is a partition of {1, . . . , d} so that the PEA holds. Because X
is elliptical, this is equivalent to T ∈ TG , or, in view of (21), to P ∈ TG . Because T̃
is a more efficient estimator of T by Theorem 2, the delta method implies that P̃ ∈ TG
obtained by using T̃ij in (21) is a more efficient estimator of P than P̂. Moreover, if P
is positive definite, it follows from Lemma 7 in the Online Supplement that the precision
matrix Ω = P−1 has the same block structure as P, that is, Ω ∈ TG . As an estimator of
−1

Ω, one may thus use Ω̃ = P̃

directly if the latter is positive definite; it then follows from

Lemma 7 that Ω̃ ∈ TG . Otherwise, P̃ can first be made positive definite using one of the
shrinkage methods described, e.g., in Rousseeuw and Molenberghs (1993); its inverse can
be further improved by averaging out the entries block-wise to obtain a matrix in TG .

25

6

SIMULATION STUDY

6.1

Performance of the clustering algorithm

In this section, we first investigate whether the true cluster structure is on the path P =
{G (d) , . . . , G (1) } returned by Algorithm 1 and how often does the criterion (19) select the
true cluster structure given the latter is on the path. We only consider the case of a
modest dimension d = 20 and focus on four different matrices Ti , i = 1, . . . , 4 displayed
in Figure 7 along with the corresponding cluster membership matrices ∆i , i = 1, . . . , 4.
The matrices T1 and T2 both have the same 4 clusters of 5 variables, but with entries of
different magnitudes so that the blocks are less well separated in T2 . The entries of T2
are moreover smaller in absolute value, which implies more variability in T̂. We therefore
expect its structure to be more difficult to identify. The structure of the matrices T3 and
T4 is more complex and formed by 8 clusters of varying size. The clusters corresponding
to T4 are the most difficult to identify.
For each of the four matrices Ti , i = 1, . . . , 4, we generated 500 datasets from multivariate Normal distributions with standard Normal margins and Kendall’s rank correlation
matrix Ti . This was replicated with three different sample sizes, n = 125, 250, 500. Examples of T̂ obtained with all 12 pairs (T, n) are shown in Figure 12 of the Online Supplement.
Note that since Algorithm 1 is ranked-based, different choices of the marginal distributions
have no impact on the results so long as they are continuous. However, the dependence
structure does matter. The simulations were thus repeated using samples from the Cauchy
copula. The results were similar and are given in Appendix C of the Online Supplement.
For each simulated sample, we applied Algorithm 1 and computed the indicator of the
presence of the true cluster membership matrix ∆ on the path P, viz.
ν1 := 1(∆ ∈ P).
26

1.0
0.5
0.0
−0.5
−1.0

(∆1 , T1 )

(∆3 , T3 )

(∆2 , T2 )

(∆4 , T4 )

Figure 7: The matrices Ti and ∆i , i = 1, 2, 3, 4 used in simulations in Section 6.
The shrinkage parameter w was held fixed through all d iterations of Algorithm 1. We
considered five values for w, 0, 0.25, 0.5, 0.75 and 1. The average value ν̄1 when w = 0.75 is
reported in Table 1. For all four matrices, the value of ν̄1 increases with the sample size. For
the matrices T1 , T2 and T3 with clearly distinct blocks, the path of Algorithm 1 very often
goes through the true structure. For the matrix T4 , it is no surprise that this almost never
happens, as the generated T̂ are very noisy and the clusters are hardly distinguishable.
The results for the remaining values of w may be found in Table 5 in Appendix C of the
Online Supplement. As expected, when the sample size is large (n = 500) the value of w
has minor impact. When the sample size is small or moderate (n = 125 or 250), w = 0
leads to poorer results, but otherwise the value of w has little impact.
Our next objective was to check whether α(i) , defined in (19), is a reasonable structure
identification tool. The selection was done with the following automated procedure.
27

Table 1: Average values of ν1 over 500 simulations with w = 0.75.
T1

T2

T3

T4

n

125

250

500

125

250

500

125

250

500

125

250

500

ν̄1

0.94

1.00

1.00

0.52

0.95

1.00

0.29

0.81

0.99

0.00

0.00

0.11

Algorithm 2.
(0) Fix α ∈ (0, 1) and set i• ← d.
For i ∈ {d − 1, . . . , 1}:
(i)

(1) let Σ̃w = Σ̃(Σ̂|τ̂ , G (i) , w) and Li = |BG (i) |;

i
h
(i)
(i)
(i)
2
(2) compute α := P χp−Li > ` τ̃ |τ̂ , Σ̃w ;
(3) if α(i) ≥ α, set i• ← i.
This simple rule returns i• so that G (i
α(i

•)

≥ α. The parameter

(i)
Σ̃w

•)

is the coarsest cluster structure such that

typically improves quickly after the first few steps. One

could therefore monitor its behaviour and decrease w as i gets smaller. In this simulation
study, we fix w = 0 through the whole process for simplicity. We run Algorithm 2 with
α = 0.05, 0.1, 0.25, 0.5, and look at the proportion of times Algorithm 2 identifies the true
structure properly, given that it is present on the path obtained with Algorithm 1. Results
for T1 , T2 and T3 are shown in Table 2; results for T4 are not shown because too few paths
contained the true structure, viz. Table 1. As the sample size increases, the percentage
gets closer to 1 − α, as one would expect. Hence, it is advantageous to choose a small value
of α once n is sufficiently large. When T = T1 , T = T2 , and T = T3 , α = 0.05 becomes
the best option when n = 125, n = 250, and n = 500, respectively. For smaller sample
28

Table 2: Proportion of the 500 samples for which Algorithm 2 identified the true structure
when it was present on the path returned by Algorithm 1.
T1

T2

T3

Level α|n

125

250

500

125

250

500

125

250

500

0.05

0.94

0.96

0.96

0.67

0.94

0.95

0.04

0.51

0.96

0.10

0.90

0.90

0.92

0.78

0.93

0.90

0.14

0.63

0.92

0.25

0.77

0.77

0.77

0.79

0.80

0.76

0.37

0.69

0.77

0.50

0.54

0.53

0.50

0.60

0.54

0.51

0.53

0.61

0.59

sizes, it is advantageous to choose a larger value of α. This phenomenon is best visible
when T = T3 and n = 125; the best option in this setting is α = 0.5.

6.2

Improved estimation of Kendall’s rank correlation matrix

Even though Algorithm 1 in combination with Algorithm 2 may not identify the true
coarsest cluster structure properly, the estimator T̃ that it returns is still considerably
more efficient than T̂. We explore this in this section by first calculating the statistic
ν2 := 1 − min{||τ̃ (j) − τ ||22 : G (j) ∈ P}/||τ̂ − τ ||22 ,
i.e., one minus the smallest squared error obtained among all members of the path P
returned by Algorithm 1, normalized by the squared error of τ̂ . To assess the performance
of Algorithm 2, we also calculate
ξ(α) := 1 − ||τ • − τ ||22 /||τ̂ − τ ||22 ,

(22)

where ∆• and τ • are the outcomes of Algorithm 2 with level α. Note that ν2 and ξ(α) are
always between 0 (worst) and 1 (best).

29

Table 3: Average values of ν2 and ξ(α) for α ∈ {0.05, 0.1, 0.25, 0.5} for 500 simulation runs
with w = 0.75.
T1

T2

T3

T4

Measure|n

125

250

500

125

250

500

125

250

500

125

250

500

ν̄2

0.63

0.65

0.65

0.68

0.79

0.79

0.50

0.59

0.63

0.50

0.51

0.56

ξ̄(0.05)

0.61

0.66

0.67

0.61

0.78

0.79

0.35

0.47

0.61

0.37

0.39

0.44

ξ̄(0.10)

0.61

0.65

0.66

0.63

0.78

0.78

0.38

0.50

0.61

0.39

0.41

0.46

ξ̄(0.25)

0.59

0.61

0.63

0.64

0.75

0.75

0.41

0.52

0.59

0.41

0.42

0.47

ξ̄(0.50)

0.53

0.55

0.54

0.61

0.70

0.69

0.43

0.51

0.53

0.41

0.42

0.45

The results for the simulation design from Section 6.1 are displayed in Table 3. First
note that the results for distinct matrices Ti are not directly comparable, because the
maximum error reduction that can possibly be achieved depends on Ti . By looking at the
row pertaining to ν̄2 , we see that the mean squared error of τ̂ can be cut substantially by
choosing the best structure on the path P returned by Algorithm 1 in each of the scenarios
¯
considered. The values of ξ(α)
suggest that Algorithm 2 often selects a structure that does
nearly as well in terms of error reduction. The reactivity to changes in α is low, which is a
sign that the gap in α(i) between reasonable and poor models is often large. Nonetheless,
one can detect a similar phenomenon as in Table 2: Algorithm 2 performs best when α
is small once the sample size is large enough. In summary, Table 3 illustrates that even
if Algorithms 1 and 2 do not select the true coarsest partition G that satisfies PEA, the
resulting estimator T̃ can considerably outperform T̂. This is particularly apparent in the
last block pertaining to T4 . Even though hardly any path returned by Algorithm 1 goes
through G in this case, Table 3 shows that the structure selected by Algorithms 1 and 2
leads to a substantial reduction of mean squared error.
The above findings suggest that selecting a structure simpler than the true one might
30

T5

T6

T7

T8

Figure 8: The matrices Ti , i = 5, 6, 7, 8 used in simulations in Section 6.2.
be beneficial when T̂ is extremely noisy because even if it introduces a small bias, it reduces
the variance considerably. To demonstrate this further, we challenge Algorithms 1 and 2
by considering d ∈ {50, 100} and four additional matrices T5 , . . . , T8 displayed in Figure 8.
T5 and T6 correspond to 50 variables and have a block structure with 10 and 19 blocks,
respectively. The matrices T7 and T8 are 100×100 noisy versions of Toeplitz matrices with
a constant decay in their entries as one moves away from the diagonal. This means that
T7 and T8 are unstructured, i.e., the only partition G that satisfies PEA is {{1}, . . . , {d}}.
We also consider rather small sample sizes, n ∈ {25, 50, 100}. Because better results are
achieved for higher values of the shrinkage parameter w in small samples, we set w = 1.
The results for 500 simulation runs are summarized in Table 4. They suggest that the
method still performs well in higher-dimensional settings, even when d is larger than n.
As already observed in Tables 2 and 3, it is better to choose a more conservative value
for α in small samples, although the sensitivity to α is small. The results are particularly
interesting for T7 and T8 , which do not have a block structure. However, they still possess
a certain structure that Algorithm 1 is able to capture and T̃ performs better than T̂ due
to a bias-variance tradeoff when n is small. Note that because K = d in both cases, there
is no gain asymptotically, so the improvement should decrease with the sample size. This
31

Table 4: Average values of ν2 and ξ(α) for α ∈ {0.05, 0.5, 0.95} for 500 simulation runs
with w = 1.
T5
Measure|n

T6

T7

T8

25

50

100

25

50

100

25

50

100

25

50

100

ν̄2

0.47

0.52

0.54

0.48

0.46

0.48

0.66

0.51

0.37

0.42

0.38

0.34

ξ̄(0.05)

0.30

0.31

0.34

0.40

0.31

0.29

0.60

0.47

0.33

0.31

0.20

0.11

ξ̄(0.50)

0.34

0.37

0.38

0.40

0.33

0.31

0.59

0.48

0.34

0.32

0.22

0.13

ξ̄(0.95)

0.37

0.40

0.44

0.41

0.34

0.33

0.58

0.47

0.32

0.33

0.23

0.15

is already apparent for T8 for the sample sizes considered here.

7

APPLICATION TO STOCK RETURNS

We collected the daily value at close of all d = 106 stocks included in the NASDAQ100
index from the beginning of 2017 up to 2017/09/16 (as reported by Google Finance on
2017/09/17). The information on the components of the index (company name, sectors
and industries) were taken from www.nasdaq.com. Our primary goal is not to cluster
together stocks whose returns have a similar distribution, but rather to simplify their
dependence structure by assuming partial exchangeability. Because the PEA is associated
with a particular clustering of the variables, we can further verify if the clusters correspond
to intuitive categories like business sectors.
As is the norm in joint modelling of financial time series, we compute the log returns for
series of stocks and apply our methodology to the series of residuals from a fitted stochastic
volatility model (Patton, 2006, 2012; Rémillard, 2013). The stochastic volatility model that
yielded satisfactory results in this application was the GARCH (1, 1) model. This produced
n = 177 residuals for each of the d = 106 stocks.
32

Note that even under the assumption that the stochastic volatility model is appropriate,
the residuals are not i.i.d. However, as shown in Rémillard (2017, Corollary 2), the empirical
copula process based on the residuals has the same asymptotic behavior as the empirical
copula process for i.i.d. observations. As a result, (13) continues to hold with the very
same asymptotic variance Σ∞ . Even if Σ may not be the exact finite-sample variance of
τ̂ based on the residuals, the result of Rémillard (2017) implies that (5) and (18) remain
true. Consequently, the methodology developed here can be used as is.
For computational purposes and because the ratio n/d = 1.67 is small, we applied
Algorithm 1 with w = 1. We then computed α(i) , i = 106, . . . , 1, as given by (19), again
with w = 1. An excerpt of the plot is shown in Figure 9; α(i) for i > 20 and i < 10 were
essentially equal to 1 and 0, respectively. Among the 106 candidate structures produced by
Algorithm 1, G (16) , G (15) and G (14) stand out as the most interesting, because they precede
drops in α(i) that are typical of questionable cluster mergers. As a matter of fact the
selection procedure is not very sensitive to the arbitrary choice of the value of α: any value
of α between 0.01 and 0.99 leads to one of G (16) , G (15) and G (14) among the 106 structures
that were returned by Algorithm 1. The reduction in dimension is striking; for instance if
we use G (15) we go from d(d − 1)/2 = 5565 parameters down to K(K + 1)/2 = 120.
The left panel in Figure 1 shows the empirical Kendall’s tau matrix T̂ when the variables
are labeled consecutively as in the NASDAQ100 index. Clearly, any pattern is very hard to
discern from this matrix. The middle panel of the same figure shows T̂ once the variables
have been relabeled according to the clusters in G (15) , i.e., so that the cluster membership
matrix is block-diagonal. Now a structure is beginning to emerge. The right panel displays
the improved estimator T̃

(15)

that takes the cluster structure into account. The benefit of

using Algorithm 1 is clearly visible here.
The clusters in G (15) are listed in Appendix D of the Online Supplement. It is comforting
33

1.0
0.8
0.6
0.4
0.2
0.0
20

18

16

14

12

10

Figure 9: The pairs (i, α(i) ), i = 20, . . . , 10 corresponding to the path returned by Algorithm
1 applied to the stock returns residuals.
that most of them have intuitive interpretations in terms of business sectors or known
affiliations, particularly when the intra-cluster correlation is large. For example, cluster 4
is almost exclusively composed of stocks corresponding to semiconductor-related companies
with intra-cluster correlation of 0.44; cluster 7 contains stocks from the biopharmaceutical
industry only with intra-cluster correlation of 0.29. Furthermore, strongly related stocks
are also clustered together: GOOG and GOOGL (.89); DISCA and DISCK (.88); FOX
and FOXA (.87); LBTYK and LBTYA (.85); LILA and LILAK (.76).

8

CONCLUSION

We have developed a new approach to identify a block structure within the matrix of
Kendall correlations. Aside from a mild partial exchangeability assumption, the method is
completely nonparametric and does not require any additional assumption on the joint or
marginal distributions of the variables (insofar as they are continuous). This contribution
has the potential to be useful in all areas where the dependence between the variables is
of interest. While goodness-of-fit tests can only be performed a posteriori, the proposed

34

method can serve as a model specification guide at the early stages of the data analysis.
We have formally shown that taking advantage of the block structure of the correlation
matrix can lead to improved inference on the correlation coefficients. Not only are the new
estimators consistent and asymptotically Normal, but their asymptotic variance is smaller
than that of the empirical Kendall coefficients. This is important, as the latter tend to be
extremely noisy when d gets large. But the improvements in the inference are not only
true asymptotically when the correlation structure is known. Our simulations have shown
that the new estimator has better mean squared error in finite samples even when the
correlation structure is not known a priori and has to be estimated from the data.
The method that we propose to identify the blocks in the correlation matrix is an
algorithm that performs an agglomerative clustering of the variables, not of observations.
While this may still resemble model-based clustering, it differs from it in many aspects.
Perhaps the most important difference is that the approach we use is independent of the
distribution of the variables and concerns only τ̂ , thus making the method model-free.
We have demonstrated that as the sample size gets large, this algorithm is based on a loss
function that will assign negligible loss to merges that agree with the true block structure
and large loss to merges that do not, which ensures that the agglomerative process will yield
a set of d potential structures that includes the true one. We have also developed a tool
that can help to identify the reasonable structures among the d proposed by the algorithm.
The asymptotic properties of the loss function suggest that this tool will be adequate for
large samples; our simulation study indicates good performance in finite samples as well.
This suggests that block-exchangeable structures can be used as general-purpose shrinkage
targets when n is small compared to d.
The loss function used in the agglomerative algorithm depends on a variance matrix Σ
that must be estimated. We have shown under the Partial Exchangeability Assumption that
35

this matrix and its inverse share a common structure property. We exploit this property
to improve the consistent plug-in estimator Σ̂ of Σ. The new estimator Σ̃w is shown to
possess the same structural properties as Σ.
Future work may follow from this proposal. The criterion α(i) is intended as a guide
and not as a formal test statistic or model selection criterion. Though such ad hoc selection
tools are common in hierarchical clustering, perhaps a more formal statistic could be of
value. For the shrinkage estimation of Σ, investigating adaptive weights w that would
diminish with n and with each iteration of Algorithm 1 may also lead to slightly improved
inference. Working with some, but not all, off-diagonal entries of Σ̂ may also improve the
performance of the procedure. Finally the already mild Partial Exchangeability Assumption
could perhaps be relaxed to allow for more general hierarchical dependence structures.

Supplementary Material
The following material is available online.
Application code: R code to replicate the application on stocks. (.R)
Appendices: Supplementary theoretical and simulation results. (.pdf) This file contains:
Estimation of Σ (Appendix A), Outstanding proofs (Appendix B), Additional results
of the simulation study (Appendix C), and Detailed clusters of the application to stock
returns (Appendix D).

36

Web-based supporting materials for
Detection of Block-Exchangeable Structure
in Large-Scale Correlation Matrices
Samuel Perreault

Thierry Duchesne

Département de mathématiques et de statistique, Université Laval

Johanna G. Nešlehová
Department of Mathematics and Statistics, McGill University
November 23, 2017

1

A

ESTIMATING Σ

Genest et al. (2011) traced back explicit formulas for the diagonal elements of Σ to Lindeberg (1927, 1929). Extending the results of Ehrenberg (1952), they then provide a formula
for the off-diagonal elements of Σ. Using this formula, given in (1) below, we define a
plug-in estimator Σ̂ of Σ. The estimator and its computation are presented in Subsection
A.1. In most cases, it is not advisable to use Σ̂ directly due to a high amount of noise in
the estimation. More so because it needs to be inverted, which is known to amplify the
estimation error when the original matrix is ill-conditioned. Fortunately, if G satisfies the
PEA, Σ has a block structure as well. The latter is described and explained in Subsection
A.2. We then use this block structure to improve the estimation of Σ by averaging entries
of Σ̂ block-wise. The resulting estimate may still contain too much noise to be useful.
Inspired by the work of Ledoit and Wolf (2004), we thus apply, in addition to the averaging
just mentioned, a simple Stein-type shrinkage procedure which depends on a parameter
w, often called the shrinkage intensity. The two shrinkage procedures are presented in
Subsection A.3.

A.1

Plug-in estimator of Σ

Let X be a random vector with continuous univariate marginals F1 , . . . , Fd and a unique
copula C, as in Section 2. Let U = (F1 (X1 ), . . . , Fd (Xd )) and recall that U has distribution
function C. For any subset {i1 , . . . , ik } ⊆ {1, . . . , d} of indices, let Ci1 ···ik denote the unique
copula of the marginal (Xi1 , . . . , Xik ) of X. As shown in Genest et al. (2011), for any
i1 6= j1 ∈ {1, . . . , d} and i2 6= j2 ∈ {1, . . . , d},

Cov(T̂i1 j1 , T̂i2 j2 ) =

4
n(n − 1)

2



n(n − 1)(n − 2)(θ1 + θ2 + θ3 + θ4 )

2

(1)

+ n(n − 1)(ϑ1 + ϑ2 ) −

2(2n − 3)
(Ti1 j1 + 1)(Ti2 j2 + 1),
n(n − 1)

where
θ1 = E(Ci1 j1 (Ui1 , Uj1 )Ci2 j2 (Ui2 , Uj2 )),

θ2 = E(C̄i1 j1 (Ui1 , Uj1 )Ci2 j2 (Ui2 , Uj2 )),

θ3 = E(Ci1 j1 (Ui1 , Uj1 )C̄i2 j2 (Ui2 , Uj2 )),

θ4 = E(C̄i1 j1 (Ui1 , Uj1 )C̄i2 j2 (Ui2 , Uj2 )),

ϑ1 = E(Ci1 j1 i2 j2 (Ui1 , Uj1 , Ui2 , Uj2 )),

(2)

ϑ2 = E(C̃i1 j1 i2 j2 (Ui1 , Uj1 , Ui2 , Uj2 )),

and C̄ denotes the survival function corresponding to C, while
C̃i1 j1 i2 j2 = Ci1 j1 − Ci1 j1 j2 − Ci1 j1 i2 + Ci1 j1 i2 j2 .
For arbitrary r, s ∈ {1, . . . , p}, Σrs = Cov(T̂ir jr , T̂is js ). From (1) and (2) and the fact that
for all i 6= j ∈ {1, . . . , d}, E(Cij (Ui , Uj )) = E(C̄ij (Ui , Uj )), we have, as n → ∞, nΣ → Σ∞ ,
where for any r, s ∈ {1, . . . , p}, the (r, s)-th entry of Σ∞ is given by
(Σ∞ )rs = 16(θ1 + θ2 + θ3 + θ4 ) − 4(Tir jr + 1)(Tis js + 1)
(3)

= 16Cov Cir ,jr (Uir , Ujr ) + C̄ir ,jr (Uir , Ujr ), Cis ,js (Uis , Ujs ) + C̄is ,js (Uis , Ujs ) .
For any i1 6= j1 ∈ {1, . . . , d} and i2 6= j2 ∈ {1, . . . , d}, a plug-in estimator of Cov(T̂i1 j1 , T̂i2 j2 )
can be defined by first replacing Ti1 j1 and Ti2 j2 by T̂i1 j1 and T̂i2 j2 , respectively. Furthermore, the quantities in (2) can be estimated as follows. For k = 1, 2, let I(k) be an n × n
matrix with entries
I(k)
rs := 1(Xrik < Xsik , Xrjk < Xsjk ).

(4)

Similarly to the plug-in estimators considered in Ben Ghorbal et al. (2009), an unbiased
estimator of θ1 is then given by
θ̂1 =

X
1
1(Xri1 < Xsi1 , Xrj1 < Xsj1 )1(Xti2 < Xsi2 , Xtj2 < Xsj2 )
n(n − 1)(n − 2) r6=s6=t
3

=

X
1
(2)
I(1)
rs Its .
n(n − 1)(n − 2) r6=s6=t

Similar formulas can be derived for the other parameters, viz.
θ̂2 =

X
1
(2)
I(1)
rs Itr ,
n(n − 1)(n − 2) r6=s6=t
θ̂4 =

X
1
(2)
I(1)
rs Ist ,
n(n − 1)(n − 2) r6=s6=t
X
(2)
I(1)
rs Irt

θ̂3 =

1
n(n − 1)(n − 2) r6=s6=t

and
ϑˆ1 :=

X
1
(2)
I(1)
rs Irs ,
n(n − 1) r6=s

ϑˆ2 :=

X
1
(2)
I(1)
rs Isr .
n(n − 1) r6=s

Given that θ̂1 , . . . , θ̂4 and ϑ̂1 , ϑ̂2 are U -statistics with square integrable kernels, they are
consistent and asymptotically Normal. These properties carry over to the resulting plug-in
estimator Σ̂ of Σ; in particular,
nΣ̂ → Σ∞

(5)

in probability, as n → ∞. Similarly, one can define a consistent plug-in estimator of Σ∞
by replacing θ1 , . . . , θ4 and Ti1 j1 and Ti2 j2 by their estimators in (3).
Finally, note that Σ̂ can be calculated efficiently using matrix products. To this end,
consider again arbitrary i1 6= j1 ∈ {1, . . . , d} and i2 6= j2 ∈ {1, . . . , d}, and for k = 1, 2,
define I(1) through (4) and set J(1) = (I(1) )> . Furthermore, let 1 be the n-dimensional
vector of ones and ◦ denote the Hadamard product. Then because the diagonal entries of
I(k) , k = 1, 2 are zero,
n(n − 1)(n − 2)

4
X
l=1

=

X h
r6=t

θ̂l =

X 

(2)

(2)

(2)

(2)

(1)
(1)
(1)
I(1)
rs Jst + Jrs Jst + Irs Ist + Jrs Ist

r6=s6=t


i
I(1) + J(1) I(2) + J(2)
rt

4






n
 
o
= 1> I(1) + J(1) I(2) + J(2) 1 − 1> I(1) + J(1) ◦ I(2) + J(2) 1.
Similarly,
n(n − 1)(ϑˆ1 + ϑˆ2 ) =

X

J(1)
rs

h

(2)

I

+J

(2)

i

n

o
(1)
(2)
(2)
=1 J ◦ I +J
1,
>

sr

r6=s

so that
n(n − 1)(n − 2)

4
X

!
θ̂l

+ n(n − 1)(ϑˆ1 + ϑˆ2 )

l=1




n

o
= 1> I(1) + J(1) I(2) + J(2) 1 − 1> J(1) ◦ I(2) + J(2) 1.

A.2

Structure of Σ and Σ−1 implied by G

Suppose that the Partial Exchangeability Assumption (PEA) holds for some partition G.
In this subsection, we describe the block structure of Σ and Σ−1 induced by the PEA.
In Subsection A.3 we then exploit this structure to derive the improved estimator Σ̃ of
Σ. To this end, let us first focus on a single entry Σrs for some arbitrary fixed r, s ∈
{1, . . . , p}. Recall that G induces the partition BG = {B1 , . . . , BL } of {1, . . . , p}, as given
in (8). Equation (1) suggests that the value of Σrs depends on the blocks in BG to which
r and s belong. To identify these blocks, let
Φ1 = {(`1 , `2 ) : 1 ≤ `1 ≤ `2 ≤ L}
be the set of all ordered pairs of block indices and define the function
φ : {1, . . . , p}2 → Φ1
(r, s) 7→ (`1 ∧ `2 , `1 ∨ `2 ) such that (r, s) ∈ B`1 × B`2 ,

5

(6)

where for any a, b ∈ R, a ∧ b = min(a, b) and a ∨ b = max(a, b). Now recall from Section
2 that (ir , jr ) is a pair of indices such that τ r = Tir jr and similarly for (is , js ). The value
of Σrs does not depend only on φ(r, s), but also on the overlap between (ir , jr ) and (is , js ).
To account for the latter, let
Φ2 = {(k1 , k2 ) : 0 ≤ k1 ≤ k2 ≤ K}
and define the function ϕ : {1, . . . , p}2 → Φ2 given by



(0, 0)
if {ir , jr } ∩ {is , js } = ∅,



ϕ(r, s) = (0, k)
if {ir , jr } ∩ {is , js } = {i}, i ∈ Gk ,




(k ∧ k , k ∨ k ) if {i , j } ∩ {i , j } = {i, j}, (i, j) ∈ (G × G ) .
1
2 1
2
r r
s s
k1
k2

(7)

(8)

Using this notation, we introduce, for any ` = (`1 , `2 ) ∈ Φ1 and k = (k1 , k2 ) ∈ Φ2 ,
C`k := {(r, s) ∈ {1, . . . , p}2 : r ≤ s, φ(r, s) = `, ϕ(r, s) = k},

(9)

and set CG := {C`k : (`, k) ∈ Φ1 × Φ2 }. Similarly to TG , we now define the set SG of
matrices with a block structure given by CG , i.e.,
SG := {S ∈ Rp×p : S symmetric and ∀(`, k) ∈ Φ1 × Φ2
(r1 , s1 ), (r2 , s2 ) ∈ C`k ⇒ Sr1 s1 = Sr2 s2 }. (10)
Finally, for each r ∈ {1, . . . , p}, ` ∈ {1, . . . , L} and k ∈ Φ2 , we shall also need the set
(r)

C`k := {s ∈ B` : ϕ(r, s) = k}.

(11)

The next proposition confirms that Σ and Σ∞ have the block structure induced by CG .
Proposition 5. Suppose that G is such that the PEA holds. Then Σ ∈ SG and Σ∞ ∈ SG .
6

Proof. Fix arbitrary (`, k) ∈ Φ1 × Φ2 and (r1 , s1 ), (r2 , s2 ) ∈ C`k . To ease the notation,
write (i1 , j1 ), (i2 , j2 ), (i3 , j3 ) and (i4 , j4 ) instead of (ir1 , jr1 ), (is1 , js1 ), (ir2 , jr2 ) and (is2 , js2 ),
respectively. To prove the claim, it suffices to show all expectations in (1) are identical
when (i1 , j1 ) is changed to (i3 , i3 ) and (i2 , j2 ) to (i4 , j4 ), respectively. Focusing on θ1 , we
want to show
E[Ci1 j1 (Ui1 , Uj1 )Ci2 j2 (Ui2 , Uj2 )] = E[Ci3 j3 (Ui3 , Uj3 ) Ci4 j4 (Ui4 , Uj4 )].

(12)

Let I be the set of unique indices from (i1 , j1 , i2 , j2 ) and Im = (im , jm ), m = 1, 2. The LHS
of (12) can then be rewritten as
Z
E(Ci1 j1 (Ui1 , Uj1 )Ci2 j2 (Ui2 , Uj2 )) =

Ci1 j1 (ui1 , uj1 )Ci2 j2 (ui2 , uj2 ) dCi1 j1 i2 j2
Z

=

CI1 (uI1 )CI2 (uI2 ) dCI .

(13)

Now define J to be the set of distinct indices from (i3 , j3 , i4 , j4 ) and Jm = (im+2 , jm+2 ),
m = 1, 2. Combining the facts that φ(r1 , s1 ) = φ(r2 , s2 ) and ϕ(r1 , s1 ) = ϕ(r2 , s2 ), we
deduce that for any k ∈ {1, . . . , K}, I and J have the same number of entries coming
from Gk , with no repetition. The same can be deduced for Im and Jm , m = 1, 2. We can
therefore use the PEA to replace CI by CJ and CIm by CJm for m = 1, 2. Consequently,
Z
Z
CI1 (uI1 )CI2 (uI2 ) dCI = CJ1 (uJ1 )CJ2 (uJ2 ) dCJ ,
thus showing that (13) is indeed equal to the RHS of (12). Equalities for the other quantities
θ2 , . . . , θ4 and ϑ1 , ϑ2 can be shown using the same technique.
Example 7. For G as given in (4) in Example 1, there are L = 6 sets in BG , viz.
BG = {B11 , B12 , B13 , B22 , B23 , B33 }≡ {B1 , . . . , B6 },

7

(14)

G1
B11

G1
B11

G2

G3

B12

B13

B22
G2

B23

B12
B22

B33
G3

B13

B23
B33

T

Sub-matrix of Σ

Figure 10: The matrix T (left) and a sub-matrix of Σ (right) from Example 7. The cells
are tinted so that, in each matrix, all entries sharing the same value are of the same colour
and colour intensity.
where for k1 , k2 ∈ {1, . . . , 3}, Bk1 k2 is as in (6). The blocks in BG are displayed in left panel
of Figure 10. In this picture, for each k1 , k2 ∈ {1, 2, 3}, the cells (ir , jr ) and (jr , ir ) for
r ∈ Bk1 k2 are tinted with the same colour, emphasizing that the entries are equal.
Given that p = d(d + 1)/2 = 9(10)/2 = 45, the 45 × 45 matrix Σ is more cumbersome to visualize. To see its structure more clearly, we vectorize T not as in (5),
but rather block by block. For instance the first 18 entries of τ are the 6 entries in T
corresponding to B11 ≡ B1 = {1, . . . , 6} followed by the 12 entries in T corresponding to
B12 ≡ B2 = {7, . . . , 18}, i.e. (τ 1 , . . . , τ 18 ) = (T1,2 , T1,3 , T1,4 , T2,3 , T2,4 , T3,4 , T1,5 , . . . , T4,7 ).
The 18×18 dimensional sub-matrix of Σ displaying the pair-wise covariances of τ̂ 1 , . . . , τ̂ 18
is showed in the right panel of Figure 10. Distinct values are depicted using different colour
and colour intensity. The colours represent distinct values of φ: for all r, s ∈ {1, . . . , 18},
the cell (r, s) is tinted red, magenta and violet if φ(r, s) equals (1, 1), (1, 2) and (2, 2), respectively. In other words, the coloured blocks are induced by B11 and B12 . Next, notice
8

that in each coloured block, the values of Σ can differ, and this is depicted through different colour intensity. This is because for any r, s ∈ {1, . . . , p}, the value of Σrs also
depends on ϕ(r, s). For example, the red block in the top left corner contains three distinct values, for if r, s ∈ {1, . . . , 6}, {ir , jr } ∩ {is , js } is either empty (red), contains one
element from G1 (light red), or contains two elements from G1 (dark red). To illustrate,
consider r = 1, 2, 6. Then (i1 , j1 ) = (1, 2), (i2 , j2 ) = (1, 3), and (i6 , j6 ) = (3, 4). Consequently, φ(1, 1) = φ(1, 2) = φ(1, 6) = (1, 1), and indeed the entries Σ11 , Σ12 and Σ16 are
tinted red. However, ϕ(1, 1) = (1, 1), ϕ(1, 2) = (0, 1) and ϕ(1, 6) = (0, 0), which is why
Σ11 , Σ12 and Σ16 are dark red, red and light red, respectively. One can indeed verify that
Σ11 6= Σ12 6= Σ16 . This shows that the block structure Σ is described by both φ and ϕ; the
sets in CG correspond to the cells above the main diagonal of Σ with the same colour and
intensity.
(r)

Finally, the right panel in Figure 10 can be used to visualize the sets C`k defined in
(r)

(11). For a given r ∈ {1, . . . , p}, the union of the sets C`k , ` ∈ {1, . . . , L}, k ∈ Φ2 may be
identified with the r-th row (or equivalently the r-th column) of Σ: the index ` determines
the colour and k the intensity in that row (column). To illustrate in the context of this
example, pick r = 1 and ` = 1, say. Then
(1)

C1(1,1) = {1},

(1)

C1(0,1) = {2, 3, 4, 5},

(1)

(1)

(1)

C1(0,0) = {6},
(1)

(1)

while C1k = ∅ for any other k ∈ Φ2 . Note that {C1(1,1) , C1(0,1) , C1(0,0) } is a partition of B11 .
The loss function ` depends on Σ through its inverse Σ−1 . The following proposition
establishes that the structure of Σ−1 is the same as that of Σ.
Proposition 6. Suppose that G is a partition for which the PEA holds. An invertible matrix
S is an element of SG if and only if S−1 ∈ SG . That is, SG is closed under inversion.
9

Proof. From the Cayley-Hamilton Theorem (Harville, 2008, p.583) S satisfies its characteristic equation
p

S +

p−1
X

cs Ss + (−1)p |S|Ip = 0,

s=1

for some known coefficients cs , s = 1, . . . , p − 1. As a consequence, the inverse of a p × p
matrix S can be represented by a linear function of its p − 1 first powers:
p

−1

S

1 X
cs Ss−1 .
=
|S| s=1

Therefore, it suffices to show that if S, Q ∈ SG and SQ is symmetric, SQ ∈ SG . To
this end, fix an arbitrary ` = (`1 , `2 ) ∈ Φ1 , k = (k1 , k2 ) ∈ Φ2 , and arbitrary pairs
(r1 , s1 ), (r2 , s2 ) ∈ C`k . To show that [SQ]r1 s1 = [SQ]r2 s2 , first note that because SQ is
symmetric by assumption, it can be assumed, without loss of generality, that r1 , r2 ∈ B`1
and s1 , s2 ∈ B`2 . For any `∗ ∈ Φ1 and any k∗ ∈ Φ2 , let S`
`∗ k ∗

values such that Srs = S

and Qrs = Q

`∗ k ∗

∗ ∗

k

and Q`

∗ ∗

k

denote the unique

whenever (r, s) ∈ C`∗ k∗ . Because BG given

by (8) is a partition of {1, . . . , p}, we can write
[SQ]r1 s1 =

p
X

Sr1 t Qts1 =

t=1

L X
X

Sr1 t Qts1 .

(15)

`=1 t∈B`

For any fixed ` ∈ {1, . . . , d} and t ∈ B` , φ(r1 , t) = (` ∧ `1 , ` ∨ `1 ) ≡ `1` and φ(s1 , t) =
(` ∧ `2 , ` ∨ `2 ) ≡ `2` . However, Sr1 t and Qts1 also depend on ϕ(r1 , t) and ϕ(s1 , t), and this
requires further partitioning of B` by means of the sets defined in (11). Specifically,
B` =

[

(r )

(s )

(C`k11 ∩ C`k12 ).

k1 ,k2 ∈Φ2
(r )

(s )

Clearly, the sets (C`k11 ∩ C`k12 ) are disjoint for distinct k1 , k2 ∈ Φ2 , and for any given
(r )

(s )

k1 , k2 ∈ Φ2 and t ∈ C`k11 ∩ C`k12 , ϕ(r1 , t) = k1 and ϕ(t, s1 ) = k2 so that Sr1 t = S`1` k1 and
10

Qts1 = Q`2` k2 . Consequently, the last expression in (15) can be rewritten as
L
X
`=1

X

X

`1` k1

S

Q

`2` k2

=

L
X

k1 ,k2 ∈Φ2 t∈C (r1 ) ∩C (s1 )
`k1
`k2

X

(r )

(s )

C`k11 ∩ C`k12 S`1` k1 Q`2` k2 .

(16)

`=1 k1 ,k2 ∈Φ2

(r )

(s )

Now for any ` ∈ {1, . . . , L} and any k1 , k2 ∈ Φ2 , Lemma 3 gives that |C`k11 ∩ C`k12 | =
(r )

(s )

(r )

(s )

|C`k21 ∩ C`k22 |. Furthermore, for any t ∈ C`k21 ∩ C`k12 , φ(r2 , t) = `1` , φ(t, s2 ) = `2` and
ϕ(r2 , t) = k1 , ϕ(t, s2 ) = k2 , so that Sr2 t = S`1` k1 , and Qts2 = Q`2` k2 . Consequently, the
RHS in (16) equals
L
X

X

(r )

(s )

C`k21 ∩ C`k22 S`1` k1 Q`2` k2 =

`=1 k1 ,k2 ∈Φ2

L
X

X

X

Sr2 t Qts2 = [SQ]r2 s2 ,

`=1 k1 ,k2 ∈Φ2 t∈C (r2 ) ∩C (s2 )
`k
`k
1

2

as claimed.
In view of Proposition 5, the following result follows directly from Proposition 6.
Corollary 1. Suppose that G is such that the PEA holds. If Σ is invertible, then Σ−1 ∈ SG .
Similarly, if Σ∞ is invertible, then Σ−1
∞ ∈ SG .

A.3

An improved estimator of Σ

Throughout this section, assume that G is a partition of {1, . . . , d} such that the PEA
holds. The empirical estimator Σ̂ defined in Section A.1 does not exploit the structural
information provided by CG . It is natural to think that, as was the case for τ , we can
improve the estimation of Σ by averaging its entries with respect to the sets in CG . But we
can do even better by exploiting the following decomposition of Σ. To this end, write
Σ=Θ−

2(2n − 3)
(τ + 1)(τ + 1)>
n(n − 1)
11

(17)

Σ

Θ

(τ + 1)(τ + 1)>

Figure 11: Submatrices of Σ, Θ and (τ + 1)(τ + 1)> from Example 8. The same vectorization of T as in Example 7 is used.
where 1 is the p-dimensional vector of ones and Θ is a p × p matrix gathering the terms
involving θ1 , . . . , θ4 and ϑ1 , ϑ2 in (1). It easily follows from the proof of Proposition 5 that
Θ ∈ SG as well as (τ + 1)(τ + 1)> ∈ SG . However, the structure of (τ + 1)(τ + 1)> is even
simpler, because the overlaps between pairs of indices described by the function ϕ need not
be taken into account. Specifically, (τ + 1)(τ + 1)> ∈ TBG ⊂ SG , where
TBG = {R ∈ Rp×p : ∀`1 , `2 ∈ {1, . . . , L}, r1 , r2 ∈ B`1 and s1 , s2 ∈ B`2 ⇒ Rr1 s1 = Rr2 s2 }.
That is, (τ + 1)(τ + 1)> possesses a block structure similar to T, but defined in accordance
with the clustering BG instead of G. In particular, (τ + 1)(τ + 1)> possesses L diagonal
blocks, as opposed to K diagonal blocks for T (counting the diagonal blocks that correspond
to clusters in G of size 1 as well). The decomposition (17) of Σ is illustrated next.
Example 8. The decomposition of Σ from Example 1 according to (17) is depicted in
Figure 11. The matrix Σ clearly inherits its structure from Θ ∈ SG and the structure of
(τ + 1)(τ + 1)> is considerably simpler.
Let Θ̂ be the plug-in empirical estimator of Θ, defined by replacing, for each (r, s) ∈
12

{1, . . . , p}2 , the parameters θ1 , . . . , θ4 and ϑ1 , ϑ2 by their empirical estimates given in Section A.1. Because Θ ∈ SG , we now define the improved estimator Θ̃, which is in SG by
construction. First, the upper triangular (including the diagonal) entries of Θ̃ are simply
the entries of Θ̂ averaged out over each block in CG . Second, because Θ̃ is symmetric its
lower triangular entries are obtained by symmetry. Furthermore, let τ̃ = τ̃ (τ̂ |G) be as in
(15) and estimate (τ + 1)(τ + 1)> by (τ̃ + 1)(τ̃ + 1)> ∈ TBG , so that even more averaging
is employed. The resulting estimator Σ is then
Σ̃ = Θ̃ −

2(2n − 3)
(τ̃ + 1)(τ̃ + 1)> .
n(n − 1)

Clearly, Σ̃ ∈ SG . Using (3), write Σ∞ = Θ∞ − 4(τ + 1)(τ + 1)> , where Θ∞ has entries
16(θ1 + · · · + θ4 ). Note that from the proof of Proposition 5, Θ∞ ∈ SG . Hence, as n → ∞,
nΘ̃ → Θ∞ element-wise in probability. Furthermore, given that τ̃ is a consistent estimator
of τ as per Theorem 2, (τ̃ + 1)(τ̃ + 1)> → (τ + 1)(τ + 1)> . Put together,
nΣ̃ → Σ∞

(18)

element-wise in probability.
When n is small compared to d, and |G| is large, not enough averaging is operated and
Σ̃ may be too noisy to be useful. We therefore apply (Steinian) shrinkage and consider
Σ̃w := Σ̃(Σ̂|τ̂ , G, w) := (1 − w)Σ̃ + wΣ̃diag ,

(19)

where Σ̃diag is the diagonal matrix whose non zero elements are the elements on the diagonal
of Σ̃ and w ∈ [0, 1] is the shrinkage intensity parameter. This type of “linear shrinkage”
follows the proposal in Devlin et al. (1975). It is easy to show that Σ̃w ∈ SG . Because Σ̃
is consistent, so is Σ̃w , as long as w → 0 as n → ∞.
Finally, note that the estimators Σ̂, Σ̃, and Σ̃w may not be positive definite, in particular
when n is small. For the methodology presented in this paper, Σ̃w needs to be invertible,
13

and this is often easily achieved by a larger value of the shrinkage parameter w when n is
small; in the data illustration and simulation study conducted in this paper no problems
with invertibility were encountered. If the estimator of Σ further needs to be positive semidefinite, one can project any of the estimators Σ̂, Σ̃, and Σ̃w to the cone of positive semidefinite matrices; this projection, say Σ̄, can be computed using the alternative direction
of multipliers algorithm as described, e.g., in Appendix A of Datta and Zou (2017). Since
the projection onto the cone of positive semidefinite matrices is a continuous mapping, Σ̄
will be consistent.

B

OUTSTANDING PROOFS

This section is divided as follows. Subsection B.1 contains six auxiliary lemmas that pertain
to the structure of SG and that are invoked in the proofs of Theorem 1, Theorem 2 and
Proposition 4. The proof of Proposition 4 is detailed in Subsection B.2 while Subsection
B.3 contains additional results on the structural properties of the inverse of correlation
matrices used in Section 5 of the paper.

B.1

Description of SG
(r)

We first present three auxiliary lemmas that pertain to the cardinality of the sets C`k
defined in (11).

Lemma 1. Let G be an arbitrary partition of {1, . . . , d} and BG as in (8). Assume that
(r)

r, s ∈ B` for some ` ∈ {1, . . . , L}. Then for all λ ∈ {1, . . . , L}, and all κ ∈ Φ2 , |Cλκ | =
(s)

|Cλκ |.
Remark 5. Before proceeding with the proof of Lemma 1, let us illustrate the claim on the
14

(r)

right panel of Figure 10. Because Σ ∈ SG , the sets C`k for a fixed r can be identified with
the r-th row of Σ; different colours and intensities correspond to different values of ` and
k, respectively; see also Example 7. Lemma 1 implies that, for example, the number of cells
with the same colour and intensity in the first six rows (corresponding to B1 ) is the same.
Proof. Fix arbitrary G, ` ∈ {1, . . . , L} and r, s ∈ B` . We will prove the assertion by
showing that there exists a bijective function h : {1, . . . , p} → {1, . . . , p} such that for all
λ ∈ {1, . . . , L} and κ ∈ Φ2 ,
(r)

t ∈ Cλκ
(r)

(s)

⇔

h(t) ∈ Cλκ ,

(20)

(s)

for then obviously |Cλκ | = |Cλκ |. To this end, first identify k1 , k2 , ∈ {1, . . . , K} such that
(ir , jr ) ∈ Gk1 × Gk2 . Note that then B` = B(k1 ∧k2 )(k1 ∨k2 ) . Because s ∈ B` by assumption,
either (is , js ) ∈ Gk1 × Gk2 or (is , js ) ∈ Gk2 × Gk1 .
First assume that (is , js ) ∈ Gk1 × Gk2 . Let π be any permutation such that
∀ i ∈ {1, . . . , d} ∀ k ∈ {1, . . . , K},

i ∈ Gk ⇔ π(i) ∈ Gk ,

(21)

and further such that
π(ir ) = is ,

π(jr ) = js .

(22)

Because ir 6= jr and is 6= js , and ir , is ∈ Gk1 , jr , js ∈ Gk2 , such a permutation always exists,
although it is generally not unique. Now define h by
h : {1, . . . , p} → {1, . . . , p}

(23)

t 7→ t∗ such that it∗ = π(it ) ∧ π(jt ) and jt∗ = π(it ) ∨ π(jt ).
First, observe that h is well defined because for any t ∈ {1, . . . , p}, it < jt , π(it ) 6= π(jt )
and hence it∗ < jt∗ , so that t∗ indeed exists. Furthermore, h is a bijection. This is because,
15

for any t∗ ∈ {1, . . . , p} and t such that it = π −1 (it∗ ) ∧ π −1 (jt∗ ) and jt = π −1 (it∗ ) ∨ π −1 (jt∗ ),
h(t) = t∗ . Furthermore, (21) implies that for any λ ∈ {1, . . . , L}, t ∈ Bλ if and only if
h(t) ∈ Bλ . To prove that h satisfies (20), it thus remains to show that for any κ ∈ Φ2 ,
ϕ(r, t) = κ if and only if ϕ(s, h(t)) = κ. This follows from the following facts, each of
which is an immediate consequence of (22):
(i) For any i ∈ {1, . . . , d}, {i} ∩ {ir , jr } = ∅ if and only if {π(i)} ∩ {is , js } = ∅.
(ii) For any i ∈ {1, . . . , d}, {i} ∩ {ir , jr } = {ir } if and only if {π(i)} ∩ {is , js } = {is }; ir
and is are elements of the same cluster Gk1 .
(iii) For any i ∈ {1, . . . , d}, {i} ∩ {ir , jr } = {jr } if and only if {π(i)} ∩ {is , js } = {js }; jr
and js are elements of the same cluster Gk2 .
This concludes the proof in the case when (is , js ) ∈ Gk1 ×Gk2 . When (is , js ) ∈ Gk2 ×Gk1 , one
can proceed analogously by constructing h from an arbitrary fixed permutation π satisfying
(21) and such that π(ir ) = js and π(jr ) = is .
Lemma 2. Let G be an arbitrary partition of {1, . . . , d} and BG as in (8). Then for any
`1 , `2 ∈ {1, . . . , L}, r ∈ B`1 , s ∈ B`2 , and k ∈ Φ2 ,
(s)

(r)

|C |
|C`2 k |
= `1 k .
|B`2 |
|B`1 |
Proof. To prove the claim, fix arbitrary `1 , `2 ∈ {1, . . . , L}, r ∈ B`1 , s ∈ B`2 , and k ∈ Φ2 .
The case `1 = `2 trivially follows from Lemma 1. For `1 6= `2 , set ` = (`1 ∧ `2 , `1 ∨ `2 ).
Using Lemma 1 again, we then have that
|C`k | =

X
t∈B`1

(t)

|C`2 k | =

X

(r)

(r)

|C`2 k | = |C`2 k | |B`1 |.

t∈B`1
(s)

Similarly, summing over elements in B`2 , |C`k | = |C`1 k | |B`2 |, which proves the claim.
16

Lemma 3. Let G be an arbitrary partition of {1, . . . , d} and BG as in (8). Assume
that r1 , r2 ∈ B`1 and s1 , s2 ∈ B`2 for some ` = (`1 , `2 ) ∈ Φ1 . Further assume that
(r1 , s1 ), (r2 , s2 ) ∈ C`k for some k ∈ Φ2 . Then for all λ ∈ {1, . . . , L}, and all κ1 , κ2 ∈ Φ2 ,
(r )

(s )

(r )

(s )

one has |Cλκ11 ∩ Cλκ12 | = |Cλκ21 ∩ Cλκ22 |, i.e.,
|{t ∈ Bλ : ϕ(r1 , t) = κ1 , ϕ(s1 , t) = κ2 }| = |{t ∈ Bλ : ϕ(r2 , t) = κ1 , ϕ(s2 , t) = κ2 }|.
Proof. To show this claim, we can proceed similarly as in the proof of Lemma 1. We will
again define a function h through (23) from a certain permutation π of (1, . . . , d) satisfying
(21). As argued in the proof of Lemma 1, h is then a well-defined bijection such that for
all λ ∈ {1, . . . , L}, t ∈ Bλ if and only if h(t) ∈ Bλ . If h is further such that for each
t ∈ {1, . . . , p} and κ1 , κ2 ∈ Φ2 ,
ϕ(r1 , t) = κ1

⇔

ϕ(r2 , h(t)) = κ1

and ϕ(s1 , t) = κ2
(r )

⇔

ϕ(s2 , h(t)) = κ2 , (24)

(s )

(r )

(s )

then it holds that for all t ∈ {1, . . . , p}, t ∈ Cλκ11 ∩ Cλκ12 if and only if h(t) ∈ Cλκ21 ∩ Cλκ22 ,
(r )

(s )

(r )

(s )

and consequently that |Cλκ11 ∩ Cλκ12 | = |Cλκ21 ∩ Cλκ22 |, as claimed.
In contrast to the proof of Lemma 1, the properties of the permutation π require a
more cumbersome case distinction. To this end, fix ` = (`1 , `2 ) ∈ Φ1 , k ∈ Φ2 , and
(r1 , s1 ), (r2 , s2 ) ∈ C`k such that r1 , r2 ∈ B`1 and s1 , s2 ∈ B`2 . Now let k11 , k12 , k21 , k22 ∈
{1, . . . , K} be such B`1 = Bk11 k12 and B`2 = Bk21 k22 . Without loss of generality, assume that
(ir1 , jr1 ), (ir2 , jr2 ) ∈ Gk11 × Gk12

and (is1 , js1 ), (is2 , js2 ) ∈ Gk21 × Gk22 .

(25)

Case I. k = ϕ(r1 , s1 ) = ϕ(r2 , s2 ) = (0, 0). Here,
|{ir1 , jr1 , is1 , js1 }| = |{ir2 , jr2 , is2 , js2 }| = 4,

(26)

and π is an arbitrary fixed permutation with the property (21) and such that
π(ir1 ) = ir2 ,

π(jr1 ) = jr2 ,

π(is1 ) = is2 ,
17

π(js1 ) = js2 .

(27)

Such a permutation exists because of (25) and (26), although it is generally not unique.
Because of (27), it holds that, for any i ∈ {1, . . . , d},
(i) {i} ∩ {ir1 , jr1 } = ∅ if and only if {π(i)} ∩ {ir2 , jr2 } = ∅, and {i} ∩ {is1 , js1 } = ∅ if and
only if {π(i)} ∩ {is2 , js2 } = ∅;
(ii) {i} ∩ {ir1 , jr1 } = {ir1 } if and only if {π(i)} ∩ {ir2 , jr2 } = {ir2 }; ir1 and ir2 are elements
of the same cluster Gk11 .
(iii) {i}∩{ir1 , jr1 } = {jr1 } if and only if {π(i)}∩{ir2 , jr2 } = {jr2 }; jr1 and jr2 are elements
of the same cluster Gk12 .
(iv) {i} ∩ {is1 , js1 } = {is1 } if and only if {π(i)} ∩ {is2 , js2 } = {is2 }; is1 and is2 are elements
of the same cluster Gk21 .
(v) {i}∩{is1 , js1 } = {js1 } if and only if {π(i)}∩{is2 , js2 } = {js2 }; js1 and js2 are elements
of the same cluster Gk22 .
Hence h fulfills (24) and the proof is complete in this case.
Case II. k = ϕ(r1 , s1 ) = ϕ(r2 , s2 ) = (0, k) for some k ∈ {1, . . . , K}. In this case,
|{ir1 , jr1 , is1 , js1 }| = |{ir2 , jr2 , is2 , js2 }| = 3.

(28)

Observe that the three distinct elements of {ir1 , jr1 , is1 , js1 } are {a1 , a2 , a3 }, say, such that
a1 ∈ {ir1 , jr1 } ∩ {is1 , js1 }, a2 ∈ {ir1 , jr1 } \ {a1 }, and a3 ∈ {is1 , js1 } \ {a1 }. Similarly,
{ir2 , jr2 , is2 , js2 } = {b1 , b2 , b3 } such that b1 ∈ {ir2 , jr2 } ∩ {is2 , js2 }, b2 ∈ {ir2 , jr2 } \ {b1 }, and
b3 ∈ {is2 , js2 } \ {b1 }. Note that then necessarily
{ir1 , jr1 } = {a1 , a2 }, {is1 , js1 } = {a1 , a3 }, {ir2 , jr2 } = {b1 , b2 }, {is2 , js2 } = {b1 , b3 }.

18

Furthermore, for m ∈ {1, 2, 3}, am and bm are members of the same cluster. Indeed, given
that k = (0, k), a1 , b1 ∈ Gk . From (25) it further follows that a2 , b2 are in Gk11 and Gk12 ,
respectively, if a2 = ir1 , b2 = ir2 , and a2 = jr1 , b2 = jr2 , respectively. If a2 = ir1 and
b2 = jr2 , then a1 = jr1 and b1 = ir2 and (25) together with the fact that a1 , b1 ∈ Gk imply
that k = k11 = k12 . Hence, a2 , b2 ∈ Gk . Similarly, if a2 = jr1 and b2 = ir2 , we also have that
a2 , b2 ∈ Gk . The verification of the fact that a3 and b3 are in the same cluster is analogous.
Now let π be any permutation with the property (21) and such that
π(a1 ) = b1 ,

π(a2 ) = b2 ,

π(a3 ) = b3 .

(29)

Such a permutation indeed exists, although it is again generally not unique; existence is
guaranteed because the am ’s are all distinct and for m ∈ {1, 2, 3}, am and bm are members
of the same cluster. Because of (29), it holds that, for any i ∈ {1, . . . , d},
(i) {i} ∩ {ir1 , jr1 } = ∅ if and only if {π(i)} ∩ {ir2 , jr2 } = ∅, and {i} ∩ {is1 , js1 } = ∅ if and
only if {π(i)} ∩ {is2 , js2 } = ∅;
(ii) {i} ∩ {ir1 , jr1 } = {a1 } if and only if {π(i)} ∩ {ir2 , jr2 } = {b1 }; a1 and b1 are elements
of the same cluster Gk .
(iii) {i} ∩ {ir1 , jr1 } = {a2 } if and only if {π(i)} ∩ {ir2 , jr2 } = {b2 }; a2 and b2 are elements
of the same cluster Gk11 , Gk12 or Gk , as the case may be.
(iv) {i} ∩ {is1 , js1 } = {a1 } if and only if {π(i)} ∩ {is2 , js2 } = {b1 }; a1 and b1 are elements
of the same cluster Gk .
(v) {i} ∩ {is1 , js1 } = {a3 } if and only if {π(i)} ∩ {is2 , js2 } = {b3 }; a3 and b3 are elements
of the same cluster Gk21 , Gk22 or Gk , as the case may be.

19

Hence h fulfills (24) and the proof is complete in this case.
Case III. k = ϕ(r1 , s1 ) = ϕ(r2 , s2 ) = (k1 , k2 ) where k1 , k2 > 0. In this case,
|{ir1 , jr1 , is1 , js1 }| = |{ir2 , jr2 , is2 , js2 }| = 2.

(30)

Write {ir1 , jr1 , is1 , js1 } = {a1 , a2 }, and {ir2 , jr2 , is2 , js2 } = {b1 , b2 }, with a1 , b1 ∈ Gk1 and
a2 , b2 ∈ Gk2 . Now let π be any permutation with the property (21) and such that
π(a1 ) = b1

π(a2 ) = b2 .

(31)

Because of (31) and the fact that {a1 , a2 } = {ir1 , jr1 } ∩ {is1 , js1 }, and {b1 , b2 } = {ir2 , jr2 } ∩
{is2 , js2 }, it holds that, for any i ∈ {1, . . . , d},
(i) {i} ∩ {ir1 , jr1 } = ∅ if and only if {π(i)} ∩ {ir2 , jr2 } = ∅, and {i} ∩ {is1 , js1 } = ∅ if and
only if {π(i)} ∩ {is2 , js2 } = ∅;
(ii) {i} ∩ {ir1 , jr1 } = {a1 } if and only if {π(i)} ∩ {ir2 , jr2 } = {b1 }; a1 and b1 are elements
of the same cluster Gk1 .
(iii) {i} ∩ {ir1 , jr1 } = {a2 } if and only if {π(i)} ∩ {ir2 , jr2 } = {b2 }; a2 and b2 are elements
of the same cluster Gk2 .
(iv) {i} ∩ {is1 , js1 } = {a1 } if and only if {π(i)} ∩ {is2 , js2 } = {b1 }; a1 and b1 are elements
of the same cluster Gk1 .
(v) {i} ∩ {is1 , js1 } = {a2 } if and only if {π(i)} ∩ {is2 , js2 } = {b2 }; a2 and b2 are elements
of the same cluster Gk2 .
Hence h fulfills (24) and the proof is complete in this case as well.
The following lemma is the cornerstone of the proof of Theorem 1.
20

Lemma 4. Let G be an arbitrary partition of {1, . . . , d}. If S ∈ SG and B is as in (9),
then B> S = B> SBB+ .
Proof. First note that for any ` ∈ {1, . . . , L} and r ∈ {1, . . . , p},
>

[B S]`r =

p
X

Bs` Ssr =

s=1

X

Ssr

(32)

s∈B`

and that [B+ ]`r = 1(r ∈ B` )|B` |−1 . Also, for any s ∈ B` ,
[BB+ ]rs = 1(r ∈ B` )|B` |−1 .

(33)

Now fix an arbitrary ` ∈ {1, . . . , L} and r ∈ {1, . . . , p} and find `∗ so that r ∈ B`∗ . Then
from (32) and (33),
[B> SBB+ ]`r =

p
X

[B> S]`s [BB+ ]sr =

s=1

p
X
X
s=1

t∈B`

=

!
Sts

1(s ∈ B`∗ )

|B`∗ |
X X Sts

s∈B`∗ t∈B`

|B`∗ |

=

X 1 X
Sts . (34)
|B
`∗ |
t∈B
s∈B ∗
`

`

Now set ` = (` ∧ `∗ , ` ∨ `∗ ) and for any k ∈ Φ2 , let S`k denote the unique value such that
(t)

Sts = S`k whenever (t, s) ∈ C`k . Because B`∗ = ∪k∈Φ2 C`∗ k for any t ∈ B` and S ∈ SG by
assumption,
X

Sts =

s∈B`∗

X X

Sts =

k∈Φ2 s∈C (t)
`∗ k

X X
k∈Φ2 s∈C (t)
`∗ k

S`k =

X

(t)

|C`∗ k |S`k .

k∈Φ2

Using this and Lemma 2, the RHS in (34) equals
X X |C (t)
X X |C (r) |
X (r)
`∗ k | `k
`k
S =
S`k =
|C`k |S`k .
|B
|B
`∗ |
`|
t∈B k∈Φ
t∈B k∈Φ
k∈Φ
`

2

`

2

21

2

(r)

Now use the fact that B` = ∪k∈Φ2 C`k and that S ∈ SG to rewrite the RHS as
X X

S`k =

X

Str ,

t∈B`

k∈Φ2 t∈C (r)
`k

which is equal to the RHS in (32), as claimed.
Next, note the following result, which is an immediate consequence of the properties of
the Moore-Penrose pseudo-inverse.
Lemma 5. Let G be an arbitrary partition of {1, . . . , d}, B as in (9) and Γ = BB+ . Then
Γ and (Ip − Γ) are idempotent, i.e., ΓΓ = Γ and (Ip − Γ)(Ip − Γ) = (Ip − Γ).
Lemma 6. Let G be an arbitrary partition of {1, . . . , d}, B as in (9) and Γ = BB+ . Then
for any S ∈ SG , ΓS = ΓSΓ = SΓ.
Proof. For arbitrary ` ∈ Φ1 and k ∈ Φ2 , let S`k denote the unique value such that
Srs = S`k whenever (r, s) ∈ C`k . Now fix arbitrary r, s ∈ {1, . . . , p} and let `1 , `2 be such
that r ∈ B`1 and s ∈ B`2 . Set ` = (`1 ∧`2 , `1 ∨`2 ). From (33) and Lemma 2 we can conclude
that
X X S`k
X |C`(s)k |
X |C`(r)ϕ |
X Sts
2
`k
1
=
=
S =
S`k = [SΓ]rs ,
[ΓS]rs =
|B
|
|B
|
|B
|
|B
|
`1
`1
`1
`2
(s)
t∈B
k∈Φ
k∈Φ
k∈Φ
`1

2

2

t∈C`

1k

2

proving that ΓS = SΓ. Furthermore, ΓS = SΓ and Lemma 5 together imply that ΓΓS =
ΓSΓ. Using the idempotence of Γ again, this simplifies to ΓS = ΓSΓ.
Remark 6. Suppose that G satisfies the PEA. Proposition 5 and Corollary 1 along with
Lemma 5 imply that (Ip − Γ)S = (Ip − Γ)S(Ip − Γ) = S(Ip − Γ) holds for both S = Σ and
S = Σ−1 .

22

B.2

Proof of Proposition 4

Because G satisfies the PEA, Γτ = τ . Therefore,
−1

−1

(τ̂ − τ̃ )> Σ̂ (τ̂ − τ̃ ) = (τ̂ − τ − Γ(τ̂ − τ ))> Σ̂ (τ̂ − τ − Γ(τ̂ − τ ))
−1

= (τ̂ − τ )> (Ip − Γ)Σ̂ (Ip − Γ)(τ̂ − τ ).
Because A 7→ A−1 is a continuous mapping on the space of nonsingular matrices (Stewart,
−1

1969), Σ̂ /n converges element-wise to Σ−1
∞ in probability as n → ∞. The asymptotic
Normality (13) implies along with Slutsky’s lemma that
−1

(τ̂ − τ̃ )> Σ̂ (τ̂ − τ̃ )

V> Σ−1
∞ V,

where V ∼ N (0p , (Ip − Γ)Σ∞ (Ip − Γ)). Now set M := (Ip − Γ)Σ∞ (Ip − Γ) and A := Σ−1
∞.
Following Lemma 5 and Remark 6, we easily get that
MA = (Ip − Γ)Σ∞ (Ip − Γ)Σ−1
∞ = (Ip − Γ)
is idempotent and has trace tr(Ip − Γ) = p − tr(Γ). An application of Theorem 8.6 in
2
Severini (2005) thus yields that V> Σ−1
∞ V is χp−tr(Γ) . The claim now follows from
p
X

L X
X
1
tr(Γ) =
Γrr =
= L.
|B
|
`
r=1
`=1 r∈B
`

B.3

Additional results for Section 5 of the paper

In this section, we consider inverses of matrices in TG . To this end, consider an arbitrary
partition G of {1, . . . , d} and introduce constraints on the diagonal entries of the matrices
in TG through the set
TG† := {R ∈ TG : for any i, j such that ∆ij = 1, Rii = Rjj }.
23

A direct consequence of this definition is that TG† ⊂ TG . Furthermore, T ∈ TG† since Tii = 1
for all i ∈ {1, . . . , d}.
Lemma 7. If R ∈ TG† is invertible, R−1 ∈ TG† .
Proof. Invoking once again the Cayley-Hamilton Theorem as in the proof of Proposition
6, it suffices to show R, Q ∈ TG† implies that RQ ∈ TG† . To this end, fix arbitrary k1 , k2 ∈
{1, . . . , K} and let (i1 , j1 ), (i2 , j2 ) ∈ Gk1 × Gk2 be such that i1 = j1 if and only if i2 = j2 .
Then for any k ∈ {1, . . . , K},
X

Ri1 s Qsj1 =

X

Ri2 s Qsj2 .

s∈Gk

s∈Gk

Therefore,
[RQ]i1 j1 =

d
X

Ri1 s Qsj1 =

s=1

K X
X

Ri1 s Qsj1 =

K X
X

Ri2 s Qsj2 = [RQ]i2 j2 ,

k=1 s∈Gk

k=1 s∈Gk

which proves the claim.

C

ADDITIONAL RESULTS OF THE SIMULATION
STUDY

In this section, we provide the omitted results for the statistics ν1 and ν2 defined in Section
6 of the paper. The results obtained with the Normal copula are given in Table 5. Figure 12
provides examples of matrices T̂ obtained with the different combinations of the factors
(T, n). The results obtained with the Cauchy copula are given in Table 6.

24

Table 5: Average values of ν1 and ν2 over the 500 simulations, as defined in Section 6,
from multivariate Normal distribution for all combinations of factors (T, n) and shrinkage
intensities w ∈ {0, 0.25, 0.5, 0.75, 1}.
w
T

T1

T2

T3

T4

n

w

0

.25

.5

.75

1

0

.25

.5

.75

1

125

0.00

0.94

0.95

0.94

0.92

0.00

0.63

0.63

0.63

0.62

250

0.02

1.00

1.00

1.00

0.99

0.01

0.65

0.65

0.65

0.65

500

1.00

1.00

1.00

1.00

1.00

0.65

0.65

0.65

0.65

0.65

125

0.00

0.41

0.50

0.52

0.50

0.00

0.64

0.67

0.68

0.68

250

0.20

0.94

0.95

0.95

0.95

0.27

0.78

0.79

0.79

0.79

500

1.00

1.00

1.00

1.00

1.00

0.79

0.79

0.79

0.79

0.79

125

0.00

0.21

0.26

0.29

0.23

0.00

0.47

0.49

0.50

0.48

250

0.07

0.76

0.81

0.81

0.75

0.08

0.58

0.59

0.59

0.57

500

0.98

0.99

0.99

0.99

0.98

0.62

0.62

0.62

0.63

0.62

125

0.00

0.00

0.00

0.00

0.00

0.00

0.43

0.48

0.50

0.51

250

0.00

0.00

0.01

0.00

0.00

0.10

0.47

0.50

0.51

0.50

500

0.03

0.08

0.10

0.11

0.08

0.49

0.53

0.55

0.56

0.53

ν̄1

ν̄2

25

∆

T

n = 125

n = 250

n = 500

T̂1

T̂2

T̂3

T̂4

Figure 12: Examples of matrices T̂ obtained from data generated using the 12 different
combinations (T, n) and the Normal copula.

26

Table 6: Average values of ν1 and ν2 over the 500 simulations, as defined in Section 6,
from a Cauchy copula with uniform (0, 1) margins for all combinations of factors (T, n)
and shrinkage intensities w ∈ {0, 0.25, 0.5, 0.75, 1}.
w
T

T1

T2

T3

T4

n

w

0

.25

.5

.75

1

0

.25

.5

.75

1

125

0.00

0.68

0.75

0.73

0.71

0.00

0.62

0.63

0.63

0.62

250

0.00

0.96

0.97

0.97

0.96

0.00

0.67

0.68

0.68

0.68

500

1.00

1.00

1.00

1.00

1.00

0.70

0.70

0.70

0.70

0.70

125

0.00

0.09

0.12

0.16

0.13

0.01

0.57

0.61

0.63

0.63

250

0.02

0.57

0.61

0.61

0.61

0.19

0.73

0.74

0.74

0.74

500

0.91

0.96

0.96

0.97

0.97

0.80

0.81

0.81

0.81

0.81

125

0.00

0.05

0.04

0.04

0.03

0.00

0.44

0.46

0.47

0.45

250

0.00

0.32

0.38

0.39

0.35

0.06

0.52

0.54

0.54

0.52

500

0.78

0.88

0.89

0.89

0.86

0.61

0.63

0.63

0.63

0.62

125

0.00

0.00

0.00

0.00

0.00

0.00

0.49

0.53

0.54

0.55

250

0.00

0.00

0.00

0.00

0.00

0.09

0.47

0.50

0.53

0.54

500

0.00

0.01

0.01

0.01

0.00

0.43

0.51

0.53

0.55

0.54

ν̄1

ν̄2

27

D

DETAILED CLUSTERS OF THE APPLICATION
TO STOCKS RETURNS

For each cluster, company names, tickers and sectors (according to nasdaq.com) are given.
A more detailed description of the clusters including the industry can be generated with
the R-code available online.
[1] "##################################################################"
[1] "######### cluster 1 ################### tau.tilde = 0.194 #######"
[1] "##################################################################"
Symbol

Name

Sector

1:

AAL

American Airlines Group, Inc.

Transportation

2:

ADP Automatic Data Processing, Inc.

Technology

3:

ESRX Express Scripts Holding Company

Health Care

4:

HAS

Hasbro, Inc. Consumer Non-Durables

5:

PAYX

Paychex, Inc.

Consumer Services

6:

SBUX

Starbucks Corporation

Consumer Services

7:

XRAY

DENTSPLY SIRONA Inc.

Health Care

[1] ""
[1] ""
[1] "##################################################################"
[1] "######### cluster 2 ################### tau.tilde = 0.246 #######"
[1] "##################################################################"
Symbol

Name

Sector

1:

AAPL

Apple Inc.

Technology

2:

ATVI

Activision Blizzard, Inc

Technology

3:

CERN

Cerner Corporation

Technology

28

4:

CHKP Check Point Software Technologies Ltd.

Technology

5:

EA

Electronic Arts Inc.

Technology

6:

EBAY

eBay Inc.

Miscellaneous

7:

EXPE

Expedia, Inc. Consumer Services

8:

HOLX

Hologic, Inc.

Health Care

9:

ISRG

Intuitive Surgical, Inc.

Health Care

10:

MU

Micron Technology, Inc.

Technology

11:

NVDA

NVIDIA Corporation

Technology

12:

IDXX

IDEXX Laboratories, Inc.

Health Care

[1] ""
[1] ""
[1] "##################################################################"
[1] "######### cluster 3 ################### tau.tilde = 0.376 #######"
[1] "##################################################################"
Symbol

Name

Sector

1:

ADBE Adobe Systems Incorporated

Technology

2:

AMZN

3:

FB

Facebook, Inc.

Technology

4:

MSFT

Microsoft Corporation

Technology

5:

NFLX

6:

PCLN

Amazon.com, Inc. Consumer Services

Netflix, Inc. Consumer Services
The Priceline Group Inc.

Miscellaneous

[1] ""
[1] ""
[1] "##################################################################"
[1] "######### cluster 4 ################### tau.tilde = 0.444 #######"
[1] "##################################################################"
Symbol

Name

29

Sector

1:

ADI

Analog Devices, Inc.

Technology

2:

AMAT

Applied Materials, Inc.

Technology

3:

AVGO

Broadcom Limited

Technology

4:

KLAC

5:

MXIM

Maxim Integrated Products, Inc.

Technology

6:

SWKS

Skyworks Solutions, Inc.

Technology

7:

TXN

Texas Instruments Incorporated

Technology

8:

LRCX

Lam Research Corporation

Technology

9:

MCHP Microchip Technology Incorporated

Technology

XLNX

Technology

10:

KLA-Tencor Corporation Capital Goods

Xilinx, Inc.

[1] ""
[1] ""
[1] "##################################################################"
[1] "######### cluster 5 ################### tau.tilde = 0.27 #######"
[1] "##################################################################"
Symbol

Name

Sector

1:

ADSK

Autodesk, Inc.

Technology

2:

CA

CA Inc.

Technology

3:

CTAS

4:

CSCO

Cisco Systems, Inc.

Technology

5:

CTXS

Citrix Systems, Inc.

Technology

6:

CTSH Cognizant Technology Solutions Corporation

Technology

7:

FISV

Fiserv, Inc.

Technology

8:

HSIC

Henry Schein, Inc.

Health Care

9:

MAR

Marriott International

Consumer Services

10:

PYPL

PayPal Holdings, Inc.

Miscellaneous

11:

QCOM

QUALCOMM Incorporated

Technology

Cintas Corporation Consumer Non-Durables

30

12:

VRSK

Verisk Analytics, Inc.

Technology

[1] ""
[1] ""
[1] "##################################################################"
[1] "######### cluster 6 ################### tau.tilde = 0.172 #######"
[1] "##################################################################"
Symbol

Name

Sector

1:

AKAM

Akamai Technologies, Inc.

Miscellaneous

2:

CSX

CSX Corporation

Transportation

3:

INTC

Intel Corporation

Technology

4:

INTU

Intuit Inc.

Technology

5:

MELI

MercadoLibre, Inc.

Miscellaneous

6:

NCLH Norwegian Cruise Line Holdings Ltd. Consumer Services

7:

SIRI

8:

SYMC

Symantec Corporation

Technology

9:

TMUS

T-Mobile US, Inc.

Public Utilities

10:

WDC

Western Digital Corporation

Technology

11:

STX

Seagate Technology PLC

Technology

12:

WYNN

Sirius XM Holdings Inc. Consumer Services

Wynn Resorts, Limited Consumer Services

[1] ""
[1] ""
[1] "##################################################################"
[1] "######### cluster 7 ################### tau.tilde = 0.294 #######"
[1] "##################################################################"
Symbol

Name

Sector

1:

ALXN

Alexion Pharmaceuticals, Inc.

Health Care

2:

AMGN

Amgen Inc.

Health Care

31

3:

BIIB

Biogen Inc.

Health Care

4:

BMRN

BioMarin Pharmaceutical Inc.

Health Care

5:

CELG

Celgene Corporation

Health Care

6:

ILMN

7:

INCY

Incyte Corporation

Health Care

8:

REGN

Regeneron Pharmaceuticals, Inc.

Health Care

9:

VRTX Vertex Pharmaceuticals Incorporated

Health Care

Illumina, Inc. Capital Goods

[1] ""
[1] ""
[1] "##################################################################"
[1] "######### cluster 8 ################### tau.tilde = 0.243 #######"
[1] "##################################################################"
Symbol

Name

Sector

1:

BIDU

Baidu, Inc.

Technology

2:

CTRP Ctrip.com International, Ltd.

Miscellaneous

3:

JD

JD.com, Inc. Consumer Services

4:

NTES

NetEase, Inc.

Miscellaneous

5:

TSLA

Tesla, Inc.

Capital Goods

[1] ""
[1] ""
[1] "##################################################################"
[1] "######### cluster 9 ################### tau.tilde = 0.12 #######"
[1] "##################################################################"
Symbol

Name

Sector

1:

CHTR

Charter Communications, Inc.

Consumer Services

2:

CMCSA

Comcast Corporation

Consumer Services

3:

COST

Costco Wholesale Corporation

Consumer Services

32

4:

DISH

DISH Network Corporation

Consumer Services

5:

FAST

Fastenal Company

Consumer Services

6:

JBHT J.B. Hunt Transport Services, Inc.

Transportation

7:

KHC

The Kraft Heinz Company Consumer Non-Durables

8:

QVCA

9:

MAT

Mattel, Inc. Consumer Non-Durables

10:

MDLZ

Mondelez International, Inc. Consumer Non-Durables

11:

MNST

Monster Beverage Corporation Consumer Non-Durables

12:

ULTA

Ulta Beauty, Inc.

Consumer Services

13:

VIAB

Viacom Inc.

Consumer Services

14:

VOD

Vodafone Group Plc

Public Utilities

15:

WBA

Walgreens Boots Alliance, Inc.

Health Care

16:

PCAR

PACCAR Inc.

Capital Goods

Liberty Interactive Corporation

Consumer Services

[1] ""
[1] ""
[1] "##################################################################"
[1] "######### cluster 10 ################### tau.tilde = 0.876 #######"
[1] "##################################################################"
Symbol

Name

Sector

1:

DISCA Discovery Communications, Inc. Consumer Services

2:

DISCK Discovery Communications, Inc. Consumer Services

[1] ""
[1] ""
[1] "##################################################################"
[1] "######### cluster 11 ################### tau.tilde = 0.338 #######"
[1] "##################################################################"
Symbol

Name

Sector

33

1:

DLTR

Dollar Tree, Inc. Consumer Services

2:

ROST

Ross Stores, Inc. Consumer Services

3:

TSCO

Tractor Supply Company Consumer Services

4:

ORLY O&#39;Reilly Automotive, Inc. Consumer Services

[1] ""
[1] ""
[1] "##################################################################"
[1] "######### cluster 12 ################### tau.tilde = 0.874 #######"
[1] "##################################################################"
Symbol

Name

Sector

1:

FOX Twenty-First Century Fox, Inc. Consumer Services

2:

FOXA Twenty-First Century Fox, Inc. Consumer Services

[1] ""
[1] ""
[1] "##################################################################"
[1] "######### cluster 13 ################### tau.tilde = 0.201 #######"
[1] "##################################################################"
Symbol
1:

Name

Sector

GILD Gilead Sciences, Inc. Health Care

2:

MYL

Mylan N.V. Health Care

3:

SHPG

Shire plc Health Care

[1] ""
[1] ""
[1] "##################################################################"
[1] "######### cluster 14 ################### tau.tilde = 0.893 #######"
[1] "##################################################################"
Symbol

Name

Sector

34

1:

GOOG Alphabet Inc. Technology

2:

GOOGL Alphabet Inc. Technology

[1] ""
[1] ""
[1] "##################################################################"
[1] "######### cluster 15 ################### tau.tilde = 0.85 #######"
[1] "##################################################################"
Symbol

Name

Sector

1:

LBTYK Liberty Global plc Consumer Services

2:

LBTYA Liberty Global plc Consumer Services

[1] ""
[1] ""
[1] "##################################################################"
[1] "######### cluster 16 ################### tau.tilde = 0.757 #######"
[1] "##################################################################"
Symbol

Name

Sector

1:

LILA Liberty Global plc Consumer Services

2:

LILAK Liberty Global plc Consumer Services

35

References
Agarwal, A., Negahban, S., and Wainwright, M. J. (2012). Noisy matrix decomposition
via convex relaxation: Optimal rates in high dimensions. The Annals of Statistics, pages
1171–1197.
Ben Ghorbal, N., Genest, C., and Nešlehová, J. (2009). On the Ghoudi, Khoudraji, and
Rivest test for extreme-value dependence. Canadian Journal of Statistics, 37:534–552.
Borkowf, C. B. (2002). Computing the nonnull asymptotic variance and the asymptotic
relative efficiency of Spearman’s rank correlation. Comput. Statist. Data Anal., 39:271–
286.
Brechmann, E. C. (2014). Hierarchical Kendall copulas: Properties and inference. Canadian
Journal of Statistics, 42:78–108.
Cai, T. T. and Zhang, L. (2017). High-dimensional Gaussian copula regression: Adaptive
estimation and statistical inference. Statistica Sinica, forthcomming.
Chandrasekaran, V., Parrilo, P. A., and Willsky, A. S. (2010). Latent variable graphical
model selection via convex optimization. In Communication, Control, and Computing
(Allerton), 2010 48th Annual Allerton Conference on, pages 1610–1613. IEEE.
Datta, A. and Zou, H. (2017). Cocolasso for high-dimensional error-in-variables regression.
Ann. Statist., forthcomming.
Devlin, S. J., Gnanadesikan, R., and Kettenring, J. R. (1975). Robust estimation and
outlier detection with correlation coefficients. Biometrika, 62(3):531–545.
Ehrenberg, A. (1952). On sampling from a population of rankers. Biometrika, 39:82–87.
36

El Maache, H. and Lepage, Y. (2003). Spearman’s rho and Kendall’s tau for multivariate data sets. In Mathematical statistics and applications: Festschrift for Constance
van Eeden, volume 42 of IMS Lecture Notes Monogr. Ser., pages 113–130. Institute of
Mathematical Statistics, Beachwood, OH.
Elton, E. J. and Gruber, M. J. (1973). Estimating the dependence structure of share prices
—implications for portfolio selection. The Journal of Finance, 28:1203–1232.
Embrechts, P., McNeil, A. J., and Straumann, D. (2002). Correlation and dependence
in risk management: Properties and pitfalls. In Risk Management: Value at Risk and
Beyond (Cambridge, 1998), pages 176–223. Cambridge University Press, Cambridge.
Engle, R. and Kelly, B. (2012). Dynamic equicorrelation. Journal of Business & Economic
Statistics, 30(2):212–228.
Fang, H.-B., Fang, K.-T., and Kotz, S. (2002). The meta-elliptical distributions with given
marginals. Journal of Multivariate Analysis, 82:1 – 16.
Fang, K.-T., Kotz, S., and Ng, K. W. (1990). Symmetric Multivariate and Related Distributions. Chapman & Hall, London.
Fang, K.-T. and Zhang, Y.-T. (1990). Generalized Multivariate Analysis. Springer Verlag
& Science Press, Beijing.
Genest, C. and Nešlehová, J. (2012). Copulas and copula models. In El-Shaarawi, A. H.
and Piegorsch, W. W., editors, Encyclopedia of Environmetrics. Wiley, Chichester, 2nd
edition.
Genest, C., Nešlehová, J., and Ben Ghorbal, N. (2011). Estimators based on Kendall’s

37

tau in multivariate copula models. Australian & New Zealand Journal of Statistics,
53:157–177.
Gregory, J. and Laurent, J.-P. (2004). In the core of correlation. Risk, 17:87–91.
Harville, D. A. (2008). Matrix Algebra From a Statistician’s Perspective. Springer Science
& Business Media.
Hoeffding, W. (1947). On the distribution of the rank correlation coefficient τ when the
variates are not independent. Biometrika, 34:183–196.
Hoeffding, W. (1948). A class of statistics with asymptotically normal distribution. The
Annals of Mathematical Statistics, 19:293–325.
Hua, L. and Joe, H. (2017). Multivariate dependence modeling based on comonotonic
factors. Journal of Multivariate Analysis, 155:317–333.
Hult, H. and Lindskog, F. (2002). Multivariate extremes, aggregation and dependence in
elliptical distributions. Advances in Applied Probability, 34:587–608.
Joe, H. (2015). Dependence Modeling With Copulas. CRC Press, Boca Raton, FL.
Krupskii, P. and Joe, H. (2013). Factor copula models for multivariate data. Journal of
Multivariate Analysis, 120:85–101.
Krupskii, P. and Joe, H. (2015). Structured factor copula models: Theory, inference
and computation. Journal of Multivariate Analysis, 138(Supplement C):53 – 73. HighDimensional Dependence and Copulas.
Kurowicka, D. and Joe, H., editors (2011). Dependence Modeling: Handbook on Vine
Copulae. World Scientific Publishing, Hackensack, NJ.
38

Ledoit, O. and Wolf, M. (2003a). Honey, I shrunk the sample covariance matrix. UPF
economics and business working paper.
Ledoit, O. and Wolf, M. (2003b). Improved estimation of the covariance matrix of stock
returns with an application to portfolio selection. Journal of Empirical Finance, 10:603–
621.
Ledoit, O. and Wolf, M. (2004). A well-conditioned estimator for large-dimensional covariance matrices. Journal of Multivariate Analysis, 88:365–411.
Lindeberg, J. (1927). Über die korrelation. In Den VI skandinaviske Matematikerkongres
i København, pages 437–446, Copenhagen, Denmark. J. Gjellerup.
Lindeberg, J. (1929). Some remarks on the mean error of the percentage of correlation.
Nordic Statistical Journal, 1:137–141.
Lindskog, F., McNeil, A. J., and Schmock, U. (2002). Kendall’s tau for elliptical distributions. In Bol, G., Nakhaeizadeh, G., Rachev, S. T., Ridder, T., and Vollmer,
K.-H., editors, Credit Risk: Measurement, Evaluation and Management, pages 149–156.
Springer.
Liu, H., Han, F., Yuan, M., Lafferty, J., and Wasserman, L. (2012). High-dimensional
semiparametric gaussian copula graphical models. The Annals of Statistics, 40:2293–
2326.
Ma, S., Xue, L., and Zou, H. (2013). Alternating direction methods for latent variable
gaussian graphical model selection. Neural Computation, 25(8):2172–2198.
Mai, J.-F. and Scherer, M. (2012). H-extendible copulas. Journal of Multivariate Analysis,
110:151—160.
39

McNeil, A. J., Frey, R., and Embrechts, P. (2015). Quantitative risk management. Princeton
Series in Finance. Princeton University Press, Princeton, NJ, revised edition. Concepts,
techniques and tools.
Nelsen, R. B. (2006). An Introduction to Copulas. Springer, New York, 2nd edition.
Patton, A. J. (2006). Modelling asymmetric exchange rate dependence. Internat. Econom.
Rev., 47:527–556.
Patton, A. J. (2012). A review of copula models for economic time series. Journal of
Multivariate Analysis, 110(Supplement C):4 – 18. Special Issue on Copula Modeling and
Dependence.
Rémillard, B. (2013). Statistical methods for financial engineering. CRC Press, Boca Raton,
FL.
Rémillard, B. (2017). Goodness-of-fit tests for copulas of multivariate time series. Econometrics, 5:13.
Rousseeuw, P. J. and Molenberghs, G. (1993). Transformation of non positive semidefinite
correlation matrices. Communications in Statistics - Theory and Methods, 22:965–984.
Schäfer, J. and Strimmer, K. (2005). A shrinkage approach to large-scale covariance matrix
estimation and implications for functional genomics. Statistical applications in genetics
and molecular biology, 4:32.
Severini, T. A. (2005). Elements of Distribution Theory. Cambridge University Press.
Sklar, A. (1959). Fonctions de répartition à n dimensions et leurs marges. Publications de
l’Institut de Statistique de l’Université de Paris, 8:229–231.
40

Stewart, G. W. (1969). On the continuity of the generalized inverse. SIAM Journal on
Applied Mathematics, 17:33–45.
Xue, L. and Zou, H. (2012). Regularized rank-based estimation of high-dimensional nonparanormal graphical models. The Annals of Statistics, 40(5):2541–2571.
Zhao, Y. and Genest, C. (2017). Inference for elliptical copula multivariate response regression. Technical report, McGill University, Montréal, Canada.

NOTATION INDEX
d – number of random variables considered: X1 , . . . , Xd ;
T – d × d Kendall τ matrix of X1 , . . . , Xd ;
τ – vectorized version of T;
p – dimension of τ , p := d(d − 1)/2;
τ̂ – usual estimator of τ ;
G – partition of {1, . . . , d} defining the block structure of T or, equivalently, the clustering of X1 , . . . , Xd ;
K – number of clusters: K = |G|;
Gk – clusters forming the partition G, that is G = {Gk , . . . , GK };
BG – partition of τ defining the blocks of T;
L – number of distinct blocks in T, L = |BG | ≤ K(K + 1)/2;
41

Bk – blocks forming the partition BG , that is BG = {B1 , . . . , BL };
TG – space of matrices clustered according to BG , in particular T ∈ TG ;
Σ – covariance matrix of τ̂ ;
Σ̂ – basic estimator of Σ;
CG – partition of Σ;
SG – space of matrices clustered according to CG , in particular Σ ∈ SG .

42

