A test for separability in covariance operators of
random surfaces
Pramita Bagchi, Holger Dette
Ruhr-Universität Bochum

arXiv:1710.08388v2 [stat.ME] 13 Nov 2017

Fakultät für Mathematik
44780 Bochum
Germany
November 15, 2017

Abstract
The assumption of separability is a simplifying and very popular assumption in
the analysis of spatio-temporal or hypersurface data structures. It is often made in
situations where the covariance structure cannot be easily estimated, for example
because of a small sample size or because of computational storage problems. In
this paper we propose a new and very simple test to validate this assumption. Our
approach is based on a measure of separability which is zero in the case of separability and positive otherwise. The measure can be estimated without calculating
the full non-separable covariance operator. We prove asymptotic normality of the
corresponding statistic with a limiting variance, which can easily be estimated from
the available data. As a consequence quantiles of the standard normal distribution
can be used to obtain critical values and the new test of separability is very easy to
implement. In particular, our approach does neither require projections on subspaces generated by the eigenfunctions of the covariance operator, nor resampling
procedures to obtain critical values nor distributional assumptions as recently used
by Aston et al. (2017) and Constantinou et al. (2017) to construct tests for separability. We investigate the finite sample performance by means of a simulation study
and also provide a comparison with the currently available methodology. Finally,
the new procedure is illustrated analyzing wind speed and temperature data.

Keywords: functional data, minimum distance, separability, space-time processes, surface data structures
AMS Subject classification: 62G10, 62G20

1

1 Introduction
Data, which is functional and multidimensional is usually called surface data and arises
in areas such as medical imaging [see Skup (2010); Worsley et al. (1996)], spectrograms
derived from audio signals or geolocalized data [see Bar-Hen et al. (2008); Rabiner and
Schafer (1978) ]. In many of these ultra high-dimensional problems a completely nonparametric estimation of the covariance operator is not possible as the number of available observations is small compared to the dimension. A common approach to obtain
reasonable estimates in this context are structural assumptions on the covariance of the
underlying process, and in recent years the assumption of separability has become very
popular, for example in the analysis of geostatistical space-time models [see Genton
(2007); Gneiting et al. (2007), among others]. Roughly speaking, this assumption allows
to write the covariance
c(s, t , s 0 , t 0 ) = E[X (s, t )X (s 0 , t 0 )]
of a (real valued) space-time process {X (s, t )}(s,t )∈S×T as a product of the space and time
covariance function, that is
c(s, t , s 0 , t 0 ) = c 1 (s, s 0 )c 2 (t , t 0 ).

(1.1)

It was pointed out by many authors that the assumption of separability yields a substantial simplification of the estimation problem and thus reduces computational costs in
the estimation of the covariance in high dimensional problems [see for example Huizenga
et al. (2002); Rougier (2017)]. Despite of its importance, there exist only a few tools to
validate the assumption of separability for surface data.
Many authors developed tests for spatio-temporal data. For example, Fuentes (2006)
proposed a test based on the spectral representation, and Lu and Zimmerman (2005);
Mitchell et al. (2005, 2006) investigated likelihood ratio tests under the assumption of a
normal distribution. Recently, Constantinou et al. (2017) derived the joint distribution
of the three statistics appearing in the likelihood ratio test and used this result to derive
the asymptotic distribution of the (log) likelihood ratio. These authors also proposed
alternative tests which are based on distances between an estimator of the covariance
under the assumption of separability and an estimator which does not use this assumption.
Aston et al. (2017) considered the problem of testing for separability in the context of
hypersurface data. These authors pointed out that many available methods require the
estimation of the full multidimensional covariance structure, which can become infeasible for high dimensional data. In order to address this issue they developed a bootstrap
test for applications, where replicates from the underlying random process are available.
To avoid estimation and storage of the full covariance finite-dimensional projections of

2

the difference between the covariance operator and a nonparametric separable approximation are proposed. In particular they suggest to project onto subspaces generated
by the eigenfunctions of the covariance operator estimated under the assumption of
separability. However, as pointed in the same references the choice of the number of
eigenfunctions onto which one should project is not trivial and the test might be sensitive with respect to this choice.
In this paper we present an alternative and very simple test for the hypothesis of
separability in hypersurface data. We consider a similar setup as in Aston et al. (2017)
and proceed in two steps. First we derive an explicit expression for the minimal distance between the covariance operator and its approximation by a separable covariance
operator, where the minimum is taken with respect to the second factor of the tensor
product. It turns out that this minimum vanishes if and only if the covariance operator
is separable. Secondly, we directly estimate the minimal distance (and not the covariance operator itself ) from the available data. As a consequence the calculation of the
test statistic does neither use an estimate of the full non-separable covariance operator
nor requires the specification of subspaces used for a projection. In the main result of
this paper we derive the asymptotic distribution of the test statistic, which is normal
(after appropriate standardization) under the null hypothesis and alternative. The limiting variance under the null hypothesis can easily be estimated and as consequence we
obtain a very simple test for the hypothesis of separability, which only requires the quantiles of the normal distribution. Moreover, in contrast to the work of Aston et al. (2017),
the test proposed here does not require resampling and - from a theoretical perspective - the limiting theorems are valid under more general and easier to verify moment
assumptions.
In Section 2 we review some basic terminology and minimize the distance between
the covariance operator and separable covariance operators with respect to the second
factor of the tensor product. This minimum distance could also be interpreted as a
measure of deviation from separability (it is zero in the case of separability and positive otherwise). In Section 3 we propose an estimator of the minimum distance and
prove asymptotic normality of a standardised version of this statistic. We also provide
a simple estimate of the limiting variance under the null hypothesis and prove its consistency. Section 4 is devoted to an investigation of the finite sample properties of the
new test and a comparison with two alternative tests for this problem, which have recently been proposed by Aston et al. (2017) and Constantinou et al. (2017). In particular
we demonstrate that - despite of its simplicity - the new procedure has very competitive
properties compared to the currently available methodology. Finally, some technical
details are deferred to the Appendix A.

3

2 Hilbert spaces and a measure of separability
We begin introducing some basic facts about Hilbert spaces, Hilbert-Schmidt operators
and tensor products. For more details we refer to the monographs of Weidmann (1980),
Dunford and Schwartz (1988) or Gohberg et al. (1990). Let H be a real separable Hilbert
space with inner-product 〈·, ·〉 and norm k · k. The space of bounded linear operators on
H is denoted by S ∞ (H ) with operator norm
T ∞ := sup kT f k.
k f k≤1

A bounded linear operator T is said to be compact if it can be written as
T=

X

s j (T )〈e j , ·〉 f j ,

j ≥1

where {e j : j ≥ 1} and { f j : j ≥ 1} are orthonormal sets of H , {s j (T ) : j ≥ 1} are the singular values of T and the series converges in the operator norm. We say that a compact
operator T belongs to the Schatten class of order p ≥ 1 and write T ∈ S p (H ) if
T p =

³X

s j (A)

´1/p

< ∞.

j ≥1

The Schatten class of order p ≥ 1 is a Banach space with norm .p and with the property
that S p (H ) ⊂ S q (H ) for p < q. In particular we are interested in Schatten classes of order
p = 1 and 2. A compact operator T is called Hilbert-Schmidt operator if T ∈ S 2 (H ) and
trace class if T ∈ S 1 (H ). The space of Hilbert-Schmidt operators S 2 (H ) is also a Hilbert
space with the Hilbert-Schmidt inner product given by
〈A, B 〉H S =

X

〈Ae j , B e j 〉

j ≥1

for each A, B ∈ S 2 (H ), where {e j : j ≥ 1} is an orthonormal basis (the inner product does
not depend on the choice of the basis).
For two real separable Hilbert spaces H1 and H2 , the tensor product of H1 and H2 ,
denoted as H := H1 ⊗ H2 , is the Hilbert space obtained by the completion of all finite
sums

N
X

u i ⊗ v j , u i ∈ H1 , v j ∈ H2 ,

i , j =1

under the inner product 〈u ⊗ v, w ⊗ z〉 = 〈u, w〉〈v, z〉, for u, w ∈ H1 and v, z ∈ H2 . For
˜ 2 is defined as the unique linear
C 1 ∈ S ∞ (H1 ) and C 2 ∈ S ∞ (H2 ), the tensor product C 1 ⊗C
operator on H1 ⊗ H2 satisfying
˜ 2 )(u ⊗ v) = C 1 u ⊗C 2 v, u ∈ H1 , v ∈ H2 .
(C 1 ⊗C
4

˜ 2 ∈ S ∞ (H ) with C 1 ⊗C
˜ 2 ∞ = C 1 ∞ C 2 ∞ . Moreover, if C 1 ∈ S p (H1 ) and
In fact C 1 ⊗C
˜ 2 ∈ S p (H ) with C 1 ⊗C
˜ 2 p = C 1 p C 2 p . For more
C 2 ∈ S p (H2 ), for p ≥ 1, then C 1 ⊗C
details we refer to Chapter 8 of Weidmann (1980). In the sequel we denote the HilbertSchmidt inner product on S 2 (H ) with H = H1 ⊗ H2 as 〈·, ·〉H S and that of S 2 (H1 ) and
S 2 (H2 ) as 〈·, ·〉S 2 (H1 ) and 〈·, ·〉S 2 (H2 ) respectively.

2.1 Measuring separability
We consider random elements X in the Hilbert space H with EkX k4 < ∞. (See Chapter
7 from Hsing and Eubank (2015)) Then the covariance operator of X is defined as C :=
E [(X − EX ) ⊗o (X − EX )], where for f , g ∈ H the operator f ⊗o g : H → H is defined by
( f ⊗o g )h = 〈h, g 〉 f for all h ∈ H .
Under the assumption EkX k4 < ∞ we have C ∈ S 2 (H ). We also assume C 2 6= 0, which
essentially means the random variable X is not degenerate. To test separability we consider the hypothesis
H0 : C = C 1 ⊗C 2 for some C 1 ∈ S 2 (H1 ) and C 2 ∈ S 2 (H2 ).

(2.1)

Our approach is based on a approximation of the operator C by a separable operator
˜ 2 with respect to the norm  · 2 . Ideally, we are looking for
C 1 ⊗C
¯
n
o
˜ 2 22 ¯¯ C 1 ∈ S 2 (H1 ),C 2 ∈ S 2 (H2 ) ,
D := inf C −C 1 ⊗C

(2.2)

such that the hypothesis of separability in (2.1) can be rewritten in terms of the distance
D, that is
H0 : D = 0 versus H1 : D > 0 .

(2.3)

However, it turns out that it is difficult to express D explicitly in terms of the covariance
operator C . For this reason we proceed in a slightly different way in two steps. First we
fix C 1 and determine
¯
n
o
˜ 2 22 ¯¯ C 2 ∈ S 2 (H2 ) ,
D(C 1 ) := inf C −C 1 ⊗C

(2.4)

˜ 2 22 with respect to second factor C 2 of the tensor
that is we are minimizing C − C 1 ⊗C
product. In particular we will show that the infimum is in fact a minimum and derive
an explicit expression for D(C 1 ) and its minimizer. Next we shows that the resulting
minimum, D(C 1 ) vanishes if and only if the hypothesis of separability holds.
For this purpose we have to introduce additional notation and have to prove several
auxiliary results. The main statement is given in Theorem 2.1 (whose formulation also

5

requires the new notation). First, consider the bounded linear operator T1 : S 2 (H ) ×
S 2 (H1 ) 7→ S 2 (H2 ) defined by
˜
T1 (A ⊗B,C
1 ) = 〈A,C 1 〉S 2 (H1 ) B

(2.5)

for all C 1 ∈ S 2 (H1 ). Similarly, let T2 : S 2 (H ) × S 2 (H2 ) → S 2 (H1 ) be the bounded linear
operator defined by
˜
T2 (A ⊗B,C
2 ) = 〈B,C 2 〉S 2 (H2 ) A

(2.6)

for all C 2 ∈ S 2 (H2 ).
Proposition 2.1. The operators T1 and T2 are well-defined, bi-linear and continuous with
˜ 〉H S ,
〈B, T1 (C ,C 1 )〉S 2 (H2 ) = 〈C ,C 1 ⊗B

(2.7)

˜ 2 〉H S .
〈A, T2 (C ,C 2 )〉S 2 (H1 ) = 〈C , A ⊗C

(2.8)

for all A,C 1 ∈ S 2 (H1 ), B,C 2 ∈ S 2 (H2 ) and C ∈ S 2 (H ).
Proof. By Lemma A.1 in Appendix A the space
D0 :=

nX
n

o
˜ i : A i ∈ S 2 (H1 ), B i ∈ S 2 (H2 ), n ∈ N
A i ⊗B

(2.9)

i =1

is a dense subset of S 2 (H1 ⊗ H2 ) (note that a similar result for the space S 1 (H1 ⊗ H2 ) has
been established in Lemma 1.6 of the supplementary material in Aston et al. (2017)). For
all B ∈ S 2 (H2 ), E ∈ D0 and C 1 ∈ S 2 (H1 ), we have
D X
E
n
〈B, T1 (E ,C 1 )〉S 2 (H2 ) = B, 〈A i ,C 1 〉S 2 (H1 ) B i

S 2 (H2 )

i =1

E
DX
n
˜ i ,C 1 ⊗B
˜
=
A i ⊗B
i =1

HS

n
X

=

〈A i ,C 1 〉S 2 (H1 ) 〈B, B i 〉S 2 (H2 )

i =1

˜ 〉H S .
= 〈E ,C 1 ⊗B

(2.10)

Using the fact that
©
ª
T1 (E ,C 1 )2 ≤ T1 (E ,C 1 )1 = sup 〈B, T1 (E ,C 1 )〉S 2 (H ) : B ∈ S 2 (H2 ), B ∞ ≤ 1 ,

(2.11)

(2.10) and the Cauchy Schwarz inequality it follows that
T1 (E ,C 1 )2 ≤ C 1 2 E 2 .

(2.12)

Therefore, for each C 1 ∈ S 2 (H1 ), we can extend T1 (.,C 1 ) continuously on S 2 (H ).
Furthermore as (2.7) holds for all C ∈ D0 and the maps C 7→ 〈B, T1 (C ,C 1 )〉S 2 (H1 ) and
˜ 〉H S are continuous for all B ∈ S 2 (H2 ) and C 1 ∈ S 2 (H1 ), we can conclude
C 7→ 〈C ,C 1 ⊗B
that (2.7) holds for all C ∈ S 2 (H ).

6

Corollary 2.1. For all C ∈ S 2 (H ), C 1 ∈ S 2 (H1 ) and C 2 ∈ S 2 (H2 ) we have
T1 (C ,C 1 )2 ≤ C 2 C 1 2 and T2 (C ,C 2 )2 ≤ C 2 C 2 2 .
Proposition 2.2. For any C ∈ S 2 (H ) and C 1 ∈ S 2 (H1 ), we have
˜ 1 (C ,C 1 )〉H S = T1 (C ,C 1 )22 .
〈C ,C 1 ⊗T
Proof. Recall the definition of the set D0 in (2.9) and let C =
A n ∈ S 2 (H1 ), B n ∈ S 2 (H2 ). We write
*

N
X

˜ 1 (C ,C 1 )〉H S = C ,C 1 ⊗
˜
〈C ,C 1 ⊗T

=

n=1

˜ n ∈ D0 , where
A n ⊗B

+
〈A n ,C 1 〉S 2 (H1 ) B n

n=1

=

PN

HS

N
X

˜ n 〉H S
〈A n ,C 1 〉S 2 (H1 ) 〈C ,C 1 ⊗B
n=1
N X
N
X

〈A n ,C 1 〉S 2 (H1 ) 〈A m ,C 1 〉S 2 (H1 ) 〈B m , B n 〉S 2 (H2 ) .

n=1 m=1

On the other hand,
*
〈T1 (C ,C 1 ), T1 (C ,C 1 )〉S 2 (H2 ) =

N
X

〈A n ,C 1 〉S 2 (H1 ) B n ,

n=1

=

N
X

N
X

+

N
X

〈A m ,C 1 〉S 2 (H1 ) B m

m=1

S 2 (H2 )

〈A n ,C 1 〉S 2 (H1 ) 〈A m ,C 1 〉S 2 (H1 ) 〈B m , B n 〉S 2 (H2 ) .

n=1 m=1

Therefore, for all C 1 ∈ S 2 (H1 ), the functions f , g : S 2 (H ) 7→ R defined by
˜ 1 (C ,C 1 )〉H S and g (C ) := T1 (C ,C 1 )22
f (C ) := 〈C ,C 1 ⊗T
are continuous and coincide on the dense subset D0 of S 2 (H ). So f (C ) = g (C ) for all
C ∈ S 2 (H ) and hence the result follows.
Theorem 2.1. For each C ∈ S 2 (H ) and C 1 ∈ S 2 (H1 ) the distance
˜ 2 2
D(C 1 ,C 2 ) = C −C 1 ⊗C
is minimized at
Ce2 =

T1 (C ,C 1 )
C 1 22

(2.13)

.

(2.14)

Moreover, for C 1 ∈ S 2 (H1 ) the minimum distance in (2.13) is given by
D(C 1 ) = C 22 −

T1 (C ,C 1 )22
C 1 22

.

˜ 2 for some C 2 ∈ S 2 (H2 ).
In particular D(C 1 ) is zero if and only if C = C 1 ⊗C
7

(2.15)

Proof. We write
˜ 2 22 =C −C 1 ⊗
˜ Ce2 22 + C 1 ⊗
˜ Ce2 −C 1 ⊗C
˜ 2 22
C −C 1 ⊗C
˜ Ce2 ,C 1 ⊗
˜ Ce2 −C 1 ⊗C
˜ 2 〉H S .
+ 2〈C −C 1 ⊗
For the last term we obtain from (2.14)
˜ Ce2 ,C 1 ⊗
˜ Ce2 −C 1 ⊗C
˜ 2 〉H S =〈C ,C 1 ⊗
˜ Ce2 〉H S − C 1 ⊗
˜ Ce2 22
〈C −C 1 ⊗
˜ 2 〉H S + 〈C 1 ⊗
˜ Ce2 ,C 1 ⊗C
˜ 2 〉H S
− 〈C ,C 1 ⊗C
=

1

˜ 1 (C ,C 1 )〉H S −
〈C ,C 1 ⊗T
2

T1 (C ,C 1 )22

C 1 2

C 1 22

˜ 2 〉H S + 〈C 2 , T1 (C ,C 1 )〉H S ,
− 〈C ,C 1 ⊗C
which is zero by (2.7) and Proposition 2.2. Therefore for all C 2 ∈ S 2 (H2 ) we have
˜ 2 22 ≥ C −C 1 ⊗
˜ Ce2 22
C −C 1 ⊗C
which proves the first assertion of Theorem 2.1.
For a proof of the representation (2.15) we substitute the operator Ce2 defined in (2.14)
for C 2 in the expression of D(C 1 ,C 2 ) and obtain
˜ Ce2 22 = 〈C −C 1 ⊗
˜ Ce2 ,C −C 1 ⊗
˜ Ce2 〉H S
D(C 1 ) =D(C 1 , Ce2 ) = C −C 1 ⊗
˜ Ce2 22 − 2〈C ,C 1 ⊗
˜ Ce2 〉H S
=C 22 + C 1 ⊗
=C 22 +

T1 (C ,C 1 )22
C 1 22

−

2
C 1 22

˜ 1 (C ,C 1 )〉H S .
〈C ,C 1 ⊗T

Then the second assertion follows from Proposition 2.2.
˜ 2 for some C 2 ∈ S(H2 ), then (2.5) implies
Now assume that C = C 1 ⊗C
T1 (C ,C 1 )22 = (〈C 1 ,C 1 〉S 2 (H1 ) )2 C 2 22 = C 1 42 C 2 22
˜ Ce2 , with Ce2 2 ≤
and therefore D(C 1 ) = 0. Conversely, if D(C 1 ) = 0, we have C = C 1 ⊗
C 2 by Corollary 2.1.

2.2 Hilbert-Schmidt integral operators
An important case for applications consists is the set H := L 2 (S × T, R) of all square integrable functions defined on S × T , where S ⊂ Rp , T ⊂ Rq are bounded measurable
sets. If the covariance operator C of the random variable X is a Hilbert-Schmidt operator it follows from Theorem 6.11 in Weidmann (1980) that there exists a kernel c ∈
¡
¢
L 2 (S × T ) × (S × T ), R such that C can be characterized as an integral operator, i.e.
Z Z
C f (s, t ) =
c(s, t , s 0 , t 0 ) f (s 0 , t 0 )d s 0 d t 0 , f ∈ L 2 (S × T, R),
S T

8

almost everywhere on S × T . Moreover the kernel is given by the covariance kernel of
£
¤
X , that is c(s, t , s 0 , t 0 ) = Cov X (s, t ), X (s 0 , t 0 ) , and the space H can be identified with the
tensor product of H1 = L 2 (S, R) and H2 = L 2 (T, R).
Similarly, if C 1 and C 2 are Hilbert-Schmidt operators on L 2 (S, R) and L 2 (T, R) respectively there exists symmetric kernels c 1 ∈ L 2 (H1 × H1 , R) and c 2 ∈ L 2 (H2 × H2 , R) such that
Z
C 1 f (s) =
c 1 (s, s 0 ) f (s 0 )d s 0 , f ∈ H1
ZS
C 2 g (t ) =
c 2 (t , t 0 )g (t 0 )d t 0 , g ∈ H2
T

almost everywhere on S and T , respectively. The following result shows that in this case
the operators T1 and T2 defined by (2.5) and (2.6), respectively, are also Hilbert-Schmidt
(or equivalently integral) operators.
¡
¢
Proposition 2.3. If C and C 1 are integral operators with kernels c ∈ L 2 (S ×T )×(S ×T ), R
and c 1 ∈ L 2 (S × S, R), then T1 (C ,C 1 ) is a an integral operator with kernel given by
Z Z
0
k(t , t ) =
c(s, t , s 0 t 0 )c 1 (s, s 0 )d sd s 0 .

(2.16)

S S

An analog result is true for the operator T2 .
Proof. By Lemma A.2 in Appendix A for a given ² > 0 there exists an integral operator C 0
with kernel c 0 such that
C −C 0 2 <
where C 0 =

²
and kc − c 0 k∞ < ²/2,
2C 1 2

PN

˜ n , and A n and B n are finite rank operators with continuous kerA n ⊗B
PN
nels a n and b n . Note that n=1
a n b n is the kernel of the operator C 0 . Let K be the integral
n=1

operator with the kernel defined in (2.16) and K 0 be the integral operator with kernel
Z Z
0
0
k (t , t ) :=
c 0 (s, t , s 0 t 0 )c 1 (s, s 0 )d sd s 0 ,
S S

then (note that K 0 is a Hilbert-Schmidt operator)
T1 (C ,C 1 ) − K 2 ≤T1 (C ,C 1 ) − T1 (C 0 ,C 1 )2 + T1 (C 0 ,C 1 ) − K 0 2 + K 0 − K 2 .

(2.17)

By (2.12) the first term is bounded by C 1 2 C − C 1 2 < ²/2. To handle the second term
note that for any f ∈ H2 ,
Ã
T1 (C 0 ,C 1 ) f (t ) =T1

N
X

!
˜ n ,C 1 f (t ) =
A n ⊗B

n=1

=

〈A n ,C 1 〉S 2 (H1 ) B n f (t )

n=1

N Z Z Z
X
n=1 T

N
X

S S

a n (s, s 0 )c 1 (s, s 0 )b n (t , t 0 ) f (t 0 )d sd s 0 d t 0
9

Z Z Z
=
=

ZT
T

S S

c 0 (s, t , s 0 , t 0 )c 1 (s, s 0 )d sd s 0 f (t 0 )d t 0

k 0 (t , t 0 ) f (t 0 )d t 0 = K 0 f (t ).

Therefore, the second term in (2.17) is zero. If |S| and |T | denote the Lebesgue measure
of the set S and T , respectively, the third term can be written as
·Z Z µZ Z
T

T

S S

0

0

0

0 0

0

(c(s, t , s , t ) − c (s, t , s t ))c 1 (s, s )d sd s

0

¶2
dtdt

0

¸1/2

≤ kc − c 0 k∞ kc 1 k2 |S|2 |T |,
which is bounded by kc 1 k2 |S|2 |T |²/2. Since the choice of ² > 0 is arbitrary, this proves
the assertion of Proposition 2.3.
Using the explicit formula for T1 described in Proposition 2.3 the minimum distance
can be expressed in terms of the corresponding kernels of the operators, that is
D(C 1 ) =

C 22 −

T1 (C ,C 1 )22
C 1 22

Z Z Z Z

=

c 2 (s, t , s 0 , t 0 )d sd s 0 d t d t 0
¤
R R £R R
0 0
0
0 2
dtdt0
S S c(s, t , s t )c 1 (s, s )d sd s
T T
R R 2
−
.
0
0
S S c 1 (s, s )d sd s
T

T

(2.18)

S S

3 Estimation and asymptotic properties
We estimate the minimum distance given in (2.15) by plugging in estimators for C and
C 1 based on a sample X 1 , X 2 , . . . , X N . The covariance operator C is estimated by
CbN :=

i
N h
1 X
(X i − X ) ⊗o (X i − X ) .
N i =1

(3.1)

For the estimation of the operator C 1 we first note that it is sufficient to estimate C 1 up
to a multiplicative constant, due to the self-normalizing form of the second term of the
minimum distance D(C 1 ). Let ∆ be any fixed element of S 2 (H2 ), note that under the null
hypothesis of separability H0 : C = C 1 ⊗C 2 we have T2 (C , ∆) = 〈C 2 , ∆〉S 2 (H2 )C 1 , and for this
choice the minimum distance in (2.15) is given by
D 0 := D(T2 (C , ∆)) = C 22 −

T1 (C , T2 (C , ∆))2
T2 (C , ∆)22

.

(3.2)

Therefore, we propose to estimate (a multiple of) the operator C 1 by
Cb1N = T2 (CbN , ∆),
10

(3.3)

and obtain the test statistic
D̂ N = CbN 22 −

T1 (CbN , T2 (CbN , ∆))2
.
T2 (CbN , ∆)2

(3.4)

2

3.1 Weak convergence
The main result of this section provides the asymptotic properties of the statistic D̂ N .
Theorem 3.1. If EkX k42 < ∞, then

p

¡
¢
b N − D 0 converges weakly to a zero mean normal
N D

distribution, where the centering term D 0 is defined in (3.2).
Proof. Observing (3.2) and (3.4) we write
Ã
!
2
2
b
b
p ¡
¢ p
T
(
C
,
T
(
C
,
∆))
T
(C
,
T
(C
,
∆))
1
2
b N − D 0 = N CbN 22 − 1 N 2 N
N D
− C 22 +
2
b
T
(C
, ∆)22
T2 (C N , ∆)2
2
(3.5)
b
and note that D N and D 0 are functions of the random variables
¡
¢T
G N = CbN 22 , T1 (CbN , T2 (CbN , ∆))2 , T2 (CbN , ∆)22 ,
¡
¢T
G = C 22 , T1 (C , T2 (C , ∆))2 , T2 (C , ∆)22 ,
respectively. Therefore, we first investigate the weak convergence of the vector

p

N (G N −

G). For this purpose we note that for A and B in S 2 (H ), the identity
A22 − B 22 = A − B 22 + 2〈A − B, B 〉H S
holds and introduce the decomposition
p

N (G N −G) = I N + I I N ,

where the random variables I N and I I N are defined by
p ¡
¢T
I N = N CbN −C 22 , T1 (CbN , T2 (CbN , ∆)) − T1 (C , T2 (C , ∆))22 , T2 (CbN , ∆) − T2 (C , ∆)22 ,
p ¡­
®
­
®
I I N =2 N CbN −C ,C H S , T1 (CbN , T2 (CbN , ∆)) − T1 (C , T2 (C , ∆)), T1 (C , T2 (C , ∆)) H S ,
­
® ¢T
T2 (CbN , ∆) − T2 (C , ∆), T2 (C , ∆) H S .
Using the linearity of T1 andT2 we further write
T1 (CbN , T2 (CbN , ∆)) − T1 (C , T2 (C , ∆)) =T1 (CbN , T2 (CbN , ∆)) − T1 (C , T2 (CbN , ∆))
+T1 (C , T2 (CbN , ∆)) − T1 (C , T2 (C , ∆))
=T1 (CbN −C , T2 (CbN , ∆)) + T1 (C , T2 (CbN −C , ∆)).

11

and obtain the representation



p
2


b
 N (C N −C )
 

2
2 
p
¢
¡
¢
¡p
1 



I N = p  T1 N (CbN −C ), T2 (CbN , ∆) + T1 C , T2 ( N (CbN −C ), ∆)  

2 
N
 ¡p
2
¢


T2 N (CbN −C ), ∆ 
2
³p
´
1
=: p F 1 N (CbN −C ), CbN ,
N


­p
®
N (CbN −C ),C H S
p
¢
¡
¢
®
 ­ ¡p

I I N =2  T1 N (CbN −C ), T2 (CbN , ∆) + T1 C , T2 ( N (CbN −C ), ∆) , T1 (C , T2 (C , ∆)) H S 
­ ¡p
¢
®
T2 N (CbN −C ), ∆ , T2 (C , ∆) H S
³p
´
=:F 2 N (CbN −C ), CbN ,

where the last equations define the functions F 1 and F 2 in an obvious manner. Note that
F := (F 1 , F 2 ) : S 2 (H ) × S 2 (H ) 7→ R6
is composition of continuous functions and hence continuous. By Proposition 5 in
p
Dauxois et al. (1982) the random variable N (CbN − C ) converges in distribution to a
centered Gaussian random element G with variance
p
Γ := lim Var( N CbN )

(3.6)

N →∞

in S 2 (H ) with respect to Hilbert-Schmidt topology. Therefore, using continuous mapping arguments we have
F

³p

´
d
N (CbN −C ), CbN → F (G ,C ) ,

and consequently
p

³p
´
³p
´
1
d
N (G N −G) = p F 1 N (CbN −C ), CbN + F 2 N (CbN −C ), CbN → F 2 (G ,C ).
N

We write


〈G ,C 〉H S





F 2 (G ,C ) = 2  〈T1 (G , T2 (C , ∆)) + T1 (C , T2 (G , ∆)) , T1 (C , T2 (C , ∆))〉H S  ,
〈T2 (G , ∆) , T2 (C , ∆)〉H S
which can be further simplified as


〈G ,C 〉H S





˜ 2 (C , ∆)〉H S + 〈C , T1 (C , T2 (C , ∆))⊗T
˜ 2 (G , ∆)〉H S  .
F 2 (G ,C ) = 2  〈G , T1 (C , T2 (C , ∆))⊗T
˜ HS
〈G , T2 (C , ∆)⊗∆〉
(3.7)
12

By Proposition 2.1 T2 (G , ∆) is a Gaussian Process in S 2 (H2 ). This fact along with Lemma
A.3 in Appendix A imply that F 2 (G ,C ) is a normal distributed random vector with mean
zero and covariance matrix, say Σ. By (3.5),
p

¡
¢ p
b N − D 0 = N ( f (G N ) − f (G)) ,
N D

where the function f : R3 7→ R is defined by f (x, y, z) = x − y/z. Therefore, using the delta
method and the fact that
¡
¢
¡
¢
P T2 (CbN , ∆)22 > 0 → P T2 (C , ∆)22 > 0 = 1
as C 2 6= 0, we finally obtain
p

¡
¢ d
¡
¢
b N − D0 →
N D
N 0, (∇ f (G))T Σ(∇ f (G))

(3.8)

as n → ∞ where ∇ f (x, y, z) = (1, −1/z, −y/z 2 )T denotes the gradient of the function f .

Remark 3.1. A sufficient condition for Theorem 3.1 to hold is EkX k42 < ∞. This is weaker
than Condition 2.1 in Aston et al. (2017), as indicated in Remark 2.2 (1) of their paper.
These authors used this assumption to prove weak convergence under the trace-norm
topology, which is required to establish Theorem 2.3 in Aston et al. (2017). In contrast the
proof of Theorem 3.1 here only requires weak convergence under the Hilbert-Schmidt
topology, which defines a weaker topology.
The following result specializes Theorem 3.1 to the null hypothesis of separability,
where a simple expression for the asymptotic variance is available.
Corollary 3.1. If EkX k42 < ∞ and the null hypothesis of separability holds, we have
p

d

b N → N (0, ν2 )
ND

as N → ∞, where the asymptotic variance is given by
ν2 =

16C 42
T2 (C , ∆)42

D

˜ T2 (C , ∆)⊗∆
˜
ΓT2 (C , ∆)⊗∆,

E
HS

(3.9)

and the operator Γ is defined in (3.6).
Proof. Under the null hypothesis of separability, we have
T2 (C , ∆) =〈C 2 , ∆〉S 2 (H2 )C 1 ,
T1 (C , T2 (C , ∆)) =〈C 2 , ∆〉S 2 (H2 ) T1 (C ,C 1 ) = 〈C 2 , ∆〉S 2 (H2 ) C 1 22C 2
13

(3.10)

and as a consequence we obtain
D 0 = C 22 −

T1 (C , T2 (C , ∆))2
T2 (C , ∆)22

= C 22 − C 1 22 C 2 22 = 0.

Therefore, the convergence to a zero mean normal distribution follows from Theorem
3.1. To calculate the asymptotic variance, note that under the null hypothesis of separability the second coordinate of F 2 (G ,C ) in (3.7) can be simplified using Proposition 2.1
and (3.10) as follows,
˜ 1 (C , T2 (C , ∆))〉H S
〈T1 (G , T2 (C , ∆)), T1 (C , T2 (C , ∆))〉H S =〈G , T2 (C , ∆)⊗T
˜ 2 〉H S
=〈C 2 , ∆〉2S 2 (H2 ) C 1 22 〈G ,C 1 ⊗C
=T2 (C , ∆)22 〈G ,C 〉H S
D
E
〈T1 (C , T2 (G , ∆)), T1 (C , T2 (C , ∆))〉H S = 〈C 1 , T2 (G , ∆)〉H S C 2 , 〈C 2 , ∆〉S 2 (H2 ) C 1 22C 2
˜ HS
=C 1 22 C 2 22 〈C 2 , ∆〉S 2 (H2 ) 〈G ,C 1 ⊗∆〉
˜ HS.
=C 22 〈G , T2 (C , ∆)⊗∆〉
The elements in the covariance matrix Σ = (Σi j ) of F 2 (G ,C ) are therefore given by
D
E
Σ11 =4 ΓC ,C

HS

E
˜
˜
ΓC ,C
ΓT2 (C , ∆)⊗∆, T2 (C , ∆)⊗∆
HS
D
E
2
2
˜
+8T2 (C , ∆)2 C 2 ΓC , T2 (C , ∆)⊗∆
HS
D
E
˜ T2 (C , ∆)⊗∆
˜
Σ33 =4 ΓT2 (C , ∆)⊗∆,
HS
D
E
D
E
2
˜
Σ12 =4T2 (C , ∆)2 ΓC ,C
+ 4C 22 ΓC , T2 (C , ∆)⊗∆
HS
D
E HS
˜
Σ13 =4 ΓC , T2 (C , ∆)⊗∆
HS
D
E
D
E
2
˜
˜ T2 (C , ∆)⊗∆
˜
Σ23 =4T2 (C , ∆)2 ΓC , T2 (C , ∆)⊗∆
+ 4C 22 ΓT2 (C , ∆)⊗∆,

Σ22 =4T2 (C , ∆)42

D

E

+ 4C 42
HS

D

HS

HS

.

Using this the asymptotic variance ν2 = (∇ f (G))T Σ(∇ f (G)) in the proof of Theorem 3.1
can be simplified by a straightforward but tedious calculation to
ν2 =

16C 42

D
E
˜
˜
ΓT
(C
,
∆)
⊗∆,
T
(C
,
∆)
⊗∆
2
2
4

T2 (C , ∆)2

HS

,

which proves our assertion.
Let u 1−α denote the (1 − α) quantile of the standard normal distribution and let ν̂2N
be a consistent estimate of the limiting variance (3.9). Then it follows from Theorem 3.1

14

and Corollary 3.1 that the test which rejects the null hypothesis of separability, whenever
ν̂
b N > u 1−α pN
D
(3.11)
N
is an asymptotic level α test. The following result shows that this test is also consistent.
Corollary 3.2. Assume that EkX k42 < ∞. If the null hypothesis of separability is not true,
p
p
b N → ∞ as N → ∞.
we have N D
Proof. By the third part of Theorem 2.1 we can write
D 0 = C 22 −

T1 (C , T2 (C , ∆))22
T2 (C , ∆)22

.

As in the proof of Theorem 2.1, the last quantity can be shown to be equal to C −
˜ 1 (C , T2 (C , ∆))22 . Note that T2 (C , ∆) ∈ S 2 (H1 ) and T1 (C , T2 (C , ∆)) ∈ S 2 (H2 ), thereT2 (C , ∆)⊗T
fore D 0 is positive if H0 is not true. Hence the result follows from Theorem 3.1.

3.2 Hilbert-Schmidt integral operators
In the remaining part of this section we concentrate on the case, where X is a random
element in H = L 2 (S × T ), where S and T are bounded intervals. In this particular scenario, we choose ∆ to be an integral operator generated by a kernel ψ(t , t 0 ). With this
choice, using the explicit formula for the operator T1 described in Proposition 2.3 the
minimum distance can be expressed in terms of the corresponding kernels, that is
Z Z Z Z
D 0 = D(T2 (C , ∆)) =
c 2 (s, t , s 0 , t 0 )d sd s 0 d t d t 0
(3.12)
T T S S
¤
R R £R R
0 0
0
0 2
dtdt0
T T
S S c(s, t , s t )c̃ 1 (s, s )d sd s
R R 2
.
−
0
0
S S c̃ 1 (s, s )d sd s
where c̃ 1 denote the kernel corresponding to the operator T2 (C , ∆), that is
Z Z
c̃ 1 (s, s 0 ) =
c(s, t , s 0 , t 0 )ψ(t , t 0 )d t d t 0 .
T

T

In this case the estimator CbN defined in (3.1) is induced by the kernel
ĉ N (s, t , s 0 , t 0 ) =

N
1 X
(X i (s, t ) − X (s, t ))(X i (s 0 , t 0 ) − X (s 0 , t 0 )) ,
N i =1

and the estimator Cb1N = T2 (CbN , ∆) defined in (3.3) is induced by the kernel
Z Z
0
ĉ 1N (s, s ) =
ĉ N (s, t , s 0 , t 0 )ψ(t , t 0 )d t d t 0 .
T

T

15

The estimator D̂ N of D 0 is calculated by plugging in ĉ N and ĉ 1N to the expression in
p
(3.12). By Corollary 3.1, under the null hypothesis of separability, the statistic N D̂ N
is asymptotically centered normal distributed with variance given by (3.9). In order to
define an estimator of this variance we introduce the following technical lemma which
is proved in Appendix A.
Lemma 3.1. Assume E||X ||42 < ∞ and that H = L 2 (S × T, R). Then the operator Γ ∈ S 2 (H )
defined by (3.6) is induced by the kernel
γ(s 1 , t 1 , s 2 , t 2 , s 3 , t 3 , s 4 , t 4 ) =k(s 1 , t 1 , s 2 , t 2 , s 3 , t 3 , s 4 , t 4 ) + c(s 1 , t 1 , s 3 , t 3 )c(s 2 , t 2 , s 4 , t 4 )
+ c(s 1 , t 1 , s 4 , t 4 )c(s 2 , t 2 , s 3 , t 3 ) − 2c(s 1 , t 1 , s 2 , t 2 )c(s 3 , t 3 , s 4 , t 4 ),
(3.13)
where the kernel k is given by
¡
k(s 1 , t 1 , s 2 , t 2 , s 3 , t 3 , s 4 , t 4 ) =E (X (s 1 , t 1 ) − µ(s 1 , t 1 ))(X (s 2 , t 2 ) − µ(s 2 , t 2 ))
¢
(X (s 3 , t 3 ) − µ(s 3 , t 3 ))(X (s 4 , t 4 ) − µ(s 4 , t 4 )) .

(3.14)

By Proposition 5 of Dauxois et al. (1982) if EkX k82 < ∞ the 4-th order kernel k defined
in (3.14) can be consistently estimated by
k̂ N (s 1 , t 1 , s 2 , t 2 , s 3 , t 3 , s 4 , t 4 ) =

N ³
1 X
(X i (s 1 , t 1 ) − X (s 1 , t 1 ))(X i (s 2 , t 2 ) − X (s 2 , t 2 ))
N i =1

´
(X i (s 3 , t 3 ) − X (s 3 , t 3 ))(X i (s 4 , t 4 ) − X (s 4 , t 4 )) ,
£
¤
i.e., k̂ N converges in probability to k as an element of L 2 (S × T )4 , R . A consistent estimator of the limiting covariance kernel γ is obtained by replacing k and c by k̂ N and ĉ N ,
respectively in (3.13), i.e.,
γ̂N (s 1 , t 1 , s 2 , t 2 , s 3 , t 3 , s 4 , t 4 ) = k̂ N (s 1 , t 1 , s 2 , t 2 , s 3 , t 3 , s 4 , t 4 ) + ĉ N (s 1 , t 1 , s 3 , t 3 )ĉ N (s 2 , t 2 , s 4 , t 4 )
+ ĉ N (s 1 , t 1 , s 4 , t 4 )ĉ N (s 2 , t 2 , s 3 , t 3 ) − 2ĉ N (s 1 , t 1 , s 2 , t 2 )ĉ N (s 3 , t 3 , s 4 , t 4 ).
(3.15)
Indeed a consistent estimator of the variance in (3.9) is given by
ν̂2 =

E
16CbN 42 D
b N T2 (CbN , ∆)⊗∆,
˜ T2 (CbN , ∆)⊗∆
˜
Γ
.
HS
T2 (CbN , ∆)4
2

b N is the operator induced by γ̂N by right integration.
where Γ

16

(3.16)

Table 1: Empirical rejection probabilities of different tests for the hypothesis of separability (level 5%). M1: the test (3.11) proposed in this paper, the different kernels are indicated
in brackets; M2: the test of Aston et al. (2017); M3: the test of Constantinou et al. (2017).
The data are generated from model (4.1), where the case β = 0 corresponds to the null
hypothesis of separability.
β

N

M1 (ψ1 )

M1 (ψ2 )

M1 (ψ3 )

M1 (ψ4 )

M2

M3

0

100

3.1 %

4.7%

3.8%

4.1%

0.4%

4.3 %

0

200

3.7%

4.2%

4.4%

4.9%

4.3%

5.3%

0

500

4.5%

5.1%

3.1%

3.2%

4.9 %

3.8%

0.3

100

33.2 %

35.1%

32.6%

33.8%

45.8%

11.3 %

0.3

200

48.3 %

46.8%

57.7%

51.9%

60.1%

22.1%

0.3

500

57.6%

63.7%

60.3%

61.8%

64.4%

35.6 %

0.5

100

59.6 %

67.5%

60.8%

63.4%

56.9 %

26.1 %

0.5

200

78.9%

80.0%

81.3%

84.6%

85.7%

41.2%

0.5

500

90 %

92.3%

90.2%

91.5%

90.1%

53.6%

0.7

100

73.8 %

75.6%

73.6%

74.3%

72.8 %

47.9%

0.7

200

87.9 %

90.2%

89.8%

88.6 %

89.4 %

58.5 %

0.7

500

98.3 %

99.0%

98.9%

99.2%

99.3%

67.9%

1

100

92.6%

95.3%

96.1%

95.4%

95.8%

61.7%

1

200

98.4%

99.9%

98.7%

100%

98.1 %

83.2%

1

500

100%

100%

100%

100%

100%

90.1%

Remark 3.2. Although the proposed estimator is based on the norm of the complete covariance kernel c, numerically we do not need to store the complete covariance kernel.
For example, we obtain for the first term of the statistic D̂ N the representation
1
CbN 22 = 2
N

Z Z Z Z ÃX
N
T

S T

S i =1

!2
(X i (s, t ) − X (s, t ))(X i (s 0 , t 0 ) − X (s 0 , t 0 ))

d sd t d s 0 d t 0

¸2
N X
N ·Z Z
1 X
(X i (s, t ) − X (s, t ))(X j (s, t ) − X (s, t ))d sd t
= 2
N i =1 j =1 T S
All other terms of the estimator in (3.4) and the variance estimator (3.16) can be represented similarly using simple matrix operations on the data matrix without storing the
full or marginal covariance kernels.

17

Table 2: Empirical rejection probabilities of different tests for the hypothesis of separability (level 5%). M1: the test (3.11) proposed in this paper, the different kernels are indicated
in brackets; M2: the test of Aston et al. (2017); M3: the test of Constantinou et al. (2017).
The data are generated from model (4.2), where the case c 0 = 1 corresponds to the null
hypothesis of separability.

c0

N

M1 (ψ1 )

M1 (ψ2 )

M1 (ψ3 )

M1 (ψ4 )

M2

M3

1

100

4.3 %

4.1%

4.0%

3.8%

5.5%

8.1 %

1

200

3.8%

3.5%

3.9%

4.2%

4.2 %

6.3%

1

500

3.2%

3.9%

5.2%

4.3%

3.1%

7.2%

3

100

38.8 %

41.3%

40.8%

41.9 %

39.1%

9.3 %

3

200

49.4 %

50.4%

52.3%

55.8%

59.7 %

13.4 %

3

500

66.7%

66.2%

65.9%

67.4%

66.8%

22.8%

5

100

60.8 %

63.1%

62.4%

61.1%

61.7 %

38.2 %

5

200

81.6 %

82.9%

83.4%

80.8%

82.7%

52.9%

5

500

94.3%

94.8%

93.7%

92.9%

94.9 %

73.5%

7

100

78.1 %

78.8%

77.7 %

73.6%

75.4%

59.4%

7

200

87.4%

89.9%

91.3%

90.2%

91.4 %

83.3%

7

500

98.8 %

99.7%

99.9%

99.4%

99.9%

91.8%

10

100

93.5%

95.1%

94.3%

93.7%

93.9%

75.2%

10

200

99.2 %

100%

99.8%

100%

100 %

91%

10

500

100 %

100%

100%

100%

100%

96.3%

4 Finite Sample Properties
In this section we study the finite sample properties of the test for the hypothesis of separability defined by (3.11) by means of a small simulation study. We also compare the
new test with the tests proposed by Aston et al. (2017) and Constantinou et al. (2017) and
illustrate potential applications with a data example. For this purpose we have implemented the new method proposed in this paper as well as the bootstrap test described
in Aston et al. (2017) and the weighted χ2 test based on the test statistic TbF as described
in Theorem 3 of Constantinou et al. (2017). The new test (3.11) depends on the choice
of the operator ∆ and we demonstrate in a simulation study that the test is not very sensitive with respect to this choice. For this purpose we consider four integral operators
with the following kernels
ψ1 (t , t 0 ) ≡ 1 , ψ2 (t , t 0 ) = |t − t 0 |,
ψ3 (t , t 0 ) = min(t , t 0 ) , ψ4 (t , t 0 ) = exp(−π(t 2 + t 02 )).
18

The test due to Aston et al. (2017) requires specification of eigen-subspace I , which was
taken to be {(1, 1)} as in their simulation studies and p-values were obtained by empirical
non-studentized bootstrap. We used the R package "covsep" (Tavakoli (2016)) to implement their method. For the tests proposed by Constantinou et al. (2017) we choose the
procedure based on the statistic TbF as in a simulation study it turned out to be the most
powerful procedure among the four methods proposed in this paper. The test requires
the specifications of the number of spatial and temporal principle components, which
were both taken to be 2.

4.1 Simulation Studies
The data were generated from a zero mean Gaussian process with two different covariance kernels. The first is the spatio-temporal covariance kernel
µ
¶
cks − s 0 k2γ
σ2
exp −
c(s, t , s , t ) =
,
(a|t − t 0 |2α + 1)τ
(a|t − t 0 |2α + 1)βγ
0

0

(4.1)

introduced by Gneiting (2002). In this covariance function, a and c are nonnegative
scaling parameters of time and space, respectively; α and γ are smoothness parameters
which take values in the interval (0, 1]; β is the separability parameter which varies in
the interval [0, 1]; σ2 > 0 is the point-wise variance; and τ ≥ βd /2, where d is the spatial
dimension. If β = 0, the covariance is separable and the space-time interaction becomes
stronger with increasing values of β. We fix γ = 1, α = 1/2, σ2 = 1, a = 1, c = 1 and τ = 1 in
the following discussion and choose different values for the parameter β specifying the
level of separability.
As a second example we consider covariance structure
c(s, t , s 0 , t 0 ) =

σ2 c 0d /2
³
´d /2
(a 02 (t − t 0 )2 + 1)1/2 a 02 (t − t 0 )2 + c 0

(4.2)

h
³ a 2 (t − t 0 )2 + 1 ´1/2 i
,
× exp − b 0 ks − s 0 k 20
a 0 (t − t 0 )2 + c 0
which was introduced by Cressie and Huang (1999). Here a 0 and b 0 are non negative
scaling parameters of time and space respectively; c 0 > 0 is the separability parameter,
σ2 > 0 is the point-wise variance and d is the spatial dimension. The covariance is separable if c 0 = 1. For the simulation study we take a 0 = 2, b 0 = 1, σ2 = 1, d = 2 and consider
different values of the parameter c 0 .
We generate data at 100 equally spaced time points in [0, 1] and 11 space points on
the grid [0, 1] × [0, 1]. The integrals are approximated by average of the function value
at grid points. The nominal significance level is taken to be 5% and empirical rejection
region are computed based on 1000 Monte-Carlo replications. The results are displayed
19

in Table 1 and 2, where we show the results for the test (3.11) proposed in this paper
(M1), the test based on projections of subspaces suggested in Aston et al. (2017) (M2)
and the test and the weighted χ2 test of Constantinou et al. (2017) (M3).
All procedures yield rather similar results under the null hypothesis, and in general
the nominal level is well approximated. Only for model (1) and sample size N = 100 the
bootstrap test of Aston et al. (2017) is very conservative, while nominal level of the test
of Constantinou et al. (2017) is slightly too large in model (2). Under the alternative this
test yields less power than the test of Aston et al. (2017) and the test (3.11) proposed in
this paper. On the other hand we do not observe substantial differences between the
test proposed in this paper and the bootstrap test suggested by Aston et al. (2017). In
particular the test (3.11) seems to be very robust with respect to different choices for the
kernel ψ. The constant kernel ψ1 yields slightly less power, but computationally it is a
little faster than the other choices.

Figure 1: Average wind speed curves for January 1961

4.2 Application to Real Data
4.2.1 Irish Wind Data
We finally apply the new method to the Irish wind data of Haslett and Raftery (1989).
The data is available at http://lib.stat.cmu.edu/datasets/wind.data and consists of daily average wind-speed at 12 meteorological stations in Ireland during the period 1961–1978. The locations of the stations can be found in Table 3 of Gneiting (2002).
We treat the data as functional observations X n (s, t ), for each month n, with N = 216. A
plot of the data for January 1961 is presented in Figure 1.

20

Figure 2: Minimum Temperature curves for January 1971

(a) Coastal Region

(b) Inland Region

Figure 3: Minimum Temperature curves (Coast and Inland) for January 1971
The data has a seasonal component. We deseasonalize the data by subtracting the
monthly mean from each curve as in Constantinou et al. (2017). The new test (3.11) is
applied to the deseasonalized data with four choices of kernels as mentioned in Section
4.1. The p-values are presented in Table 3. This suggests departure from separability,
which agrees with the findings of both Gneiting (2002) and Constantinou et al. (2017).
Table 3: P-values of the test (3.11) (with different kernels) for the Irish Wind data
Kernels

ψ1

ψ2

ψ3

ψ4

p-value

0.012

0.008

0.003

0.005

4.2.2 Spanish Temperature Data
We also analyze the Spanish Temperature Data used in Cuesta-Albertos and FebreroBande (2010). The data consists of daily minimum temperature of 15 meteorological
sites of northern part of Spain during the period 1971–200. Each of these stations are
located at different provinces covering 6 coastal and 9 inland areas of Spain. As before
21

we treat them as functional data X n (s, t ), for each month n with N = 360. The temperature of January 1971 is presented in Figure 2. As indicated by Cuesta-Albertos and
Febrero-Bande (2010), the data from inland and coastal region behave quite differently,
the coastal area being more stable. This is also evident from Figure 3, showing temperature curve for these two regions separately.
The data is deseasonalized before analysis by subtracting the monthly means from
each curve. The new test (3.11) is then applied to both the whole data and to the data
restricted to the two regions with four choices of kernel as before. All the tests yield
p-values very close to 0 (< 0.001). This indicates a strong non-separability in the data.
Note that Cuesta-Albertos and Febrero-Bande (2010) also found strong space-time interaction, even after considering the effect of the regions. Therefore the violation of
separability assumption is expected in this scenario.
Acknowledgments. This work has been supported in part by the Collaborative Research
Center “Statistical modelling of nonlinear dynamic processes” (SFB 823, Teilprojekt A1,
C1) of the German Research Foundation (DFG). The authors would like to thank Martina
Stein, who typed parts of this paper with considerable technical expertise. The authors
are also grateful to Juan Cuesta Albertos and to Shahin Tavakoli for very helpful discussions.

References
Aston, J. A. D., Pigoli, D., and Tavakoli, S. (2017). Tests for separability in nonparametric
covariance operators of random surfaces. The Annals of Statistics, 45(4):1431–1461.
Bar-Hen, A., Bel, L., and Cheddadi, R. (2008). Spatio-temporal functional regression on
paleoecological data. Functional and Operatorial Statistics, pages 54–56.
Constantinou, P., Kokoszka, P., and Reimherr, M. (2017). Testing separability of spacetime functional processes. Biometrika, 104(2):425–437.
Cressie, N. and Huang, H.-C. (1999).
stationary covariance functions.

Classes of nonseparable, spatio-temporal

Journal of the American Statistical Association,

94(448):1330–1339.
Cuesta-Albertos, J. and Febrero-Bande, M. (2010). A simple multiway anova for functional data. Test, 19(3):537–557.
Dauxois, J., Pousse, A., and Romain, Y. (1982). Asymptotic theory for the principal component analysis of a vector random function: some applications to statistical inference. Journal of multivariate analysis, 12(1):136–154.
22

Dunford, N. and Schwartz, J. T. (1988). Linear Operators. Vols. I, II, III. Wiley Classics
Library. John Wiley & Sons, Inc., New York.
Fuentes, M. (2006). Testing for separability of spatial-temporal covariance functions.
Journal of Statistical Planning and Inference, 136(2):447–466.
Genton, M. G. (2007). Separable approximations of space-time covariance matrices.
Environmetrics, 18:681–695.
Gneiting, T. (2002). Nonseparable, stationary covariance functions for space–time data.
Journal of the American Statistical Association, 97(458):590–600.
Gneiting, T., Genton, M. G., and Guttorp, P. (2007). Geostatistical space-time models,
stationarity, separability and full symmetry. In Finkenstadt, B., Held, L., and Isham,
V., editors, Statistical Methods for Spatio-Temporal Systems, pages 151–175. Chapman
and Hall/CRC, Boca Raton.
Gohberg, I., Goldberg, S., and Kaashoek, M. A. (1990). Classes of Linear Operators. Vol. I.
Operator Theory: Advances and Applications 49. Birkhäuser, Basel.
Haslett, J. and Raftery, A. E. (1989). Space-time modelling with long-memory dependence: Assessing ireland’s wind power resource. Applied Statistics, pages 1–50.
Hsing, T. and Eubank, R. (2015). Theoretical Foundations of Functional Data Analysis,
with an Introduction to linear Operators. Wiley, New York.
Huizenga, H. M., de Munck, J. C., Waldorf, L. J., and Grasman, R. P. P. P. (2002). Spatiotemporal EEG/MEG source analysis based on a parametric noise covariance model. IEEE
Transactions on Biomedical Engineering, 49(6):533–539.
Lifshits, M. (2012). Lectures on gaussian processes. In Lectures on Gaussian Processes,
pages 1–117. Springer.
Lu, N. and Zimmerman, D. L. (2005). The likelihood ratio test for a separable covariance
matrix. Statistics & Probability Letters, 73:449–457.
Mitchell, M. W., Genton, M. G., and Gumpertz, M. L. (2005). Testing for separability of
space-time covariances. Environmetrics, 16:819–831.
Mitchell, M. W., Genton, M. G., and Gumpertz, M. L. (2006). A likelihood ratio test for
separability of covariances. Journal of Multivariate Analysis, 97(5):1025–1043.
Rabiner, L. R. and Schafer, R. W. (1978). Digital Processing of Speech Signals 100. PrenticeHall, Englewood Cliffs.
23

Rougier, J. (2017).

A representation theorem for stochastic processes with sepa-

rable covariance functions, and its implications for emulation.

arXiv preprint

arXiv:1702.05599.
Skup, M. (2010). Longitudinal fmri analysis: A review of methods. Statistics and its
interface, 3(2):235.
Tavakoli, S. (2016).

covsep: Tests for Determining if the Covariance Structure of 2-

Dimensional Data is Separable. R package version 1.0.0.
Weidmann, J. (1980). Linear operators in Hilbert spaces. Graduate texts in mathematics.
Springer-Verlag.
Worsley, K. J., Marrett, S., Neelin, P., Vandal, A. C., Friston, K. J., and Evans, A. C. (1996).
A unified statistical approach for determining significant signals in images of cerebral
activation. Human Brain Mapping, 4:58–73.

A Some technical results
A.1 Properties of Tensor Product Hilbert Spaces
In this section we present some auxiliary results which have been used in the proofs of
the main results. The following two results are generalizations of Lemma 1.6 and Lemma
1.7 in the online supplement of (Aston et al., 2017).
Lemma A.1. The set
o
nX
n
˜ i : A i ∈ S 2 (H1 ), B i ∈ S 2 (H2 ), are finite rank operators , n ∈ N
A i ⊗B
D0 :=

(4.3)

i =1

is dense in S 2 (H1 ⊗ H2 ).
P
Proof. Let T ∈ S 2 (H1 ⊗ H2 ), then T = j ≥1 λ j u j ⊗o v j where (λ j ) j ≥1 are the singular valP
ues of T , with j λ2j < ∞; and {u j } j ≥1 and {v j } j ≥1 are orthonormal subsets of H1 ⊗ H2 .
P
The vectors u j and v j ’s can be further decomposed as u j = l ≥1 s j ,l e l(1) ⊗ e l(2) and v j =
P
(1)
(2)
l ≥1 t j ,l f l ⊗ f l , where the series converge in the norm of H1 ⊗ H2 . Let
uM
j =

M
X
l =1

s j ,l e l(1) ⊗ e l(2) and v M
j =

M
X
l =1

Given 1 > ² > 0, choose N large enough such that T −
we write
 



 

N
N
X
X



M
M
λ j u j ⊗o v j 
λ j u j ⊗o v j  ≤ T −
T −




j =1
j =1

2

2

24

t j ,l f l(1) ⊗ f l(2) .

PN

j =1 λ j u j ⊗o v j 2

≤ ²/2. With this



X
h
i
N

M
M
M 
+  λ j (u j − u j ) ⊗o v j + u j ⊗o (v j − v j ) 
 j =1


2

≤ ²/2 +

N
X
j =1

"




M
M 
λ j (u j − u M
)
⊗
v
+
u
⊗
(v
−
v
)
o j
o
j
j
j
j 

2

N
X

≤ ²/2 +

j =1

#1/2 "
λ2j
"

≤ ²/2 + T 2


N 
X

M
M 
)
⊗
v
+
u
⊗
(v
−
v
)
(u j − u M
o j
o
j
j
j
j 

#1/2

2

j =1

°° ° °
°°
°i
N h°
X
°
°° ° ° M°°
M °
)
v
+
−
v
)
°(u j − u M
°
°u
°
°(v
j
j
j
j
j °

#1/2
.

j =1

Choose M ≥ 1 such that
(
ku j − u M
j k ≤ min

²2

)
and kv j − v M
j k≤

,1
2

12N T 2

²2
12N T 22

.

k, with this choice we have
k ≤ ku j k + ku j − u M
As ku j k = 1 and ku M
j
j




N
X

M
M
λ j u j ⊗o v j  ≤ ²
T −


j =1
2

The proof is now complete by noting that

PN

M
j =1 λ j u j

⊗o v M
is a finite rank operator.
j

Lemma A.2. Let C ∈ S 2 (L 2 (S × T, R)), where T and S are compact subsets of Rp and Rq
respectively, be an integral operator with symmetric continuous kernel c. For any ² > 0,
PN
˜ n , where A n : L 2 (S, R) → L 2 (S, R), B n : L 2 (T, R) → L 2 (T, R)
there exists an C 0 = n=1
A n ⊗B
are finite rank operators with continuous kernels a n and b n respectively, such that
(a) C −C 0 2 ≤ ²,
(b) sups,s 0 ∈S,t ,t 0 ∈T |c(s, t , s 0 , t 0 )−c 0 (s, t , s 0 , t 0 )| ≤ ², where c 0 is the kernel of the operator C 0 .
Proof. By Mercer’s Theorem, there exists continuous orthonormal functions {u n }n≥1 ⊂
L 2 (S × T, R) and λn is the summable sequence of positive eigenvalues, such that
X
c(s, t , s 0 , t 0 ) =
λn u n (s, t )u n (s 0 , t 0 )
n≥1

where the convergence is absolute and uniform.
Let Un be the integral operator with kernel u n and C N :=

PN

n=1 λn Un

⊗o Un . Denote

the kernel of C N as c N . As u N can be approximated by sums of tensor products of conP
(1)
(2)
(1)
(2)
M
tinuous functions, let u N
:= lM=1 f n,l
⊗ f n,l
, where f n,l
∈ L 2 (S, R), f n,l
∈ L 2 (T, R) are continuous, and choose M such that
(
kUn −UnM k ≤ min
25

²2
12N κC 22

)
,κ ,

where κ = maxn=1,...,N kUn k∞ . Writing C N ,M :=

PN

M
M
n=1 λn Un ⊗2 Un , and denoting the ker-

nel by c N ,M , we have
N

kc − c

N ,M

N
X

k∞ ≤

£
¤
λn kUn −UnM k∞ kUn k∞ + kUnM k∞ kUn −UnM k∞

n=1

An application of the Cauchy-Schwartz inequality along with the choice of M gives an
upper bound to the last quantity to be ²/2. Similar calculations show that
C N −C N ,M 2 ≤ ²/2.
Finally, as C N ,M is indeed a finite sum of tensor products of finite rank integral operators,
we have the desired result.
We conclude this section with a simple result about Gaussian processes on a Hilbert
space. For this purpose recall that a random element G on a real separable Hilbert space
H is said to be Gaussian with mean µ ∈ H and covariance operator Γ : H 7→ H if for
all x ∈ H , the random variable 〈G , x〉 has a normal distribution with mean 〈µ, x〉 and
variance 〈Γx, x〉. (See Section 1.3 of Lifshits (2012) for more details)
Lemma A.3. Let H1 and H2 be two real separable Hilbert spaces and G be a Gaussian
˜ is a Gaussian process in
process on S 2 (H2 ). Then for all A ∈ S 2 (H1 ), the process A ⊗G
S 2 (H1 ⊗ H2 ).
˜ 〉 has a normal
Proof. We will show for any T ∈ S 2 (H1 ⊗ H2 ), the random variable 〈T, A ⊗G
distribution. By Lemma A.1 and continuity of inner-product it is enough to show the
PN
˜ n then
result for T ∈ D0 . Therefore let T = n=1
A n ⊗B
˜ 〉=
〈T, A ⊗G

N
X

〈A n , A〉S 2 (H1 ) 〈B n , G 〉S 2 (H2 )

n=1

which is sum of normal random variables and hence normal.

A.2 Proof of Lemma 3.1:
Proof. Without loss of generality we assume E(X ) = µ = 0. To obtain the kernel γ, we
need to calculate the limit of
N Cov (ĉ N (s 1 , t 1 , s 2 , t 2 ), ĉ N (s 3 , t 3 , s 4 , t 4 )) .
To this end we write
E(ĉ N (s 1 , t 1 , s 2 , t 2 )ĉ N (s 3 , t 3 , s 4 , t 4 )) =

1
E((H X (s 1 , t 1 ))T (H X (s 2 , t 2 ))(H X (s 3 , t 3 ))T (H X (s 4 , t 4 )))
N2
26

with X (s, t ) = (X 1 (s, t ), X 2 (s, t ), · · · , X N (s, t ))T and H = I − N1 J where I is the N ×N identity
matrix and J is the N × N matrix with (J )i j = 1 for all i , j . Using the fact that H T H = H ,
we have
1
E(X (s 1 , t 1 )T H X (s 2 , t 2 )X (s 3 , t 3 )T H X (s 4 , t 4 ))
N2
´
X
1 ³X
= 2E
h i j X i (s 1 , t 1 )X j (s 2 , t 2 ) h kl X k (s 3 , t 3 )X k (s 4 , t 4 )
N
i,j
k,l
X
¡
¢
1
= 2 h i i h j j E X i (s 1 , t 1 )X i (s 2 , t 2 )X j (s 3 , t 3 )X j (s 4 , t 4 )
N i,j
¡
¢
1 X
+ 2 h i2j E X i (s 1 , t 1 )X i (s 3 , t 3 )X j (s 2 , t 2 )X j (s 4 , t 4 )
N i,j
¡
¢
1 X
+ 2 h i2j E X i (s 1 , t 1 )X i (s 4 , t 4 )X j (s 2 , t 2 )X j (s 3 , t 3 )
N i,j
¡
¢
1 X
+ 2 h i2i E X i (s 1 , t 1 )X i (s 2 , t 2 )X i (s 3 , t 3 )X i (s 4 , t 4 ) .
N i

E(ĉ N (s 1 , t 1 , s 2 , t 2 )ĉ N (s 3 , t 3 , s 4 , t 4 )) =

The last equality is due to the symmetry of H , the independence of the X i and the fact
that E(X i (s, t )) = 0. Further simplifying the last quantity we have
E(ĉ N (s 1 , t 1 , s 2 , t 2 )ĉ N (s 3 , t 3 , s 4 , t 4 )) =

(N − 1)2
c(s 1 , t 1 , s 2 , t 2 )c(s 3 , t 3 , s 4 , t 4 )
N2
(N − 1)
+
c(s 1 , t 1 , s 3 , t 3 )c(s 2 , t 2 , s 4 , t 4 )
N2
(N − 1)
c(s 1 , t 1 , s 4 , t 4 )c(s 2 , t 2 , s 3 , t 3 )
+
N2
µ
¶
(N − 1)2
1
+
k(s 1 , t 1 , s 2 , t 2 , s 3 , t 3 , s 4 , t 4 ) + O
.
N3
N2

Therefore we finally obtain
N Cov (ĉ N (s 1 , t 1 , s 2 , t 2 ), ĉ N (s 3 , t 3 , s 4 , t 4 )) =

1 − 2N
c(s 1 , t 1 , s 2 , t 2 )c(s 3 , t 3 , s 4 , t 4 )
N
(N − 1)
+
c(s 1 , t 1 , s 3 , t 3 )c(s 2 , t 2 , s 4 , t 4 )
N
(N − 1)
+
c(s 1 , t 1 , s 4 , t 4 )c(s 2 , t 2 , s 3 , t 3 )
N
µ ¶
1
(N − 1)2
+
k(s 1 , t 1 , s 2 , t 2 , s 3 , t 3 , s 4 , t 4 ) + O
,
2
N
N

and taking the limit of the last expression as N → ∞, gives the desired result.

27

