INFERENCE FOR THE CROSS-COVARIANCE OPERATOR OF
STATIONARY FUNCTIONAL TIME SERIES

arXiv:1704.00315v1 [math.ST] 2 Apr 2017

GREGORY RICE AND MARCO SHUM

Abstract. When considering two or more time series of functions or curves, for instance those derived from densely observed intraday stock price data of several companies, the empirical cross-covariance operator is of fundamental importance due to its role
in functional lagged regression and exploratory data analysis. Despite its relevance, statistical procedures for measuring the significance of such estimators are undeveloped.
We present methodology based on a functional central limit theorem for conducting
statistical inference for the cross-covariance operator estimated between two stationary, weakly dependent, functional time series. Specifically, we consider testing the null
hypothesis that two series possess a specified cross-covariance structure at a given lag.
Since this test assumes that the series are jointly stationary, we also develop a changepoint detection procedure to validate this assumption, which is of independent interest.
The most imposing technical hurdle in implementing the proposed tests involves estimating the spectrum of a high dimensional spectral density operator at frequency
zero. We propose a simple dimension reduction procedure based on functional PCA to
achieve this, which is shown to perform well in a small simulation study. We illustrate
the proposed methodology with an application to densely observed intraday price data
of stocks listed on the NYSE.

1. Introduction
Functional time series analysis (FTSA) has grown substantially in the last decade and
a half in order to provide methodology for functional data objects that are obtained
sequentially over time. Perhaps the most typical way such data arises is when dense
records of continuous time processes are segmented into collections of curves in some
natural way. For example, high frequency records of pollution levels may be segmented
to form daily pollution curves, or tick-by-tick asset price data may be used to construct
daily intraday price or return curves; see [4] and [28]. Other examples include sequentially observed curves that describe physical phenomena, as arise in functional magnetic
resonance imaging or in the observation of DNA minicircles; see [1] and [33]. We refer
the reader to [39] and [15] for overviews of the field of functional data analysis, and to
[8] and [22] for an introduction to FTSA.
Most of the developments in these two fields focus on analyzing functional data obtained
from a single source, e.g. intraday price curves derived from a single asset, or in comparing functional data from several independent populations. To give a few examples that
are related to this work, [31], [34], and [37] develop methods for performing inference for
the covariance operator of functional data. Using a self normalization approach, a two
sample test for the second order structure with functional time series data that allows
for some dependence across the populations is developed in [50].
However, in many situations of interest, functional data are obtained simultaneously from
two or more sources, e.g. intraday price curves derived from several assets. In such cases
1

2

GREGORY RICE AND MARCO SHUM

one often wishes to quantify the potentially complex dependence relationships between
such curves, and one way of achieving this is through the empirical cross-covariance
operator. Although the notion of the cross-covariance operator between random elements
in a Hilbert space was put forward over forty years ago in [5], statistical methodology for
estimating and performing further inference for the cross-covariance structure between
collections of curves seems quite new.
Measuring the cross-covariance between collections of curves has received some attention
in the context of multivariate longitudinal and functional data. Under the assumption
that the given data is a simple random sample of multivariate longitudinal data, [14],
[41], and [51] develop measures of cross-covariance between longitudinal data sources,
including measures based on canonical correlation analysis and principal component
analysis. Principal component analysis of multivariate functional data is also studied in
[11] and [12], and in [35] an analog of the covariance matrix for multivariate functional
data using Fréchet integration is defined. In each of these cases, the potential effect of
temporal dependence among the functional units is not considered.
In the context of bivariate functional time series, the cross-covariance operator and its
lagged versions are arguably of greater importance. Methods for lagged functional time
series regression, which have recently been put forward in [20] and [36], are naturally
based on the Fourier transform of the lagged cross-covariance operators. Moreover, the
initial exploratory analysis of any such series would typically begin by considering the
sequence of lagged cross-covariance operators to try to gain insight into the relationship
between the series.
Despite the apparent utility of the cross-covariance operator of functional time series,
statistical inference for it has not yet been considered, to the best of our knowledge.
In Chapter 4 of the seminal work of Bosq [8], a central limit theorem is given for the
covariance and auto-covariance operators of functional time series that may be represented as linear processes. A portmanteau-type test for independence of two functional
time series is developed in [26] based on the norms of cross-covariance operators at long
lags, but their test assumes under the null hypothesis that the individual series are independent, and is hence not suitable for quantifying the significance of estimates of the
cross-covariance between curves.
Additionally, when the data are obtained as bivariate functional time series, it is of
special interest to know if the cross-covariance structure changes during the observation
period. Being able to test for such a feature 1) helps validate the assumption of joint
stationarity needed to apply inferential procedures for the cross-covariance operator,
2) is of use for determining if the regression function changes in the functional lagged
regression problem, and 3) may be of independent interest since the presence and location
of such a change point may signify an important event. This problem has also not been
addressed, although several authors have considered analogous problems in the context
of finite dimensional time series; we refer to [13], [47], [48], and [2]. Change point analysis
for the mean of functional time series has been considered recently in [42] and [10].
In this paper, we consider two hypothesis testing problems: 1) tests for a specified
cross-covariance structure between two functional time series, e.g. that the two series
are uncorrelated at a given lag, and 2) change point tests for the covariance structure
within a given sample. Two varieties of test statistics are proposed in each of these
settings that are based on either the standard L2 distance or dimension reduction based
methods using a suitable principal component basis. These two classes of statistics

3

possess complimentary advantages, which we detail by means of a theoretical result.
The asymptotic properties of each test statistic are established assuming a general weak
dependence condition similar to the one introduced in [21], which includes nonlinear time
series, as well as the majority of functional time series models studied to date, under
mild regularity conditions.
This methodology is primarily motivated by a basic observation of Brillinger [9] that
inference for the covariance and/or cross-covariance of time series can be made by performing inference for the mean of a suitably constructed series. This idea was utilized
in [18] in order to conduct non-parametric inference for the cross-covariance matrix of
finite dimensional time series.
The rest of the paper is organized as follows. Section 2 contains the main assumptions
and notation of the paper, as well as asymptotic results for cross-covariance function
estimators under these assumptions. The results developed in Section 2 are utilized to
develop hypothesis tests in Section 3. Some details about the practical implementation of
the methods developed in Section 3 are provided in Section 4, which include methods for
overcoming the technical challenge of estimating the eigen-elements of a high-dimensional
spectral density operator at frequency zero that arise in the limiting distribution of
the test statistics. These testing and estimation procedures were studied by means of
Monte Carlo simulation, the results of which we present in Section 5. We illustrate our
methodology with an application to cumulative intraday return curves derived from the
1-minute resolution price of Microsoft and Exxon Mobile stock listed on the New York
stock exchange from the year 2001 in Section 6. All technical derivations and proofs are
provided in the appendices following these sections.
2. Asymptotic properties of cross-covariance function estimates
Before we proceed, we introduce a bit of notation. Let h·, ·id denote the standard inner
product on the space L2 [0, 1]d of real valued square integrable functions defined on [0, 1]d ,
1/2
and let k·kd = h·, ·id . We write f for the function f (t) when it does not cause confusion.
R
R1
We use the notation to denote 0 .
In this section, we suppose that {(Xi (t), Yi (s)), t, s ∈ [0, 1]}i∈Z is a jointly stationary
sequence of real valued stochastic processes whose sample paths are in L2 [0, 1] from
which we have observed a sample of length T , {(X1 (t), Y1 (s)), ..., (XT (t), YT (s))}. For
instance, (Xi (t), Yi (s)) could be used to denote the price of stock X and stock Y on day
i at intraday times t and s normalized to the unit interval. We let µX (t) = EX0 (t), and
µY (s) = EY0 (s), and define
CXY (t, s) = cov(X0 (t), Y0 (s)) = E[(X0 (t) − µX (t))(Y0 (s) − µY (s))],
to be the cross-covariance function (or kernel) between {Xi } and {Yi } at lag zero. CXY
defines the cross-covariance operator cXY : L2 [0, 1] → L2 [0, 1] via
Z
cXY (f )(t) =

CXY (t, s)f (s)ds.

This relationship implies that we may conduct inference for cXY by conducting inference
for the function CXY . Based on the sample, we may estimate CXY with

4

GREGORY RICE AND MARCO SHUM

bT xc
1 X
ĈXY (t, s, x) =
(Xi (t) − X̄(t))(Yi (s) − Ȳ (s)),
T i=1

which denotes the partial sample estimate of xCXY based on the first x proportion of
the sample, where
T
T
1X
1X
Xi (t), and Ȳ (t) =
Yi (t).
X̄(t) =
T i=1
T i=1

Under mild regularity conditions on the process {(Xi (t), Yi (s))}, which are implied by
the Assumption 2.1 below, ĈXY (t, s, 1) is a consistent estimator of CXY (t, s) in L2 [0, 1]2 .
The motivation for considering now the partial sample estimates of CXY is due to our
application to change point testing for the cross-covariance developed below. We note
that when interested in studying the cross-covariance of the series at some fixed lag that
is different from zero, say ` > 0, one can apply all of the below methods to the sample
{(X`+1 (t), Y1 (s)), ..., (XT (t), YT −` (s))} of length T − ` with only superficial changes to
the proofs and notation, and so we otherwise make no special mention of this case.
In order to derive the asymptotic properties of ĈXY , we make use of the following
assumption that imposes stationarity, weak dependence, and moment conditions on the
functional time series.
Assumption 2.1. (a) There exists a measurable function gXY : S ∞ → L2 [0, 1]×L2 [0, 1],
where S is a measurable space, and a sequence of independent and identically distributed
(iid) innovations {i , : i ∈ Z} taking values in S such that (Xi , Yi ) = gXY (i , i−1 , . . .).
(b) For all m ≥ 1, the m-dependent sequence (Xi,m , Yi,m ) = g(i , . . . , i−m+1 , ∗i−m,m , ∗i−m−1,m , . . .)
with ∗i,m being independent copies of i,0 , and {ε∗i,m : i ∈ Z} independent of {i : i ∈ Z},
satisfies for some p > 4,
EkXi − Xi,m kp

1/p

= O(m−α ), and EkYi − Yi,m kp

1/p

= O(m−β )

where α, β > 1.
Assumption 2.1(a) implies that {(Xi , Yi )} is a jointly stationary sequence of Bernoulli
shifts in L2 [0, 1] × L2 [0, 1] that is driven by an underlying iid innovation sequence. The
space of functional time series models contained within this class is quite large, including
the functional ARMA and GARCH processes; see [8] and [3]. Condition (b) defines a
type of Lp -m-approximability condition along the lines of [21], which can often be easily
verified when a time series model for the observations is given. The rate condition on the
decay of these coefficients, which is somewhat stronger than the main condition studied
in [21], is used in order to show that certain smoothed periodogram type spectral density
operator estimates defined below are consistent. These assumptions could be replaced
by mixing conditions and functional versions of cummulant summability conditions as
presented in [32] and [49], which are more comparable to the assumptions in [18].
Theorem 2.1. Under Assumption 2.1, there exists a sequence of Gaussian processes,
{ΓT (t, s, x), t, s, x ∈ [0, 1]}T ∈N , defined on the same probability space as {(Xi , Yi )}, that
satisfy

5

(2.1)


2
Z Z √ 
bT xc
CXY (t, s) − ΓT (t, s, x) dtds = oP (1),
sup
T ĈXY (t, s, x) −
T
0≤x≤1

where EΓT (t, s, x) = 0, and
cov(ΓT (t, s, x), ΓT (u, v, y)) = min(x, y)D(t, s, u, v),
where D(t, s, u, v) is the long run covariance function of the sequence {(Xi (t)−µX (t))(Yi (s)−
µY (s))}, namely
D(t, s, u, v) =

∞
X

cov((X0 (t) − µX (t))(Y0 (s) − µY (s)), (X` (u) − µX (u))(Y` (v) − µY (v)).

`=−∞

Theorem 2.1 provides a Skorokhod-Dudley-Wichura type characterization of an invariance principle for ĈXY that can be utilized to establish the asymptotic properties of
continuous functionals of ĈXY , and we consider several such statistics in Section 3 below in order to carry out hypothesis testing for CXY . The proof of Theorem 2.1 is
given in Appendix A. The function D describes the asymptotic covariance function of
√
T ĈXY (t, s, 1). D naturally defines a Hilbert-Schmidt integral operator, d : L2 [0, 1]2 →
2
L [0, 1]2 , given by
ZZ
d(f )(t, s) =
D(t, s, u, v)f (u, v)dudv,
which further defines an orthonormal basis of eigenfunctions ϕi in L2 [0, 1]2 , and a nonnegative sequence of eigenvalues λ1 ≥ λ2 ≥ ... ≥ 0 satisfying
(2.2)

d(ϕi )(t, s) = λi ϕi (t, s).

We define these quantities here as they appear in the limiting distributions and definitions
for the test statistics considered below.
3. Inference for the cross-covariance function
Theorem 2.1 points to some asymptotically validated methods to measure the significance
of estimates of CXY . For instance, we may wish to test based on the estimate ĈXY (·, ·, 1)
H0,1 : CXY = C0 versus HA,1 : CXY 6= C0 ,
where equality is understood in the L2 [0, 1]2 sense, and C0 is a given function of interest.
This null function might be determined from historical data, or taken to be zero in order
to test for zero cross-covariance between Xi and Yi at a given lag. Since the hypothesis
H0,1 is well posed only when the sequence {(Xi , Yi )} is, at least weakly, jointly stationary,
it is also of interest to determine whether or not this assumption is valid. We frame this
as a second hypothesis test of the time homogeneity of the cross-covariance against the
“at most one” change point in the cross-covariance alternative:
(1)

(2)

(T )

H0,2 : CXY = CXY = · · · = CXY , versus

6

GREGORY RICE AND MARCO SHUM
(1)

(k∗ )

(k∗ +1)

HA,2 : C1 = CXY = · · · = CXY 6= CXY

(T )

= · · · = CXY = C2 , for some k ∗ = bT θc, θ ∈ (0, 1),

(i)

where CXY (t, s) = cov(Xi (t), Yi (s)). We proceed by developing test statistics for each
of these hypotheses. In order to test H0,1 , we first consider a statistic based on the
normalized L2 distance of ĈXY to C0 :
FT = T kĈXY (·, ·, 1) − C0 k22 .
Corollary 3.1. Under Assumption 2.1 and H0,1 ,
D

FT −→

(3.1)

∞
X

λi Ni2 , as T → ∞,

i=1

where {λi , i ≥ 1} are defined in (2.2), and {Ni , i ≥ 1} are independent standard normal
random variables.
Corollary 3.1 shows that a test of asymptotic size α of H0,1 may be obtained by comparing
FT to the 1 − α quantile of the limit distribution given in (3.1), which depends on the
unknown eigenvalues of the operator d. In order to estimate these quantiles, one can
estimate a suitably large number of eigenvalues λi using an estimate of d, and then
continue by using these estimates to approximate the limiting distribution via Monte
Carlo simulation. The details of this implementation are discussed in Section 4 below,
including how to obtain consistent estimates of a finite number of the λ0i s.
The fact that the limiting distribution of FT is nonpivotal though encourages one to consider alternate test statistics based on projecting ĈXY into finite dimensional subspaces
of L2 [0, 1]2 . A natural choice of the finite dimensional space to choose is the one spanned
by the eigenbasis generated by d. In fact, it is a fairly straightforward calculation to
show that for any positive integer p, under Assumption 2.1, the inner products
√
h T (ĈXY (·, ·, 1) − CXY ), ϕi i2 , 1 ≤ i ≤ p
are asymptotically independent, and hence projecting into the directions of ϕi has the
effect of partitioning the centered estimator ĈXY into approximately mean zero and independent components. This is similar to the motivation provided for dynamical principal
component analysis of Brillinger [9], which has been studied in the context of functional
time series data in [19], and [32].
In this direction, let
FT,p =

√
p
X
h T (ĈXY (·, ·, 1) − C0 ), ϕ̂i i2
2

i=1

λ̂i

,

where ϕ̂i , and λ̂i , 1 ≤ i ≤ p are consistent estimates of ϕi and λi , 1 ≤ i ≤ p, i.e. they
satisfy
(3.2)

max kϕ̂i − ĉi ϕi k2 = oP (1), and max |λ̂i − λi | = oP (1),

1≤i≤p

1≤i≤p

where ĉi = sign(hϕi , ϕ̂i i). We discuss in Section 4 how to obtain consistent estimates of
λi and ϕi (t, s).

7

Corollary 3.2. Under Assumption 2.1 and H0,1 ,
D

FT,p −→ χ2 (p), as T → ∞,
where χ2 (p) denotes a chi-squared random variable with p degrees of freedom.
We now turn to the consistency and power properties of FT and FT,p . The following
result shows that both statistics diverge at rate T under HA,1 , so long as the difference
CXY − C0 is not orthogonal to the first p elements of the principal component basis
{ϕi , i ≥ 1}.
Theorem 3.1. Under Assumption 2.1, and HA,1 ,
p

FT,p P X hCXY − C0 , ϕi i22
FT P
.
−→ kCXY − C0 k2 , and
−→
T
T
λ
i
`=1
√
Moreover, if C0 = C0,T satisfies k T (CXY − C0 ) − CA k2 → 0 as T → ∞ for some
element CA of L2 [0, 1]2 , then
(3.3)

D

FT −→

kCA k22

+2

∞
X

1/2
λi hCA , ϕi iNi

+

∞
X

λi Ni2 ,

i=1

i=1

and
D

FT,p −→

p
X
hCA , ϕi i2
2

λi

i=1

+2

p
X
hCA , ϕi i2 Ni
i=1

1/2
λi

+

p
X

Ni2 .

i=1

This “local-alternative” result provides insight into the complimentary strengths and
weaknesses of both the norm based and dimension reduction based test statistics. Evidently if CXY − C0 is orthogonal to the first p principal components of d, then the
test statistic FT,p is not expected to have more than trivial power. Additionally in this
case, the norm based test statistic FT is expected to have improved power over other
alternatives in which CXY − C0 is of the same magnitude, since the second term on the
right hand side (3.3), which has mean zero, will have a smaller variance in this case.
Conversely, if CXY − C0 is contained in the subspace spanned by the first p principal
components of d, then FT,p is expected to be more powerful, since then this statistic
effectively defines a uniformly most powerful test of H0,1 , assuming the data is Gaussian,
in the p-dimensional subspace spanned by ϕ1 , ..., ϕp .
In order to test H0,2 versus HA,2 , we define analogously to FT and FT,p ,
ZT = T sup kĈXY (·, ·, x) − xĈXY (·, ·, 1)k22 .
0≤x≤1

and
ZT,p = T sup
0≤x≤1

p
X
hĈXY (·, ·, x) − xĈXY (·, ·, 1), ϕ̂i i2
2

i=1

λ̂i

.

ZT and ZT,p are each maximally selected CUSUM type statistics based on comparing the
partial sample estimates of CXY to the estimator from the whole sample. The following

8

GREGORY RICE AND MARCO SHUM

corollaries of Theorem 2.1 quantify the large sample behavior of these statistics under
H0,2 .
Corollary 3.3. Under Assumption 2.1 and H0,2 ,
∞
X
D
λi Bi2 (x), as T → ∞,
ZT −→ sup
0≤x≤1

i=1

where {λi , i ≥ 1} are defined in (2.2), and {Bi (x), i ≥ 1, x ∈ [0, 1]} are iid standard
Brownian bridges on [0, 1]. If in addition (3.2) holds, then
p
X

D

ZT,p −→ sup
0≤x≤1

Bi2 (x), as T → ∞,

i=1

It follows then that a test of H0,2 with asymptotic level α is obtained by rejecting if
ZT or ZT,p are larger than the 1 − α quantiles of their limiting distributions detailed
in Corollary 3.3. These limiting distributions may again be approximated using Monte
Carlo simulation. The statistics ZT and ZT,p diverge under HA,2 in conjunction with
some mild ergodicity assumptions, which we explain in Section A.3 in the appendix.
4. Implementation
Implementing the testing procedures outlined above requires the estimation of the eigenvalues, and, in case of the dimension reduction based test statistics FT,p and ZT,p , the
eigenfunctions of d. We first develop methodology for estimating d and its spectrum,
and then describe some numerical methods for carrying out this estimation.
4.1. Estimation of D, d, ϕi and λi . As d is simply a scalar multiple of the spectral
density operator at frequency zero of the stationary sequence {(Xi (t) − µX (t))(Yi (s) −
µY (s))} in L2 [0, 1]2 , as defined in [32], it may be naturally estimated with a smoothed
periodogram type estimator. Let
(4.1)

∞
X

 
`
D̂T (t, s, u, v) =
Wb
γ̂` (t, s, u, v),
h
`=−∞

where h is a bandwidth parameter satisfying,
h = h(T ) → ∞,

h

→ 0 as T → ∞,
T 1/2
and, with X̄j (t) = Xj (t) − X̄(t) and Ȳj (s) similarly defined,

T −`



1 X



X̄
(t)
Ȳ
(s)
−
Ĉ
(t,
s,
1)
X̄
(u)
Ȳ
(v)
−
Ĉ
(u,
v,
1)
,
`≥0
j
j
XY
j+`
j+`
XY


 T j=1
γ̂` (t, s, u, v) =

T



1 X 



X̄j (t)Ȳj (s) − ĈXY (t, s, 1) X̄j+` (u)Ȳj+` (v) − ĈXY (u, v, 1) , ` < 0.

 T
j=1−`
(4.2)

We take the function Wb to be a symmetric and continuous weight function with bounded
support of order b; see Chapter 7 of [38]. D̂T then defines an estimator of d by

9

dˆT (f )(t, s) =

ZZ
D̂T (t, s, u, v)f (u, v)dudv,

which further defines estimates of the eigenvalues and eigenfunctions of d satisfying
(4.3)

dˆT (ϕ̂i )(t, s) = λ̂i ϕ̂i (t, s).

In the simulations and application below, we take the weight function Wb to be the
simple Bartlett weight function, W1 (x) = (1 − |x|)1(|x| < 1), which is of order 1, and
h = dN 1/5 e. Under Assumption 2.1, we have the following result:
Theorem 4.1. Under Assumption 2.1, kD − D̂T k4 = oP (1).
In order to obtain consistent estimates of the first p elements of the spectrum of d, we
assume for the sake of simplicity that the first p eigenvalues of d are distinct, which
implies that the corresponding eigenspaces of the first p eigenvalues are one dimensional.
Assumption 4.1. We assume that there exists an integer p ≥ 1 satisfying that
λ1 > · · · > λp > λp+1 ≥ 0
where {λi , i ≥ 1} are defined in (2.2).
Assumption 4.1 could be relaxed by utilizing some of the ideas presented in [40], but
we do not pursue adapting those here. Under this assumption, the following result is
implied by Theorem 4.1 and the results in Section 6.1 of [17].
Corollary 4.1. Under Assumption 2.1 and 4.1, (3.2) holds.
4.2. Numerical implementation. Although so far we have presented results as if the
functional data at hand were observed on their entire domains, in practice the data will
consist of only discrete observations of the underlying functions. Let Xi (tj ) and Yi (tj ),
1 ≤ j ≤ R, denote the observed values of the functions Xi (t) and Yi (t), observed at the
common points {t1 , . . . , tR }. We assume here that each functional observation Xi (t) and
Yi (s) are observed at common points in their domains, as this matches our simulations
and data example below, although this could be easily relaxed.
It is straightforward to estimate the test statistics FT , FT,p , ZT , and ZT,p from the discrete
data and simple Riemann sum approximations to the inner products and norms, so long
as the eigenvalue and eigenfunction estimates λ̂i and ϕ̂i are given. However it is less
clear how to estimate the eigen-elements satisfying (4.3) from the discrete data. One
could in principle estimate the eigenvalues and eigenfunctions (λi , ϕi (t, s)), 1 ≤ i ≤ p
of d by calculating the spectrum of the 4-way tensor D̂T (ti , tj , tk , t` ), 1 ≤ i, j, k, ` ≤
R of dimension R4 , and employing linear interpolation to complete the eigenfunction
estimates, but this becomes computationally infeasible for even moderate values of R
due to the sheer dimension of the tensor, and the fact that this tensor is typically dense.
What we propose instead is a dimension reduction based approach involving the functional principal components (fPC’s) of the individual series {Xi } and {Yi }. Let
"
#
T
1X
(X` (ti ) − X̄(ti ))(X` (tj ) − X̄(tj )) : 1 ≤ i, j ≤ R
b
cX :=
T `=1

10

GREGORY RICE AND MARCO SHUM

be the sample covariance matrix of the discretized observations of the X sample, and define b
cY similarly. Calculating the spectrum of b
cX and b
cY yields eigenvalues {νX,1 , . . . , νX,R },
{νY,1 , . . . , νY,R } and eigenvectors,
{θX,1 , . . . , θX,R }, {θY,1 , . . . , θY,R }, the latter of which,
√
when multiplied by R, yield discrete approximations to the fPC’s of the sequences
{Xi } and {Yi }, θX,i (tj ), θY,i (tj ), 1 ≤ i, j ≤ R, respectively.
We then use the product basis of L2 [0, 1]4 generated by these functions to reduce the
dimension of D̂T . The projections of D̂T onto the product basis of the first q elements
each of θX,i (t), θY,i (t) may be stored in a 4-way tensor with q 4 elements, M, via
Z
Z
Mijkr = · · · D̂T (t, s, u, v)θX,i (t)θY,j (s)θX,k (u)θY,r (v)dtdsdudv.
Each of these elements can be estimated with a simple Riemann sum approximation.
The cutoff q must be selected by the user. We suggest using the total variance explained
approach for this: we take q so that θX,i (t) and θY,i (t), 1 ≤ i ≤ q, explains at least v%
of the variation in each series, where v is close to but strictly smaller than 1.
The eigenvalues {λ̂i , 1 ≤ i ≤ q 2 } and eigen-arrays {Φ̂i ∈ Rq×q , 1 ≤ i ≤ q 2 } of M satisfy
MΦ̂` = λ̂` Φ̂` , 1 ≤ ` ≤ q 2 , where MΦ̂` [i, j] =

q
X

Mijkr Φ̂` [k, r].

k,r=1

.
This eigenvalue problem may be solved numerically by solving for the eigenvalues/vectors
of a q 2 by q 2 square matrix that is “tiled” with the cross-sections of M, as is implemented
in svd.tensor function in the the package tensorA in R; see [45]. The eigenfunctions
ϕi (t, s) may then be estimated with
2

ϕ̂i (t, s) =

q
X

Φ̂i [j, k]θX,j (t)θY,k (s).

j,k=1

We use these estimates for the

λ̂0i s

and ϕ̂0i s in the simulations and application below.

5. Simulation Study
5.1. Outline. We now present the results of a small simulation experiment which aimed
to study the finite sample properties of the test statistics introduced above. All of the
simulations reported below were done using the R language ([43]). The number of potential experimental settings that could be considered to study Theorem 2.1 and Corollaries
3.1-3.3 is enormous. For the sake of brevity, we mainly focused on demonstrating that for
a somewhat rich class of data generating processes (DGP’s) exhibiting serial correlation,
the tests based on FT , FT,p , ZT , and ZT,p hold their size well, and that the eigenvalue
and eigenfunction estimation procedure explained in Section 4.2 is adequate for such hypothesis testing problems. Towards this goal, we considered the following basic structure
for generating synthetic data depending on the parameter α ∈ [0, 1]:
(5.1)

Xi (t) = αεc,i (t) + (1 − α)εx,i (t), Yi (t) = αεc,i (t) + (1 − α)εy,i (t), 1 ≤ i ≤ T

where εc,i (t), εx,i (t), and εy,i (t) are mutually independent sequences that satisfy either of
two models given below. In this case, the functional series Xi (t) and Yi (t) are correlated

11

through their common dependence on εc,i (t), with the strength of this dependence controlled by α. If α = 0 for instance, then the two sequences Xi and Yi are independent.
We took the error sequences to satisfy either:
IID: εc,i (t) = Wc,i (t), εx,i (t) = Wx,i (t), and εy,i (t) = Wy,i (t), where {Wc,i (t)},
{Wx,i (t)}, {Wy,i (t)} are mutually independent sequences of IID standard Brownian motions.
R
R
FAR(1): εc,iR(t) = Φ(t, s)εc,i (s)ds+Wc,i (t), εx,i (t) = Φ(t, s)εx,i (s)ds+Wx,i (t),
and εy,i (t) = Φ(t, s)εy,i (s)ds+Wy,i (t), where {Wc,i (t)}, {Wx,i (t)}, and {Wy,i (t)}
are defined above, and Φ(t, s) = min(t, s).
The sequence {(Xi , Yi )} satisfying (5.1) with error sequences satisfying either of the above
models satisfy Assumption 2.1. The choice of the Brownian motions for the innovation
sequences is partially motivated by our application to intraday returns data below, see
Figure 6.1. For each setting of T and α, the data was generated on an equally spaced
grid on [0, 1] with R = 100 points, and the simulation was repeated 1000 times for each
setting in order to calculate the empirical size and power curves presented below. We
chose q to be 3 in the estimation of procedure for for the eigenvalues and eigenfunctions
in Section 4, which in the vernacular of fPCA corresponds to a total variance explained
of around 93% on average for each sequence Xi and Yi . We considered larger values
of q, and found that it did not have much of an effect on the results, although the
computational time increases substantially as q increases. We similarly took n = 3 in
the definition of FT,p and ZT,p .
To study the size of the test of H0,1 based on FT and FT,3 , we generated data according
to (5.1) for each setting of the error sequence with α = 0, and performed tests with
nominal levels of 10%, 5%, and 1% of H0,1 : CXY = 0. The rejection rates from 1000
simulations from each statistic are reported in Table 5.1. One thing that was clear
based on these simulations was that the tests based on either statistic tended to be
somewhat oversized. This issue improves as T increases, as predicted by Corollaries
3.1 and 3.2, and all tests achieved quite good size once T ≥ 300. The presence of
serial correlation pronounces the size inflation of the tests, however this is fairly well
controlled by the simple weight function-bandwidth smoothed periodogram estimator.
This could be improved further by increasing the bandwidth at the cost of increased
computational time, according to unreported simulations. The fact that the empirical
levels are quite close to nominal and improve with increasing T is indicative that the
eigenvalue/eigenfunction estimation described in Section 4.2 is performing adequately
for this application. We aim in future work to improve the size properties of the test
under serial correlation by implementing a data driven bandwidth selection routine, and
to consider alternative methods for estimating the eigenvalues/eigenfunctions of d.
We also studied the empirical power of FT and FT,p by testing H0,1 : CXY = 0 with
data following (5.1) and α increasing from 0 to 1. The results of this simulation are
reported in the case FAR(1) errors as power curves in Figure 5.1. We observed that
both tests were powerful for large enough values of α, although it is difficult to quantify
their power without a comparable test. Interestingly, the simple norm based test based
on FT possessed better power than FT,3 in all the examples that we considered.
Regarding the test of H0,2 based on ZT and ZT,3 , we only present the results of an
empirical size study. In this case we applied each statistic to data generated according
to (5.1) with α = 0 and α = 0.5 for T = 100, 300 and 1000. The empirical size from

12

GREGORY RICE AND MARCO SHUM

Statistic: FT
IID
T

10%

5%

FAR(1)
1%

50 0.133 0.073 0.013
100 0.129 0.061 0.015
300 0.114 0.049 0.008
Statistic:
IID
T

10%

5%

1%

10%

5%

1%

0.155 0.087 0.018
0.148 0.069 0.020
0.139 0.077 0.018
FT,3
FAR(1)
10%

5%

1%

50 0.142 0.078 0.026
0.171 0.096 0.025
100 0.118 0.073 0.019
0.189 0.087 0.034
300 0.121 0.041 0.012
0.134 0.076 0.021
Table 5.1. Empirical sizes with nominal levels of 10%, 5%, and 1% for
a test of H0,1 : CXY = 0 where the data was generated according to (5.1)
with α = 0.

Figure 5.1. Power curves as a function of α ranging from 0 to 0.8 at
increments of 0.1, for a test H0,1 : CXY = 0 with level of 5% applied to
(Xi , Yi ) following (5.1) with FAR(1) errors using FT (left hand panel), and
FT,3 (right hand panel).
1000 independent simulations for each setting is presented in Table 5.2. We saw that in
general the change point tests based on ZT and ZT,3 had quite good size, even in the
presence of serial correlation. Changing α had no effect on the size of these tests, as
expected.

13

Statistic: ZT
IID
T

α

100

0
0.5
300 0
0.5
1000 0
0.5

T

α

10%
0.093
0.125
0.102
0.115
0.133
0.123

10%

5%

1%

0.048 0.007
0.066 0.020
0.047 0.008
0.060 0.020
0.071 0.006
0.058 0.012
Statistic: ZT,3
IID
5%

1%

FAR(1)
10%

5%

1%

0.149
0.094
0.122
0.121
0.146
0.139

0.087
0.046
0.079
0.057
0.063
0.075

0.031
0.006
0.025
0.015
0.010
0.012

FAR(1)
10%

5%

1%

100

0 0.108 0.060 0.016
0.166 0.098 0.027
0.5 0.114 0.070 0.025
0.132 0.058 0.016
300 0 0.132 0.073 0.021
0.132 0.067 0.011
0.5 0.121 0.056 0.015
0.104 0.062 0.013
1000 0 0.143 0.082 0.024
0.102 0.058 0.012
0.5 0.126 0.067 0.014
0.094 0.052 0.010
Table 5.2. Empirical sizes with nominal levels of 10%, 5%, and 1% for
(T )
(1)
a test of H0,2 : CXY = . . . = CXY where the data was generated according
to (5.1) with α = 0 and α = 0.5.

6. Application to cumulative intraday returns
A natural example of functional time series data are those derived from densely recorded
asset price data, such as intraday stock price data. Recently there has been an upsurge
in quantitative research focused on analyzing the information contained within curves
constructed from such data; we refer the reader to [6], [46], [16], [30], and [29]. Price
curves associated with popular companies are commonly displayed at websites like yahoo.com/finance, and the “patterns” in them, or lack thereof, are of apparent interest
to day traders.
The specific data that we consider was obtained from www.nasdaq.com, and consists of
1 minute resolution closing prices of a single share of Microsoft (ticker MSFT), and Exxon
Mobile (ticker XOM) stock from January 2nd to December 31st, 2001, which comprises
data from 248 trading days (T = 248) with R = 389 observations per day. We aimed
to apply the methods introduced in Sections 3 to study the cross-covariance structure
between these two companies stock prices on the intraday scale. Let PM,i (tj ), and PX,i (tj )
i = 1, . . . , T, j = 1, . . . , R, denote the prices of Microsoft and Exxon mobile stock on day
i at intraday time tj , respectively. The first three price curves of each series constructed
from the raw price data and linear interpolation are displayed in the left hand panel of
Figure 6.1. The functional time series of price curves are evidently nonstationary due to
frequent level shifts and volatility, and hence we considered the following transformation
of these curves akin to taking the log returns for scalar price data:

14

GREGORY RICE AND MARCO SHUM

Figure 6.1. The left hand panel displays three intraday price curves
derived from the one-minute resolution closing prices of a single share of
XOM and MSFT. The right hand panel displays the corresponding CIDR
curves.
Definition 1. Suppose Pi (tj ), i = 1, . . . , T, j = 1, . . . , R, is the price of a financial asset
at time tj on day i. The functions
ri (tj ) = 100[ln Pi (tj ) − ln Pi (t1 )],

j = 1, 2, . . . , R,

i = 1, . . . , T,

are called the cumulative intraday returns (CIDR’s).
Since the logarithm is increasing, the CIDR curves have nearly the same shape as the
original daily price curves, but the assumption of stationarity is much more plausible for
these curves. According to their definition, the CIDR’s always start from zero, so level
stationarity is enforced, and taking the logarithm helps reduce potential scale inflation.
The stationarity of CIDR curves derived from intraday stock price data was argued
empirically in [25].
Let Xi (t) = rM,i (t) and Yi (t) = rX,i (t), 1 ≤ i ≤ T , denote the CIDR curves derived from
the Microsoft and Exxon Mobile stock price data and linear interpolation, respectively.
The first three CIDR curves of each series are plotted in the right hand panel of Figure
6.1.
Before estimating and measuring the significance of the cross-covariance between these
functional time series of CIDR’s, we first tested for its homogeneity within the sample, a
test of H0,2 , using the test statistic ZT . For this analysis, we took q = 3 for the method
outline in Section 4.2, which corresponded to approximately 94% variance explained in
each series. A “CUSUM chart” of c(t) = T kĈXY (·, ·, t/T ) − (t/T )ĈXY (·, ·, 1)k22 versus
t is given in Figure 6.2 with the corresponding 10%, 5% and 1% significance levels of
the estimated limiting distribution in Corollary 3.1. The statistic ZT far exceeded the
1% level, which indicates that there is strong evidence that the covariance relationship
is heterogenous within the sample. Also apparent in the plot, the largest difference

15

between the partial sample cross-covariance estimate occurs on March 22nd, 2001. We
segmented the data into two sub-samples before and after this point of lengths T1 = 56
and T2 = 192, respectively, and again tested for the homogeneity of the cross-covariance
within each sub-sample. In both sub-samples the homogeneity could not be rejected
with any significance. Interestingly, the sample including the intraday XOM and MSFT
returns data before and after the terrorist attacks on September 11, 2001 did not seem
to exhibit a change point in the cross-covariance relationship. The date March 22nd,
2001 does however seem to be in the proximity of some fairly important events relative
to the world oil market. On March 17th, OPEC had announced that they were cutting
oil production by 4%, and just three days later the largest off shore oil rig in the world,
Petrobras 36, sank off of the coast of Brazil.
In order to measure the significance of the estimates of the cross-covariance before and
after the change point, we applied a test of H0,1 : CXY = 0 to each sub-sample. The null
hypothesis was strongly rejected in both cases with p-values smaller than 0.000. The
cross-covariance surfaces are displayed in Figure 6.3, from which we can see that the
shape of the cross-covariance is quite different before and after the initial change point:
the surface goes from being negative to predominantly positive.

16

GREGORY RICE AND MARCO SHUM

Figure 6.2.
The top panel displays the “CUSUM chart” of
kĈXY (·, ·, t/T ) − (t/T )ĈXY (·, ·, 1)k22 versus t calculated from the XOM and
MSFT price data from 2001. The maximum value ZT is significant to the
0.01 level, and is achieved at the value t corresponding to March 22nd.
The lower left and lower right panels show similar CUSUM charts for the
sub samples of data before and after this point, respectively.
References
[1] J. Aston and C. Kirch. Detecting and estimating epidemic changes in dependent functional data.
Journal of Multivariate Analysis, 109:204–220, 2012.
[2] A. Aue, S. Hörmann, L. Horváth, and M. Reimherr. Break detection in the covariance structure of
multivariate time series models. The Annals of Statistics, 37:4046–4087, 2009.

17

Figure 6.3. The left and right hand panels display the estimator ĈXY
before and after the date 3/22/2001, respectively. Both of these surface
estimates were measured to be significantly different from zero, as the test
of H0,1 : CXY = 0 based on FT yielded p-values less than 0.000 in both
cases.
[3] A. Aue, L. Horváth, and D. F. Pellatt. Functional generalized autoregressive conditional heteroskedasticity. Technical report, University of California, Davies, 2015.
[4] A. Aue, D. D. Norinho, and S. Hörmann. On the prediction of stationary functional time series.
Journal of the American Statistical Association, 110:378–392, 2014.
[5] C. Baker. Joint measures and cross-covariance operators. Transactions of the american mathematical society, 186:273–289, 1973.
[6] O. E. Barndorff-Nielsen and N. Shephard. Econometric analysis of realized covariance: High frequency based covariance, regression and correlation in financial economics. Econometrica, 72:885–
925, 2004.
[7] I. Berkes, L. Horváth, and G. Rice. Weak invariance principles for sums of dependent random
functions. Stochastic Processes and Their Applications, 123:385–403, 2013.
[8] D. Bosq. Linear Processes in Function Spaces. Springer, New York, 2000.
[9] D. R. Brillinger. Time Series: Data Analysis and Theory. Holt, New York, 1975.
[10] B. Buccia and M. Wendler. Change-point detection and bootstrap for hilbert space valued random
fields. Journal of Multivariate Analysis, 155:344–368, 2017.
[11] J-M. Chiou, Y-T. Chen, and Y-F. Yang. Multivariate functional principal component analysis: A
normalization approach. Statistica Sinica, 24:15711596, 2014.
[12] J-M. Chiou and H-G. Müller. Linear manifold modeling of multivariate functional data. Journal of
the Royal Statistical Society Series B, 76:605–626, 2013.
[13] H. Dette, W. Wu, and Z. Zhou. Change point analysis of second order characteristics in nonstationary time series. Technical Report, 2015.
[14] J. Dubin and H-G. Müller. Dynamical correlation for multivariate longitudinal data. Journal of the
American Statistical Association, 100:872–881, 2005.
[15] F. Ferraty and P. Vieu. Nonparametric Functional Data Analysis: Theory and Practice. Springer,
New York, 2006.
[16] R. Gabrys, L. Horváth, and P. Kokoszka. Tests for error correlation in the functional linear model.
Journal of the American Statistical Association, 105:1113–1125, 2010.

18

GREGORY RICE AND MARCO SHUM

[17] I. Gohberg, S. Golberg, and M. A. Kaashoek. Classes of Linear Operators, volume 49 of Operator
Theory: Advances and Applications. Birkhaüser, 1990.
[18] K. El Himdi, R. Roy, and P. Duchesne. Tests for non–correlation of two multivariate time series:
a nonparametric approach. IMS Lecture Notes Monogr. Ser., Mathematical Statistics and Applications: Festschrift for Constance van Eeden, 42:397–416, 2003.
[19] S. Hörmann, L. Kidziński, and M. Hallin. Dynamic functional principal components. Journal of the
Royal Statistical Society, Ser. B, 77:319–348, 2015.
[20] S. Hörmann, L. Kidziński, and P. Kokoszka. Estimation in functional lagged regression. Journal of
Time Series Analysis, 36:541–561, 2015.
[21] S. Hörmann and P. Kokoszka. Weakly dependent functional data. The Annals of Statistics, 38:1845–
1884, 2010.
[22] S. Hörmann and P. Kokoszka. Functional time series. In C. R. Rao and T. Subba Rao, editors,
Time Series, volume 30 of Handbook of Statistics. Elsevier, 2012.
[23] L. Horváth, M. Hušková, and G. Rice. Testing independence for functional data. Journal of Multivariate Analysis, 117:100–119, 2013.
[24] L. Horváth, P. Kokoszka, and R. Reeder. Estimation of the mean of functional time series and a
two sample problem. Journal of the Royal Statistical Society (B), 74, 2012.
[25] L. Horváth, P. Kokoszka, and G. Rice. Testing stationarity of functional time series. Journal of
Econometrics, 179:66–82, 2014.
[26] L. Horváth and G. Rice. Testing for independence between two functional time series. Journal of
Econometrics, 189:371–382, 2015.
[27] M. Jirak. On weak invariance principles for sums of dependent random functionals. Statistics and
Probability Letters, 83:2291–2296, 2013.
[28] V. Kargin and A. Onatski. Curve forecasting by functional autoregression. Journal of Multivariate
Analysis, 99:2508–2526, 2008.
[29] P. Kokoszka and M. Reimherr. Asymptotic normality of the principal components of functional
time series. Stochastic Processes and their Applications, 123:1546–1562, 2013.
[30] H-G. Müller, R. Sen, and U. Stadtmüller. Functional data analysis for volatility. Journal of Econometrics, 165:233–245, 2011.
[31] V. M. Panaretos, D. Kraus, and J. H. Maddocks. Second-order comparison of Gaussian random
functions and the geometry of DNA minicircles. Journal of the American Statistical Association,
105:670–682, 2010.
[32] V. M. Panaretos and S. Tavakoli. Fourier analysis of stationary time series in function space. Annals
of Statistics, 41(2):568–603, 2012.
[33] V. M. Panaretos and S. Tavakoli. Detecting and localizing differences in functional time series
dynamics: A case study in molecular biophysics. Journal of the American Statistical Association,
111(515):1020–1035, 2016.
[34] E. Paparoditis and T. Sapatinas. Bootstrap-based k sample testing for functional data. Biometrika,
103:727–733, 2016.
[35] A. Petersen and H-G. Müller. Fréchet integration and adaptive metric selection for interpretable
covariance of multivariate functional data. Biometrika, 103:103–120, 2016.
[36] T. Pham and V. Panaretos. Methodology and convergence rates for functional time series regression.
Technical Report, 2017.
[37] D. Pigoli, J. Aston, I. Dryden, and P. Secchi. Distances and inference for covariance operators.
Biometrika, 101:409–422, 2014.
[38] M. B. Priestley. Spectral Analysis and Time Series. Academic Press, 1981.
[39] J. O. Ramsay and B. W. Silverman. Functional Data Analysis. Springer, 2005.
[40] M. Reimherr. Functional regression with repeated eigenvalues. Statistics and Probability Letters,
107:62–70, 2015.
[41] N. Serban, A-M. Staicu, and R-J. Carroll. Multilevel cross-dependent binary longitudinal data.
Biometrics, 69:903–913, 2013.
[42] O. Sharipov, J. Tewes, and M. Wendler. Sequential block bootstrap in a hilbert space with application to change point analysis. Canadian Journal of Statistics, 44:300–322, 2016.

19

[43] R Core Team. R: A Language and Environment for Statistical Computing., 2015. R package version
≥ 2.2.0.
[44] L. Torgovitski. Darling-Erdös-type CUSUM-procedure for functional data II. Working paper, University of Cologne, 2016.
[45] K. Gerald van den Boogaart. tensorA: Advanced tensors arithmetic with named indices, 2010. R
package version ≥ 2.2.0.
[46] Y. Wang and J. Zou. Vast volatility matrix estimation for high–frequency financial data. The Annals
of Statistics, 38:953–978, 2010.
[47] D. Wied, M. Arnold, N. Bissantz, and D. Ziggel. A new fluctuation test for constant variances with
applications to finance. Metrika, 75:1111–1127, 2012.
[48] D. Wied, W. Krämer, and H. Dehling. Testing for a change in correlation at an unknown point in
time using an extended functional delta method. Econometric Theory, 28:570–589, 2012.
[49] X. Zhang. White noise testing and model diagnostic checking for functional time series. Journal of
Econometrics, 194:76–95, 2016.
[50] X. Zhang and X. Shao. Two sample inference for the second-order property of temporally dependent
functional data. Bernoulli, 21:909–929, 2015.
[51] L. Zhou, J. Huang, and R-J. Carroll. Joint modelling of paired sparse functional data using principal
components. Biometrika, 95:601–619, 2008.

Appendix A. Proofs of technical results
A.1. Proof of Theorem 2.1. In order to prove Theorem 2.1, we may assume without
loss of generality that µX = µY = 0. Let
bT xc
1 X
Xi (t)Yi (s).
C̃XY (t, s, x) =
T i=1

Lemma A.1. Under Assumption 2.1,
 
1
sup kĈXY (·, ·, x) − C̃XY (·, ·, x)k2 = OP
.
T
0≤x≤1
Proof. According to the definitions of ĈXY and C̃XY and the triangle inequality,
(A.1)
bT xc
1 X
kĈXY (·, ·, x) − C̃XY (·, ·, x)k2 =
(Xi − X̄) ⊗ (Yi − Ȳ ) − Xi ⊗ Yi
T i=1

2

bT xc

(A.2)

=

1 X
(−X̄ ⊗ Yi − Xi ⊗ Ȳ + X̄ ⊗ Ȳ )
T i=1
bT xc

bT xc

1 X
1 X
≤
X̄ ⊗ Yi +
Xi ⊗ Ȳ
T i=1
T i=1

(A.3)

2

(A.4)

= G1 (x) + G2 (x) + G3 (x).

For the term G1 (x), we have that

(A.5)

G1 (x) =

bT xc
bT xc
X
1 X
1
X̄ ⊗ Yi = kX̄k1
Yi .
T i=1
T
i=1
2

1

2

2

bT xc
1 X
+
X̄ ⊗ Ȳ
T i=1

2

20

GREGORY RICE AND MARCO SHUM

Assumption 2.1 implies that both the series Xi and Yi satisfy the conditions of Theorem
3.3 of Berkes et al [7], from which it follows that

kX̄k1 = OP

1
√
T

bT xc


, and sup
0≤x≤1

X
i=1

Yi = OP

√ 
T .

1

This combined with (A.5) implies that sup0≤x≤1 G1 (x) = OP (1/T ). Parrel arguments
show that sup0≤x≤1 G2 (x) = OP (1/T ) and sup0≤x≤1 G3 (x) = OP (1/T ), from which the
result follows in light of (A.1).

Lemma A.2. Under Assumption 2.1 and if q = p/2 with p defined in Assumption 2.1,
then
(EkXi ⊗ Yi − Xi,m ⊗ Yi,m kq2 )1/q = O(m−γ ),
where γ = min(α, β), and Xi,m and Yi,m are defined in Assumption 2.1.
Proof. We have according to the triangle inequalities in L2 [0, 1]2 and for (E(·)q )1/q that
(A.6)
(EkXi ⊗ Yi − Xi,m ⊗ Yi,m kq2 )1/q = (EkXi ⊗ Yi − Xi ⊗ Yi,m + Xi,m ⊗ Yi − Xi,m ⊗ Yi,m kq2 )1/q
≤ (E(kXi ⊗ Yi − Xi ⊗ Yi,m k2 + kXi,m ⊗ Yi − Xi,m ⊗ Yi,m k2 )q )1/q
≤ (EkXi ⊗ Yi − Xi ⊗ Yi,m kq2 )1/q + (EkXi,m ⊗ Yi − Xi,m ⊗ Yi,m kq2 )1/q
= (EkXi kq1 kYi − Yi,m kq1 )1/q + (EkYi kq1 kXi − Xi,m kq1 )1/q .
According to the Cauchy-Schwarz inequality and stationarity, we have that
(A.7)
1/2
1/2
1/2
1/2
EkXi kq1 kYi − Yi,m kq1 ≤ (EkXi k2q
E(kYi − Yi,m k2q
= (EkX0 k2q
E(kY0 − Y0,m k2q
.
1 )
1 )
1 )
1 )
One obtains a similar bound with the roles of X and Y swapped, which combined with
the last line of (A.6) we get that (with p = 2q)
(EkXi ⊗Yi −Xi,m ⊗Yi,m kq2 )1/q ≤ (EkX0 kp1 )1/p (EkY0 −Y0,m kp1 )1/p +(EkY0 kp1 )1/p (EkX0 −X0,m kp1 )1/p .
Since both (EkX0 kp1 )1/p and (EkY0 kp1 )1/p are finite, and according to Assumption 2.1
1/p
1/p
(EkXi − Xi,m kp
= O(m−α ), and (EkYi − Yi,m kp
= O(m−β ), we obtain that
(EkXi ⊗ Yi − Xi,m ⊗ Yi,m kq2 )1/q = O(m−γ ),
where γ = min(α, β) as needed.

Proof of Theorem 2.1. According to Lemma A.1, it is enough to prove Theorem 2.1
for the process C̃XY (t, s, x). C̃XY (t, s, x) − (bT xc/T )CXY (t, s) is the partial sum of
the random functions ξi (t, s) = Xi (t)Yi (s) − EX0 (t)Y0 (s) which form a mean zero and
stationary sequence in L2 [0, 1]2 . Let ξi,m (t, s) = Xi,m (t)Yi,m (s) − EX0 (t)Y0 (s), denote an
m-dependent approximation to ξi . It follows directly from A.2 that
(Ekξi − ξi,m kq2 )1/q = O(m−γ ),

21

where γ > 1, and this implies that
∞
X

(Ekξi − ξi,m kq2 )1/q < ∞.

m=1

It now follows from Theorem 1 of [27] that if
√
ΞT (t, s, x) := T



bT xc
CXY (t, s)
C̃XY (t, s, x) −
T



bT xc
1 X
=√
ξi (t, s),
T i=1

then
ZZ
(A.8)

sup

(ΞT (t, s, x) − ΓT (t, s, x))2 dtds = oP (1),

0≤x≤1

where ΓT has the properties attributed in Theorem 2.1, which establishes the theorem.

To see how Corollaries 3.1-3.3 follow from this result, we note that for all T ,
D

{ΓT (t, s, x) : t, s, x ∈ [0, 1]} = {Γ0 (t, s, x) : t, s, x ∈ [0, 1]},
where
Γ0 (t, s, x) =

∞
X

1/2

λ` W` (x)ϕ` (t, s),

`=1

with {W` (t), t ∈ [0, 1], ` ∈ N} being iid standard Brownian motions, and (λ` ϕ` ) being
the eigenvalues and eigenfunctions of the operator d. This follows from simply calculating
the mean and covariance function of Γ0 and using Mercer’s theorem.
Proof of Theorem 3.1. The first part of the theorem follows directly from the ergodic
theorem in Hilbert spaces; see for example Appendix A of [23]. In order to get the second
√
L2 [0,1]2
result, we have using the assumption that T (CXY − C0 ) → CA and Theorem 2.1
that, in case of FT ,
FT = T kĈXY (·, ·, 1) − CXY + CXY − C0 k22
= T h(ĈXY (·, ·, 1) − CXY ) + (CXY − C0 ), (ĈXY (·, ·, 1) − CXY ) + (CXY − C0 )i2
√
√
= T kĈXY (·, ·, 1) − CXY k22 + 2h T (ĈXY (·, ·, 1) − CXY ), T (CXY − C0 )i2 + T kCXY − C0 k22
∞
∞
X
X
D
1/2
2
→
λi Ni + 2
λi hCA , ϕi iNi + kCA k22 , (T → ∞),
i=1

i=1

as needed. The result for FT,p follows similarly.



22

GREGORY RICE AND MARCO SHUM

A.2. Proof of Theorem 4.1. Let
∞
X

 
`
γ̃` (t, s, u, v),
D̃T (t, s, u, v) =
Wb
h
`=−∞

(A.9)

with h defined in (4.2), Wb being a bounded, symmetric, and continuous weight function
of order b, and

γ̃` (t, s, u, v) =


T −`

1X



(Xj (t)Yj (s) − EX0 (t)Y0 (s)) (Xj+` (t)Yj+` (s) − EX0 (t)Y0 (s)) ,


 T j=1

`≥0


T

1 X



(Xj (t)Yj (s) − EX0 (t)Y0 (s)) (Xj+` (t)Yj+` (s) − EX0 (t)Y0 (s)) , ` < 0.

 T
j=1−`

Lemma A.3. Under Assumption 2.1,
kD̂ − D̃k4 = oP (1).
Proof. According to the definition of γ̂` , when ` ≥ 0, we have that
(A.10)
T −`
1X
(X̄j (t)Ȳj (s) − ĈXY (t, s, 1))(X̄j+` (u)Ȳj+` (v) − ĈXY (u, v, 1))
γ̂` (t, s, u, v) =
T j=1
T −`
1X
(Xj (t)Yj (s) − EX0 (t)Y0 (s)) − X̄(t)Yj (s) + X̄(t)Ȳ (s) − Xj (t)Ȳ (s)
=
T j=1
!

+ (EX0 (t)X0 (s) − Ĉ(t, s, 1)) × Xj+` (u)Yj+` (v) − EX0 (u)Y0 (v) − X̄(u)Yj+` (v)
!
+ X̄(u)Ȳ (v) − Xj+` (u)Ȳ (v) + (EX0 (u)X0 (v) − Ĉ(u, v, 1))
= γ̃(t, s, u, v) +

24
X

Rp,` (t, s, u, v),

p=1

where the terms Rp,` represent the remaining 24 terms in the definition of γ̂` obtained by
completing the multiplication that results in 25 terms, the first of which corresponds to
γ̃` . We now aim to show that for all p and `, EkRp,` k4 = O(T −1/2 ). For R1,` , we have by
the triangle inequality, a couple of applications of the Cauchy-Schwarz inequality, and
the assumed stationarity that
(A.11)
T −`
X
1
EkR1,` k4 = E
(Xj (t)Yj (s) − EX0 (t)Y0 (s))X̄(u)Yj+` (v)
T
j=1

4

23
T −`
1X
≤
Ek(Xj (t)Yj (s) − EX0 (t)Y0 (s))X̄(u)Yj+` (v)k4
T j=1

=

T −`
1X
Ek(Xj (t)Yj (s) − EX0 (t)Y0 (s))k2 kX̄k1 kYj+` k1
T j=1

T −`
1X
≤
(Ek(Xj (t)Yj (s) − EX0 (t)Y0 (s))k22 )1/2 (E(kX̄k1 kYj+` k1 )2 )1/2
T j=1
T −`
1X
(Ek(Xj (t)Yj (s) − EX0 (t)Y0 (s))k22 )1/2 (EkX̄k41 )1/4 (EkYj+` k41 )1/4
≤
T j=1

T −`
(Ek(X0 (t)Y0 (s) − EX0 (t)Y0 (s))k22 )1/2 (EkX̄k41 )1/4 (EkY0 k41 )1/4 .
T
It follows from Proposition 3.1 of [44] that under Assumption 2.1, (EkX̄k41 )1/4 = O(T −1/2 ),
and also according to Assumption 2.1, (EkY0 k41 )1/4 = O(1) and (Ek(X0 (t)Y0 (s)−EX0 (t)Y0 (s))k22 )1/2 =
O(1). This implies with (A.11) that EkR1,` k4 = O(T −1/2 ). Similar arguments using the
fact that (EkȲ k41 )1/4 = O(T −1/2 ) and kEX0 (t)Y0 (s) − ĈXY (t, s, 1)k2 = O(T −1/2 ) can
be applied to get that EkRp,` k = O(T −1/2 ), 2 ≤ p ≤ 24. This implies via the triangle
inequality and (A.10) that
=

Ekγ̂` − γ̃` k4 = O(T −1/2 ),
from which it follows that

(A.12)

E

∞
X
`=1

Wb

!
!
∞
X
`
`
Wb
(γ̂` − γ̃` ) ≤
Ek(γ̂` − γ̃` )k4 = O(hT −1/2 ) = o(1),
h
h
`=1
4

by (4.2) and the fact that Wb has bounded support. A parallel argument may be used
to show that for ` < 0, that Ekγ̂` − γ̃` k4 = O(T −1/2 ), from which we obtain that
(A.13)

E

−1
X
`=−∞

Wb

!
`
(γ̂` − γ̃` ) = O(hT −1/2 ) = o(1).
h
4

(A.12) and (A.13) together with the triangle inequality imply that
!
!
∞
−1
X
X
`
`
(γ̂` − γ̃` ) +E
Wb
(γ̂` − γ̃` ) = o(1),
EkD̂ − D̃k4 ≤ E
Wb
h
h
`=1
`=−∞
4

and this implies the result with Markov’s inequality.

4



Lemma A.4. Under Assumption 2.1,
kD̃ − Dk4 = oP (1).
Proof. γ̃` is the autocovariance estimator based on a sample of size T at lag ` of the
mean zero random functions {ξi (t, s), i ∈ Z}, where, recalling from above, ξi (t, s) =
Xi (t)Yi (s) − EX0 (t)Y0 (s). With ξi,m (t, s) = Xi,m (t)Yi,m (s) − EX0 (t)Y0 (s), we have from
A.2 and Lyapounov’s inequality that

24

GREGORY RICE AND MARCO SHUM

lim m(Ekξi − ξi,m k22 )1/2 = 0,

m→∞

and

∞
X

(Ekξi − ξi,m k22 )1/2 < ∞.

m=1

This shows that ξi satisfy the conditions of Theorem 2 from [24], which implies the result.

Proof of Theorem 4.1. The Theorem follows directly from Lemmas A.3 and A.4.



A.3. Consistency of ZT and ZT,p . In order for the tests based on ZT and ZT,p to be
consistent, we only require that the series before and after the change is ergodic, which
we quantify by the following assumption.
Assumption A.1. Under HA,2 , we assume that
∗

k
1 X
(Xi (t) − µX (t))(Yi (s) − µY (s)) − C1 (t, s) = oP (1),
k ∗ i=1
2

and

T
X
1
(Xi (t) − µX (t))(Yi (s) − µY (s)) − C2 (t, s) = oP (1).
T − k ∗ i=k∗ +1
2

Theorem A.1. Under HA,2 and Assumption A.1 holds, then
P

ZT −→ ∞.
If in addition C1 and C2 in HA,2 satisfy hC1 − C2 , ϕi i2 6= 0 for some 1 ≤ i ≤ p, then
P

ZT,p −→ ∞.
Assumption A.1 implies that the covariance estimators of the data before and after
the change are consistent with their underlying population quantities. This would hold
if, for example, the function gXY which defines how the series (Xi , Yi ) depends on the
(1)
(2)
underlying innovation sequence were to change from gXY to gXY in such a way that HA,2
held, since then the series before and after the change would still be ergodic. The proof
of Theorem A.1 is straightforward, so the details are omitted.
Department of Statistics and Actuarial Science, University of Waterloo, Waterloo,
ON, Canada
Department of Statistics and Actuarial Science, University of Waterloo, Waterloo,
ON, Canada

