Inference for the cross-covariance operator of stationary functional time series

TL;DR

This paper develops statistical inference methods for the cross-covariance operator of stationary functional time series, including hypothesis testing and change-point detection, with practical application to stock price data.

Contribution

It introduces a functional central limit theorem-based approach for inference on the cross-covariance operator and proposes a spectral density estimation method using functional PCA.

Findings

Method performs well in simulations

Effective change-point detection procedure

Application to NYSE stock data demonstrates utility

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 change-point 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.

Tables

Table 1. Table 5.1. Empirical sizes with nominal levels of 10%, 5%, and 1% for a test of H 0 , 1 : C X ​ Y = 0 : subscript 𝐻 0 1 subscript 𝐶 𝑋 𝑌 0 H_{0,1}:\ C_{XY}=0 where the data was generated according to ( 5.1 ) with α = 0 𝛼 0 \alpha=0 .

Statistic: $F_T$
	IID				FAR(1)
T	10%	5%	1%		10%	5%	1%
50	0.133	0.073	0.013		0.155	0.087	0.018
100	0.129	0.061	0.015		0.148	0.069	0.020
300	0.114	0.049	0.008		0.139	0.077	0.018
Statistic: $F_{T,3}$
	IID				FAR(1)
T	10%	5%	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 2. Table 5.2. Empirical sizes with nominal levels of 10%, 5%, and 1% for a test of H 0 , 2 : C X ​ Y ( 1 ) = … = C X ​ Y ( T ) : subscript 𝐻 0 2 superscript subscript 𝐶 𝑋 𝑌 1 … superscript subscript 𝐶 𝑋 𝑌 𝑇 H_{0,2}:\ C_{XY}^{(1)}=\ldots=C_{XY}^{(T)} where the data was generated according to ( 5.1 ) with α = 0 𝛼 0 \alpha=0 and α = 0.5 𝛼 0.5 \alpha=0.5 .

Statistic: $Z_T$
		IID				FAR(1)
T	$\alpha$	10%	5%	1%		10%	5%	1%
100	0	0.093	0.048	0.007		0.149	0.087	0.031
	0.5	0.125	0.066	0.020		0.094	0.046	0.006
300	0	0.102	0.047	0.008		0.122	0.079	0.025
	0.5	0.115	0.060	0.020		0.121	0.057	0.015
1000	0	0.133	0.071	0.006		0.146	0.063	0.010
	0.5	0.123	0.058	0.012		0.139	0.075	0.012
Statistic: $Z_{T,3}$
		IID				FAR(1)
T	$\alpha$	10%	5%	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

Equations

$$C_{XY}(t,s)=\text{mbox}cov(X_0(t),Y_0(s))=E[(X_0(t)-{\mu }_X(t))(Y_0(s)-{\mu }_Y(s))],$$

$$c_{XY}(f)(t)=\int C_{XY}(t,s)f(s)ds.$$

$${\hat{C}}_{XY}(t,s,x)=\frac{1}{T}\sum _{i=1}^{\lfloor Tx\rfloor }(X_i(t)-\bar{X}(t))(Y_i(s)-\bar{Y}(s)),$$

$$\bar{X}(t)=\frac{1}{T}\sum _{i=1}^T{X}_i(t),\text{mbox}and\bar{Y}(t)=\frac{1}{T}\sum _{i=1}^T{Y}_i(t).$$

$$\big{(}E\|X_{i}-X_{i,m}\|^{p}\big{)}^{1/p}=O(m^{-\alpha}),\mbox{ and }\big{(}E\|Y_{i}-Y_{i,m}\|^{p}\big{)}^{1/p}=O(m^{-\beta})$$

$$\underset{0\le x\le 1}{\sup }\int \int {\left(\sqrt{T}\left({\hat{C}}_{XY}(t,s,x)-\frac{\lfloor Tx\rfloor }{T}{C}_{XY}(t,s)\right)-{\Gamma }_T(t,s,x)\right)}^{2}dtds=o_P(1),$$

$$\text{mbox}cov({\Gamma }_T(t,s,x),{\Gamma }_T(u,v,y))=\min (x,y)D(t,s,u,v),$$

$$D(t,s,u,v)=\sum _{\ell =-\infty }^{\infty }\text{mbox}cov((X_0(t)-{\mu }_X(t))(Y_0(s)-{\mu }_Y(s)),(X_{\ell }(u)-{\mu }_X(u))(Y_{\ell }(v)-{\mu }_Y(v)).$$

$$d(f)(t,s)=\int \int D(t,s,u,v)f(u,v)dudv,$$

$$d({\phi }_i)(t,s)={\lambda }_i{\phi }_i(t,s).$$

$$H_{0,1}:C_{XY}=C_0\text{mbox}versus{H}_{A,1}:C_{XY}\ne C_0,$$

$$H_{0,2}:C_{XY}^{(1)}=C_{XY}^{(2)}=\cdots =C_{XY}^{(T)},\text{mbox}versus$$

$$H_{A,2}:C_1=C_{XY}^{(1)}=\cdots =C_{XY}^{(k^{\ast })}\ne C_{XY}^{(k^{\ast }+1)}=\cdots =C_{XY}^{(T)}=C_2,\text{mbox}forsome{k}^{\ast }=\lfloor T\theta \rfloor ,\theta \in (0,1),$$

$$F_T=T\Vert {\hat{C}}_{XY}(\cdot ,\cdot ,1)-C_0{\Vert }_2^2.$$

$$F_T\stackrel{D}{⟶}\sum _{i=1}^{\infty }{\lambda }_i{\mathcal{N}}_i^2,\text{mbox}asT\to \infty ,$$

$$⟨\sqrt{T}({\hat{C}}_{XY}(\cdot ,\cdot ,1)-C_{XY}),{\phi }_i{⟩}_2,1\le i\le p$$

$$F_{T,p}=\sum _{i=1}^p\frac{⟨\sqrt{T}({\hat{C}}_{XY}(\cdot ,\cdot ,1)-C_0),{\hat{\phi }}_i{⟩}_2^2}{{\hat{\lambda }}_i},$$

$$\underset{1\le i\le p}{\max }\Vert {\hat{\phi }}_i-{\hat{c}}_i{\phi }_i{\Vert }_2=o_P(1),\text{mbox}and\underset{1\le i\le p}{\max }|{\hat{\lambda }}_i-{\lambda }_i|=o_P(1),$$

$$F_{T,p}\stackrel{D}{⟶}{\chi }^2(p),\text{mbox}asT\to \infty ,$$

$$\frac{F_T}{T}\stackrel{P}{⟶}\Vert C_{XY}-C_0{\Vert }^2,\text{mbox}and\frac{F_{T,p}}{T}\stackrel{P}{⟶}\sum _{\ell =1}^p\frac{⟨C_{XY}-C_0,{\phi }_i{⟩}_2^2}{{\lambda }_i}.$$

$$F_T\stackrel{D}{⟶}\Vert C_A{\Vert }_2^2+2\sum _{i=1}^{\infty }{\lambda }_i^{1/2}⟨C_A,{\phi }_i⟩{\mathcal{N}}_i+\sum _{i=1}^{\infty }{\lambda }_i{\mathcal{N}}_i^2,$$

$$F_{T,p}\stackrel{D}{⟶}\sum _{i=1}^p\frac{⟨C_A,{\phi }_i{⟩}_2^2}{{\lambda }_i}+2\sum _{i=1}^p\frac{⟨C_A,{\phi }_i{⟩}_2{\mathcal{N}}_i}{{\lambda }_i^{1/2}}+\sum _{i=1}^p{\mathcal{N}}_i^2.$$

$$Z_T=T\underset{0\le x\le 1}{\sup }\Vert {\hat{C}}_{XY}(\cdot ,\cdot ,x)-x{\hat{C}}_{XY}(\cdot ,\cdot ,1){\Vert }_2^2.$$

$$Z_{T,p}=T\underset{0\le x\le 1}{\sup }\sum _{i=1}^p\frac{⟨{\hat{C}}_{XY}(\cdot ,\cdot ,x)-x{\hat{C}}_{XY}(\cdot ,\cdot ,1),{\hat{\phi }}_i{⟩}_2^2}{{\hat{\lambda }}_i}.$$

$$Z_T\stackrel{D}{⟶}\underset{0\le x\le 1}{\sup }\sum _{i=1}^{\infty }{\lambda }_i{B}_i^2(x),\text{mbox}asT\to \infty ,$$

$$Z_{T,p}\stackrel{D}{⟶}\underset{0\le x\le 1}{\sup }\sum _{i=1}^p{B}_i^2(x),\text{mbox}asT\to \infty ,$$

$${\hat{D}}_T(t,s,u,v)=\sum _{\ell =-\infty }^{\infty }{W}_b\left(\frac{\ell }{h}\right){\hat{\gamma }}_{\ell }(t,s,u,v),$$

$$h=h(T)\to \infty ,\frac{h}{T^{1/2}}\to 0\text{mbox}asT\to \infty ,$$

$$\displaystyle\hat{\gamma}_{\ell}(t,s,u,v)=\left\{\begin{array}[]{lr}\displaystyle\frac{1}{T}\sum_{j=1}^{T-\ell}\left(\bar{X}_{j}(t)\bar{Y}_{j}(s)-\hat{C}_{XY}(t,s,1)\right)\left(\bar{X}_{j+\ell}(u)\bar{Y}_{j+\ell}(v)-\hat{C}_{XY}(u,v,1)\right),&\ell\geq 0\vspace{.3cm}\\ \displaystyle\frac{1}{T}\sum_{j=1-\ell}^{T}\left(\bar{X}_{j}(t)\bar{Y}_{j}(s)-\hat{C}_{XY}(t,s,1)\right)\left(\bar{X}_{j+\ell}(u)\bar{Y}_{j+\ell}(v)-\hat{C}_{XY}(u,v,1)\right),&\ell<0.\end{array}\right.$$

$${\hat{d}}_T(f)(t,s)=\int \int {\hat{D}}_T(t,s,u,v)f(u,v)dudv,$$

Full text

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 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 Fré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 $L^{2}$ distance or dimension reduction based methods using a suitable principal component basis. These two classes of statistics 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 $\langle\cdot,\cdot\rangle_{d}$ denote the standard inner product on the space $L^{2}[0,1]^{d}$ of real valued square integrable functions defined on $[0,1]^{d}$, and let $\|\cdot\|_{d}=\langle\cdot,\cdot\rangle_{d}^{1/2}$. We write $f$ for the function $f(t)$ when it does not cause confusion. We use the notation $\int$ to denote $\int_{0}^{1}$.

In this section, we suppose that $\{(X_{i}(t),Y_{i}(s)),\;t,s\in[0,1]\}_{i\in{\mathbb{Z}}}$ is a jointly stationary sequence of real valued stochastic processes whose sample paths are in $L^{2}[0,1]$ from which we have observed a sample of length $T$, $\{(X_{1}(t),Y_{1}(s)),...,(X_{T}(t),Y_{T}(s))\}$. For instance, $(X_{i}(t),Y_{i}(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 $\mu_{X}(t)=EX_{0}(t)$, and $\mu_{Y}(s)=EY_{0}(s)$, and define

[TABLE]

to be the cross-covariance function (or kernel) between $\{X_{i}\}$ and $\{Y_{i}\}$ at lag zero. $C_{XY}$ defines the cross-covariance operator $c_{XY}:L^{2}[0,1]\to L^{2}[0,1]$ via

[TABLE]

This relationship implies that we may conduct inference for $c_{XY}$ by conducting inference for the function $C_{XY}$. Based on the sample, we may estimate $C_{XY}$ with

[TABLE]

which denotes the partial sample estimate of $xC_{XY}$ based on the first $x$ proportion of the sample, where

[TABLE]

Under mild regularity conditions on the process $\{(X_{i}(t),Y_{i}(s))\}$, which are implied by the Assumption 2.1 below, $\hat{C}_{XY}(t,s,1)$ is a consistent estimator of $C_{XY}(t,s)$ in $L^{2}[0,1]^{2}$. The motivation for considering now the partial sample estimates of $C_{XY}$ 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 $\ell>0$, one can apply all of the below methods to the sample $\{(X_{\ell+1}(t),Y_{1}(s)),...,(X_{T}(t),Y_{T-\ell}(s))\}$ of length $T-\ell$ 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 $\hat{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 $g_{XY}\colon S^{\infty}\to L^{2}[0,1]\times L^{2}[0,1]$, where $S$ is a measurable space, and a sequence of independent and identically distributed (iid) innovations $\{\epsilon_{i},\colon i\in\mathbb{Z}\}$ taking values in $S$ such that $(X_{i},Y_{i})=g_{XY}(\epsilon_{i},\epsilon_{i-1},\ldots)$.

(b) For all $m\geq 1$, the $m$-dependent sequence $(X_{i,m},Y_{i,m})=g(\epsilon_{i},\ldots,\epsilon_{i-m+1},\epsilon^{*}_{i-m,m},\epsilon^{*}_{i-m-1,m},\ldots)$ with $\epsilon^{*}_{i,m}$ being independent copies of $\epsilon_{i,0}$, and $\{\varepsilon_{i,m}^{*}\colon i\in\mathbb{Z}\}$ independent of $\{\epsilon_{i}\colon i\in\mathbb{Z}\}$, satisfies for some $p>4$,

[TABLE]

where $\alpha,\beta>1$.

Assumption 2.1(a) implies that $\{(X_{i},Y_{i})\}$ is a jointly stationary sequence of Bernoulli shifts in $L^{2}[0,1]\times L^{2}[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 $L^{p}$-$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, $\{\Gamma_{T}(t,s,x),\;t,s,x\in[0,1]\}_{T\in{\mathbb{N}}}$, defined on the same probability space as $\{(X_{i},Y_{i})\},$ that satisfy

[TABLE]

where $E\Gamma_{T}(t,s,x)=0$, and

[TABLE]

where $D(t,s,u,v)$ is the long run covariance function of the sequence $\{(X_{i}(t)-\mu_{X}(t))(Y_{i}(s)-\mu_{Y}(s))\}$, namely

[TABLE]

Theorem 2.1 provides a Skorokhod-Dudley-Wichura type characterization of an invariance principle for $\hat{C}_{XY}$ that can be utilized to establish the asymptotic properties of continuous functionals of $\hat{C}_{XY}$, and we consider several such statistics in Section 3 below in order to carry out hypothesis testing for $C_{XY}$. The proof of Theorem 2.1 is given in Appendix A. The function $D$ describes the asymptotic covariance function of $\sqrt{T}\hat{C}_{XY}(t,s,1)$. $D$ naturally defines a Hilbert-Schmidt integral operator, $d:L^{2}[0,1]^{2}\to L^{2}[0,1]^{2}$, given by

[TABLE]

which further defines an orthonormal basis of eigenfunctions $\varphi_{i}$ in $L^{2}[0,1]^{2}$, and a nonnegative sequence of eigenvalues $\lambda_{1}\geq\lambda_{2}\geq...\geq 0$ satisfying

[TABLE]

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 $C_{XY}$. For instance, we may wish to test based on the estimate $\hat{C}_{XY}(\cdot,\cdot,1)$

[TABLE]

where equality is understood in the $L^{2}[0,1]^{2}$ sense, and $C_{0}$ 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 $X_{i}$ and $Y_{i}$ at a given lag. Since the hypothesis $H_{0,1}$ is well posed only when the sequence $\{(X_{i},Y_{i})\}$ 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:

[TABLE]

[TABLE]

where $C_{XY}^{(i)}(t,s)=\mbox{cov}(X_{i}(t),Y_{i}(s))$. We proceed by developing test statistics for each of these hypotheses. In order to test $H_{0,1}$, we first consider a statistic based on the normalized $L^{2}$ distance of $\hat{C}_{XY}$ to $C_{0}$:

[TABLE]

Corollary 3.1**.**

Under Assumption 2.1 and $H_{0,1}$,

[TABLE]

where $\{\lambda_{i},\;i\geq 1\}$ are defined in (2.2), and $\{\mathcal{N}_{i}$, $i\geq 1\}$ are independent standard normal random variables.

Corollary 3.1 shows that a test of asymptotic size $\alpha$ of $H_{0,1}$ may be obtained by comparing $F_{T}$ to the $1-\alpha$ 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 $\lambda_{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 $\lambda_{i}^{\prime}s$.

The fact that the limiting distribution of $F_{T}$ is nonpivotal though encourages one to consider alternate test statistics based on projecting $\hat{C}_{XY}$ into finite dimensional subspaces of $L^{2}[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

[TABLE]

are asymptotically independent, and hence projecting into the directions of $\varphi_{i}$ has the effect of partitioning the centered estimator $\hat{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

[TABLE]

where $\hat{\varphi}_{i}$, and $\hat{\lambda}_{i}$, $1\leq i\leq p$ are consistent estimates of $\varphi_{i}$ and $\lambda_{i}$, $1\leq i\leq p$, i.e. they satisfy

[TABLE]

where $\hat{c}_{i}=sign(\langle\varphi_{i},\hat{\varphi}_{i}\rangle)$. We discuss in Section 4 how to obtain consistent estimates of $\lambda_{i}$ and $\varphi_{i}(t,s)$.

Corollary 3.2**.**

Under Assumption 2.1 and $H_{0,1}$,

[TABLE]

where $\chi^{2}(p)$ denotes a chi-squared random variable with $p$ degrees of freedom.

We now turn to the consistency and power properties of $F_{T}$ and $F_{T,p}$. The following result shows that both statistics diverge at rate $T$ under $H_{A,1}$, so long as the difference $C_{XY}-C_{0}$ is not orthogonal to the first $p$ elements of the principal component basis $\{\varphi_{i},\;i\geq 1\}$.

Theorem 3.1**.**

Under Assumption 2.1, and $H_{A,1}$,

[TABLE]

Moreover, if $C_{0}=C_{0,T}$ satisfies $\|\sqrt{T}(C_{XY}-C_{0})-C_{A}\|_{2}\to 0$ as $T\to\infty$ for some element $C_{A}$ of $L^{2}[0,1]^{2}$, then

[TABLE]

and

[TABLE]

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 $C_{XY}-C_{0}$ is orthogonal to the first $p$ principal components of $d$, then the test statistic $F_{T,p}$ is not expected to have more than trivial power. Additionally in this case, the norm based test statistic $F_{T}$ is expected to have improved power over other alternatives in which $C_{XY}-C_{0}$ 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 $C_{XY}-C_{0}$ is contained in the subspace spanned by the first $p$ principal components of $d$, then $F_{T,p}$ is expected to be more powerful, since then this statistic effectively defines a uniformly most powerful test of $H_{0,1}$, assuming the data is Gaussian, in the $p$-dimensional subspace spanned by $\varphi_{1},...,\varphi_{p}$.

In order to test $H_{0,2}$ versus $H_{A,2}$, we define analogously to $F_{T}$ and $F_{T,p}$,

[TABLE]

and

[TABLE]

$Z_{T}$ and $Z_{T,p}$ are each maximally selected CUSUM type statistics based on comparing the partial sample estimates of $C_{XY}$ to the estimator from the whole sample. The following corollaries of Theorem 2.1 quantify the large sample behavior of these statistics under $H_{0,2}$.

Corollary 3.3**.**

Under Assumption 2.1 and $H_{0,2}$,

[TABLE]

where $\{\lambda_{i},\;i\geq 1\}$ are defined in (2.2), and $\{B_{i}(x)$, $i\geq 1,\;x\in[0,1]\}$ are iid standard Brownian bridges on $[0,1]$. If in addition (3.2) holds, then

[TABLE]

It follows then that a test of $H_{0,2}$ with asymptotic level $\alpha$ is obtained by rejecting if $Z_{T}$ or $Z_{T,p}$ are larger than the $1-\alpha$ quantiles of their limiting distributions detailed in Corollary 3.3. These limiting distributions may again be approximated using Monte Carlo simulation. The statistics $Z_{T}$ and $Z_{T,p}$ diverge under $H_{A,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 $F_{T,p}$ and $Z_{T,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$, $\varphi_{i}$ and $\lambda_{i}$

As $d$ is simply a scalar multiple of the spectral density operator at frequency zero of the stationary sequence $\{(X_{i}(t)-\mu_{X}(t))(Y_{i}(s)-\mu_{Y}(s))\}$ in $L^{2}[0,1]^{2}$, as defined in [32], it may be naturally estimated with a smoothed periodogram type estimator. Let

[TABLE]

where $h$ is a bandwidth parameter satisfying,

[TABLE]

and, with $\bar{X}_{j}(t)=X_{j}(t)-\bar{X}(t)$ and $\bar{Y}_{j}(s)$ similarly defined,

[TABLE]

We take the function $W_{b}$ to be a symmetric and continuous weight function with bounded support of order $b$; see Chapter 7 of [38]. $\hat{D}_{T}$ then defines an estimator of $d$ by

[TABLE]

which further defines estimates of the eigenvalues and eigenfunctions of $d$ satisfying

[TABLE]

In the simulations and application below, we take the weight function $W_{b}$ to be the simple Bartlett weight function, $W_{1}(x)=(1-|x|)\mathds{1}(|x|<1)$, which is of order 1, and $h=\lceil N^{1/5}\rceil$. Under Assumption 2.1, we have the following result:

Theorem 4.1**.**

Under Assumption 2.1, $\|D-\hat{D}_{T}\|_{4}=o_{P}(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\geq 1$ satisfying that

[TABLE]

where $\{\lambda_{i},\;i\geq 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 $X_{i}(t_{j})$ and $Y_{i}(t_{j})$, $1\leq j\leq R$, denote the observed values of the functions $X_{i}(t)$ and $Y_{i}(t)$, observed at the common points $\{t_{1},\ldots,t_{R}\}$. We assume here that each functional observation $X_{i}(t)$ and $Y_{i}(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 $F_{T}$, $F_{T,p}$, $Z_{T}$, and $Z_{T,p}$ from the discrete data and simple Riemann sum approximations to the inner products and norms, so long as the eigenvalue and eigenfunction estimates $\hat{\lambda}_{i}$ and $\hat{\varphi}_{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 $(\lambda_{i},\varphi_{i}(t,s)),\;\;1\leq i\leq p$ of $d$ by calculating the spectrum of the 4-way tensor $\hat{D}_{T}(t_{i},t_{j},t_{k},t_{\ell})$, $1\leq i,j,k,\ell\leq R$ of dimension $R^{4}$, 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 $\{X_{i}\}$ and $\{Y_{i}\}$. Let

[TABLE]

be the sample covariance matrix of the discretized observations of the $X$ sample, and define $\widehat{c}_{Y}$ similarly. Calculating the spectrum of $\widehat{c}_{X}$ and $\widehat{c}_{Y}$ yields eigenvalues $\{\nu_{X,1},\ldots,\nu_{X,R}\}$, $\{\nu_{Y,1},\ldots,\nu_{Y,R}\}$ and eigenvectors, $\{\theta_{X,1},\ldots,\theta_{X,R}\}$, $\{\theta_{Y,1},\ldots,\theta_{Y,R}\}$, the latter of which, when multiplied by $\sqrt{R}$, yield discrete approximations to the fPC’s of the sequences $\{X_{i}\}$ and $\{Y_{i}\}$, $\theta_{X,i}(t_{j})$, $\theta_{Y,i}(t_{j})$, $1\leq i,j\leq R$, respectively.

We then use the product basis of $L^{2}[0,1]^{4}$ generated by these functions to reduce the dimension of $\hat{D}_{T}$. The projections of $\hat{D}_{T}$ onto the product basis of the first $q$ elements each of $\theta_{X,i}(t)$, $\theta_{Y,i}(t)$ may be stored in a 4-way tensor with $q^{4}$ elements, $\mathcal{M}$, via

[TABLE]

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 $\theta_{X,i}(t)$ and $\theta_{Y,i}(t)$, $1\leq i\leq q$, explains at least $v\%$ of the variation in each series, where $v$ is close to but strictly smaller than 1.

The eigenvalues $\{\hat{\lambda}_{i},1\leq i\leq q^{2}\}$ and eigen-arrays $\{\hat{{\bf\Phi}}_{i}\in\mathbb{R}^{q\times q},1\leq i\leq q^{2}\}$ of $\mathcal{M}$ satisfy

[TABLE]

.

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 $\mathcal{M}$, as is implemented in svd.tensor function in the the package tensorA in R; see [45]. The eigenfunctions $\varphi_{i}(t,s)$ may then be estimated with

[TABLE]

We use these estimates for the $\hat{\lambda}_{i}^{\prime}s$ and $\hat{\varphi}_{i}^{\prime}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 $F_{T}$, $F_{T,p}$, $Z_{T}$, and $Z_{T,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 $\alpha\in[0,1]$:

[TABLE]

where $\varepsilon_{c,i}(t)$, $\varepsilon_{x,i}(t)$, and $\varepsilon_{y,i}(t)$ are mutually independent sequences that satisfy either of two models given below. In this case, the functional series $X_{i}(t)$ and $Y_{i}(t)$ are correlated through their common dependence on $\varepsilon_{c,i}(t)$, with the strength of this dependence controlled by $\alpha$. If $\alpha=0$ for instance, then the two sequences $X_{i}$ and $Y_{i}$ are independent. We took the error sequences to satisfy either:

• IID: $\varepsilon_{c,i}(t)=W_{c,i}(t)$, $\varepsilon_{x,i}(t)=W_{x,i}(t)$, and $\varepsilon_{y,i}(t)=W_{y,i}(t)$, where $\{W_{c,i}(t)\}$, $\{W_{x,i}(t)\}$, $\{W_{y,i}(t)\}$ are mutually independent sequences of IID standard Brownian motions.
• FAR(1): $\varepsilon_{c,i}(t)=\int\Phi(t,s)\varepsilon_{c,i}(s)ds+W_{c,i}(t)$, $\varepsilon_{x,i}(t)=\int\Phi(t,s)\varepsilon_{x,i}(s)ds+W_{x,i}(t)$, and $\varepsilon_{y,i}(t)=\int\Phi(t,s)\varepsilon_{y,i}(s)ds+W_{y,i}(t)$, where $\{W_{c,i}(t)\}$, $\{W_{x,i}(t)\}$, and $\{W_{y,i}(t)\}$ are defined above, and $\Phi(t,s)=\min(t,s)$.

The sequence $\{(X_{i},Y_{i})\}$ 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 $\alpha$, 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 $X_{i}$ and $Y_{i}$. 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 $F_{T,p}$ and $Z_{T,p}$.

To study the size of the test of $H_{0,1}$ based on $F_{T}$ and $F_{T,3}$, we generated data according to (5.1) for each setting of the error sequence with $\alpha=0$, and performed tests with nominal levels of 10%, 5%, and 1% of $H_{0,1}:C_{XY}=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\geq 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 $F_{T}$ and $F_{T,p}$ by testing $H_{0,1}:C_{XY}=0$ with data following (5.1) and $\alpha$ increasing from [math] 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 $\alpha$, although it is difficult to quantify their power without a comparable test. Interestingly, the simple norm based test based on $F_{T}$ possessed better power than $F_{T,3}$ in all the examples that we considered.

Regarding the test of $H_{0,2}$ based on $Z_{T}$ and $Z_{T,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 $\alpha=0$ and $\alpha=0.5$ for $T=100,300$ and $1000$. The empirical size from 1000 independent simulations for each setting is presented in Table 5.2. We saw that in general the change point tests based on $Z_{T}$ and $Z_{T,3}$ had quite good size, even in the presence of serial correlation. Changing $\alpha$ had no effect on the size of these tests, as expected.

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 $P_{M,i}(t_{j}),$ and $P_{X,i}(t_{j})$ $i=1,\ldots,T,j=1,\ldots,R$, denote the prices of Microsoft and Exxon mobile stock on day $i$ at intraday time $t_{j}$, 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:

Definition 1**.**

Suppose $P_{i}(t_{j}),i=1,\ldots,T,j=1,\ldots,R$, is the price of a financial asset at time $t_{j}$ on day $i$. The functions

[TABLE]

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 $X_{i}(t)=r_{M,i}(t)$ and $Y_{i}(t)=r_{X,i}(t)$, $1\leq i\leq 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 $H_{0,2}$, using the test statistic $Z_{T}$. 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\|\hat{C}_{XY}(\cdot,\cdot,t/T)-(t/T)\hat{C}_{XY}(\cdot,\cdot,1)\|_{2}^{2}$ 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 $Z_{T}$ 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 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 $T_{1}=56$ and $T_{2}=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 $H_{0,1}:C_{XY}=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.

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 $\mu_{X}=\mu_{Y}=0$. Let

[TABLE]

Lemma A.1**.**

Under Assumption 2.1,

[TABLE]

Proof.

According to the definitions of $\hat{C}_{XY}$ and $\tilde{C}_{XY}$ and the triangle inequality,

[TABLE]

For the term $G_{1}(x)$, we have that

[TABLE]

Assumption 2.1 implies that both the series $X_{i}$ and $Y_{i}$ satisfy the conditions of Theorem 3.3 of Berkes et al [7], from which it follows that

[TABLE]

This combined with (A.5) implies that $\sup_{0\leq x\leq 1}G_{1}(x)=O_{P}(1/T)$. Parrel arguments show that $\sup_{0\leq x\leq 1}G_{2}(x)=O_{P}(1/T)$ and $\sup_{0\leq x\leq 1}G_{3}(x)=O_{P}(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

[TABLE]

where $\gamma=\min(\alpha,\beta)$, and $X_{i,m}$ and $Y_{i,m}$ are defined in Assumption 2.1.

Proof.

We have according to the triangle inequalities in $L^{2}[0,1]^{2}$ and for $(E(\cdot)^{q})^{1/q}$ that

[TABLE]

According to the Cauchy-Schwarz inequality and stationarity, we have that

[TABLE]

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$)

[TABLE]

Since both $(E\|X_{0}\|_{1}^{p})^{1/p}$ and $(E\|Y_{0}\|_{1}^{p})^{1/p}$ are finite, and according to Assumption 2.1 (E\|X_{i}-X_{i,m}\|^{p}\big{)}^{1/p}=O(m^{-\alpha}),\mbox{ and }(E\|Y_{i}-Y_{i,m}\|^{p}\big{)}^{1/p}=O(m^{-\beta}), we obtain that

[TABLE]

where $\gamma=\min(\alpha,\beta)$ as needed.

∎

Proof of Theorem 2.1.

According to Lemma A.1, it is enough to prove Theorem 2.1 for the process $\tilde{C}_{XY}(t,s,x)$. $\tilde{C}_{XY}(t,s,x)-(\lfloor Tx\rfloor/T)C_{XY}(t,s)$ is the partial sum of the random functions $\xi_{i}(t,s)=X_{i}(t)Y_{i}(s)-EX_{0}(t)Y_{0}(s)$ which form a mean zero and stationary sequence in $L^{2}[0,1]^{2}$. Let $\xi_{i,m}(t,s)=X_{i,m}(t)Y_{i,m}(s)-EX_{0}(t)Y_{0}(s)$, denote an $m$-dependent approximation to $\xi_{i}$. It follows directly from A.2 that

[TABLE]

where $\gamma>1$, and this implies that

[TABLE]

It now follows from Theorem 1 of [27] that if

[TABLE]

then

[TABLE]

where $\Gamma_{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$,

[TABLE]

where

[TABLE]

with $\{W_{\ell}(t),\;t\in[0,1],\;\ell\in\mathbb{N}\}$ being iid standard Brownian motions, and $(\lambda_{\ell}$ $\varphi_{\ell})$ being the eigenvalues and eigenfunctions of the operator $d$. This follows from simply calculating the mean and covariance function of $\Gamma_{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 result, we have using the assumption that $\sqrt{T}(C_{XY}-C_{0})\stackrel{{\scriptstyle L^{2}[0,1]^{2}}}{{\to}}C_{A}$ and Theorem 2.1 that, in case of $F_{T}$,

[TABLE]

as needed. The result for $F_{T,p}$ follows similarly. ∎

A.2. Proof of Theorem 4.1

Let

[TABLE]

with $h$ defined in (4.2), $W_{b}$ being a bounded, symmetric, and continuous weight function of order $b$, and

[TABLE]

Lemma A.3**.**

Under Assumption 2.1,

[TABLE]

Proof.

According to the definition of $\hat{\gamma}_{\ell}$, when $\ell\geq 0$, we have that

[TABLE]

where the terms $R_{p,\ell}$ represent the remaining 24 terms in the definition of $\hat{\gamma}_{\ell}$ obtained by completing the multiplication that results in 25 terms, the first of which corresponds to $\tilde{\gamma}_{\ell}$. We now aim to show that for all $p$ and $\ell$, $E\|R_{p,\ell}\|_{4}=O(T^{-1/2})$. For $R_{1,\ell}$, we have by the triangle inequality, a couple of applications of the Cauchy-Schwarz inequality, and the assumed stationarity that

[TABLE]

It follows from Proposition 3.1 of [44] that under Assumption 2.1, $(E\|\bar{X}\|_{1}^{4})^{1/4}=O(T^{-1/2})$, and also according to Assumption 2.1, $(E\|Y_{0}\|_{1}^{4})^{1/4}=O(1)$ and $(E\|(X_{0}(t)Y_{0}(s)-EX_{0}(t)Y_{0}(s))\|_{2}^{2})^{1/2}=O(1)$. This implies with (A.11) that $E\|R_{1,\ell}\|_{4}=O(T^{-1/2})$. Similar arguments using the fact that $(E\|\bar{Y}\|_{1}^{4})^{1/4}=O(T^{-1/2})$ and $\|EX_{0}(t)Y_{0}(s)-\hat{C}_{XY}(t,s,1)\|_{2}=O(T^{-1/2})$ can be applied to get that $E\|R_{p,\ell}\|=O(T^{-1/2})$, $2\leq p\leq 24$. This implies via the triangle inequality and (A.10) that

[TABLE]

from which it follows that

[TABLE]

by (4.2) and the fact that $W_{b}$ has bounded support. A parallel argument may be used to show that for $\ell<0$, that $E\|\hat{\gamma}_{\ell}-\tilde{\gamma}_{\ell}\|_{4}=O(T^{-1/2}),$ from which we obtain that

[TABLE]

(A.12) and (A.13) together with the triangle inequality imply that

[TABLE]

and this implies the result with Markov’s inequality. ∎

Lemma A.4**.**

Under Assumption 2.1,

[TABLE]

Proof.

$\tilde{\gamma}_{\ell}$ is the autocovariance estimator based on a sample of size $T$ at lag $\ell$ of the mean zero random functions $\{\xi_{i}(t,s),\;i\in{\mathbb{Z}}\}$, where, recalling from above, $\xi_{i}(t,s)=X_{i}(t)Y_{i}(s)-EX_{0}(t)Y_{0}(s)$. With $\xi_{i,m}(t,s)=X_{i,m}(t)Y_{i,m}(s)-EX_{0}(t)Y_{0}(s)$, we have from A.2 and Lyapounov’s inequality that

[TABLE]

and

[TABLE]

This shows that $\xi_{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 $Z_{T}$ and $Z_{T,p}$

In order for the tests based on $Z_{T}$ and $Z_{T,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 $H_{A,2}$, we assume that

[TABLE]

and

[TABLE]

Theorem A.1**.**

Under $H_{A,2}$ and Assumption A.1 holds, then

[TABLE]

If in addition $C_{1}$ and $C_{2}$ in $H_{A,2}$ satisfy $\langle C_{1}-C_{2},\varphi_{i}\rangle_{2}\neq 0$ for some $1\leq i\leq p$, then

[TABLE]

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 $g_{XY}$ which defines how the series $(X_{i},Y_{i})$ depends on the underlying innovation sequence were to change from $g_{XY}^{(1)}$ to $g_{XY}^{(2)}$ in such a way that $H_{A,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.