On changepoint detection in functional data using empirical energy distance

On changepoint detection in functional data using empirical energy distance

\email{[email removed]} \email{[email removed]} \email{[email removed]} \address{$^{1}$Department of Mathematics, Drexel University, Philadelphia, PA 19104 USA } \address{$^2$Department of Mathematics, University of Utah, Salt Lake City, UT 84112--0090 USA } \address{$^3$Department of Economics, Finance and Accounting; School of Business and Economics, University of Leicester, Leicester, U.K.; Department of Economics and Management, University of Pavia, Pavia, Italy.} \address{$^*$Research supported in part under NSF grant DMS-2309570.} \subjclass[2020]{60F17}

\doublespacing

Introduction

The analysis of datasets where the data are observed as functions, rather than scalars or vectors, has been investigated in numerous contributions over the past few years. Functional Data Analysis (FDA) has become ubiquitous in virtually all applied sciences, in the case where data are genuinely functional in nature, and also when a parsimonious description of the data is called for. FDA appears naturally in the analysis of economic and financial data; examples include analysing the term structure of interest rates, where, for each time period, the observed maturities are the discrete approximation of the continuum of maturities (hays); and modelling intraday return density trajectories (bathia2010). In climate science, it is typical to model temperatures - which are recorded at a high frequency basis, e.g. several times per day - at a lower frequency (e.g. yearly), with the intra-period data representing the discretised functional observations - see, for example, horvath:kokoszka:2012 and king2018functional. In medical imaging, several datasets arise that can be modelled as possibly multi-dimensional and functional valued ( sangalli). See also ramsay2002applied for further examples.

On account of the huge relevance of the topic, inferential theory for FDA has been studied in many contributions. The literature has developed useful dimension reduction tools such as the functional version of Principal Components (hall2006properties), and the full-blown estimation theory for linear regression, dynamic models and also nonlinear models such as the functional version of ARCH and GARCH (horvath:kokoszka:2012 ). However, the validity of inferential theory often hinges on having some stability in the structure of the data, such as the constancy of the mean function, or of the whole distribution. Hence, testing for the possible presence of changepoints (in the mean, in higher order moments, or even in the whole distribution) is of paramount importance.

Changepoint detection is well-studied in the context of scalar or vector-valued time series, and we refer, inter alia, to casini2019 for a useful review containing several examples and applications. In contrast, changepoint analysis in functional data has received only limited attention. berkes2009detecting propose a CUSUM-based test statistic to detect changepoints in the mean of independent functional-valued observations; hormann:kokoszka:2010, zhang2011testing, aston and aue2018detecting consider extensions to deal with dependent data, which typically occur in a time series context. These contributions, broadly speaking, are based on the unweighted CUSUM process, and it is possible to show that, in this case, $N^{1/2}$ periods away from the beginning/end of the sample (where $N$ is the sample size). On the other hand, detection of early/late occurring breaks is very important, due to its implications on the ability to assess timely whether a model which has been valid so far is still appropriate e.g. for forecasting.

Main contributions of this paper

In this paper, we bridge the gaps discussed above, by proposing a procedure to detect changepoints for serially dependent, possibly multivariate functional-valued time series, allowing for breaks to occur close to the sample endpoints. We consider very general changes, which could occur in various functions/functionals of the data, including the mean, higher order moments, and in general functions which completely characterize the underlying distribution such as the characteristic function. Specifically, we develop a novel family of weighted test statistics based on the notion of energy distance (see e.g. szekely:rizzo:2005; szekely:rizzo:2017; and baringhaus:franz:2004). The energy distance is a metric designed to measure the distance between the distributions of two independent random vectors (say $X$ and $Y$), defined as

where $X^{\prime },Y^{\prime }$ are independent copies of $X$ and $Y$, respectively, $\left\vert \cdot \right\vert $ is the Euclidean norm, and $ 0<\eta <2$; it can be shown (see Theorem 2 in szekely:rizzo:2005) that $\mathcal{E}_{\eta }(X,Y)=0$ if and only if $X$ and $Y$ have the same distribution. Empirical energy distances have been used by matteson:james:2014 and biau:bleakley:mason:2016 to study distributional changepoint problems for a sequence of independent, vector-valued time series; chakraborty2021high extend the theory to the case of high-dimensional sequences.

Taking $\eta =2$ in ((ref)) leads to a statistic suitable for testing equality of expectations rather than equality of distributions ( szekely:rizzo:2005). Hence, we consider an empirical version of $ \mathcal{E}_{2}(X,Y)$, constructed at every point in the sample $1\leq k\leq N$, comparing the sample average before and after $k$ in a (conceptually) similar way to the CUSUM\ process. We then consider weighted versions of the empirical energy process, with weights designed to boost the value taken by the process when $k$ is close to the beginning/end of the sample. Changepoint detection can thus be based on the maximally selected weighted empirical energy process. The resulting tests have nontrivial power versus breaks occurring (much) closer to the sample endpoints than $N^{1/2}$ periods, while still having power versus mid-sample breaks. In Section (ref) we show that the limiting distribution of our test statistics contains the integral of the square of a Gaussian process which depends - in a highly nontrivial way - on nuisance parameters. Hence, in order to compute critical values, we propose a method based on the Karhunen-Lo\`{e}ve (KL henceforth) expansion, which appears to be easier to use than e.g. the bootstrap (see e.g., albeit in a different context, inoue2001). Our theory is stated for the general case of multivariate functional time series whose argument can also be multivariate, which is relevant in several applications of FDA, including shape analysis (kenobi2010shape) and medical imaging (kurtek2010novel).

For the sake of clarity, our presentation focuses mainly on detecting changes in the mean of functional observations. However, our approach can be readily applied to consider different changepoint problems; in Section (ref) we discuss how our tests can be used to detect distributional changes, by applying it to the empirical characteristic function. Testing for changes in the distribution is arguably of great importance; as inoue2001 puts it, \textquotedblleft \lbrack ...]\ stability of distribution, moments, or parameters is essential to the proofs of asymptotic properties of the maximum likelihood method, generalized method of moments, and nonparametric method. Consequently, instability can affect estimation and inference.\textquotedblright\ (p. 156). Contributions on this topic often require independence assumptions, and are relatively scarce even in the case of scalar or vector-valued observations: in addition to the papers by matteson:james:2014 and others referred to above, other approaches include inoue2001, who uses the unweighted CUSUM process based on the empirical distribution function; antoch2008, who use a combination of rank statistics; and huskova:meintanis:2006, who use the empirical characteristic function for scalar observations.

The remainder of the paper is organised as follows. We present our test statistics in Section (ref). We study its asymptotic theory in Section (ref): we derive the weak limit under the null in Section (ref); we study power, estimation of the breakdate, and binary segmentation in Section (ref); we offer a methodology to compute critical values in Section (ref). We extend our approach to detecting changes in the distribution of the data is in Section (ref) . In Section (ref), we report a comprehensive simulation exercise; an empirical application to intraday returns is in Section (ref). Section (ref) concludes. Further Monte Carlo evidence, an empirical application to temperature data, lemmas and proofs are relegated to the Supplement.

NOTATION. Henceforth, $\mathcal{T}$ denotes a compact subset of $\mathbb{R} ^{d}$, and $\{x(t),~t\in \mathcal{T}\}$ is a square integrable function; whenever convenient, we write $x$ in place of $x(t)$ or in place of $ \{x(t),~t\in \mathcal{T}\}$. For any $r\geq 1$, given two square integrable $ \mathbb{R}^{r}$-valued functions $\{x(t),~t\in \mathcal{T}\}$ and $ \{y(t),~t\in \mathcal{T}\}$, we define the inner product $\langle x,y\rangle =\int_{\mathcal{T}}x^{\top }(t)y(t)dt$, where \textquotedblleft $^{\top }$ \textquotedblright\ is the usual transpose; and we define the $L^{2}$-norm $ \Vert x\Vert =\sqrt{\langle x,x\rangle }$, writing $x=y$ if $\Vert x-y\Vert =0$. When unambiguous, we write $\{a_{\ell }\}$ to denote a given sequence $ \{a_{\ell },-\infty <\ell <\infty \}$. We also write the symbol $\int $ in place of $\int_{\mathcal{T}}$. We use: \textquotedblleft $\underset{\mathcal{ D}[0,1]}{\overset{w}{\longrightarrow }}$\textquotedblright\ to denote weak convergence in $\mathcal{D}[0,1]$; \textquotedblleft $\overset{{\mathcal{D}}} {\rightarrow }$\textquotedblright\ to denote convergence in distribution; \textquotedblleft $\overset{\mathcal{P}}{\rightarrow }$\textquotedblright\ for convergence in probability; \textquotedblleft a.s.\textquotedblright\ for \textquotedblleft almost surely\textquotedblright ; \textquotedblleft $ \overset{{\mathcal{D}}}{=}$\textquotedblright\ for equality in distribution; $\lfloor \cdot \rfloor $ for the integer value function; and $\left\vert \cdot \right\vert $ to denote the Euclidean norm of a vector, or the Frobenius norm of a matrix. Other relevant notation is introduced further in the paper.

The test statistics: definition, assumptions and asymptotics

We consider a sequence of $\mathbb{R}^{r}$-valued functional observations of the form

where for each $i$, $\mu _{i}$ and $\epsilon _{i}$ are $\mathbb{R}^{r}$ -valued square integrable functions. We aim to test

against the $R$-change alternative:

for $i=1,\ldots ,R$, with the convention that $k_{0}=1$ and $k_{R+1}=N$.

As discussed in the introduction, a possible way of detecting changes is based on the energy distance defined in ((ref)). Our approach is based on weighted functionals of $\{V_{N}(k),\,2\leq k\leq N-2\}$, defined as the empirical version of the energy distance\footnote{ See also sejdinovic:sriperumbudur:gretton:fukumizu:2013 for further discussion on generalizations of the energy distance.} calculated for $\eta =2$

Throughout the paper, we assume that the sequence $\{X_{i}\}$ is weakly dependent:

Assumption (ref) states that the sequence $\{\epsilon _{\ell },~-\infty <\ell <\infty \}$ is stationary and ergodic, and it can be approximated by a sequence with finite-order dependence (see hormann:kokoszka:2010). The Bernoulli shift representation in Assumption (ref) is widely employed in the analysis of scalar time series, where it can be verified in the most commonly used DGPs in econometrics and statistics. As far as functional-valued data are concerned, examples when Assumption (ref) hold include linear processes in Hilbert spaces ( horvath:kokoszka:2012), and a large class of non-linear processes, including functional ARCH and GARCH models ( aue:horvath:pellatt:2017) and bilinear models ( hormann:kokoszka:2010).

Asymptotics

We will consider the following statistics:

where $0\leq \alpha <1$ and $u=k/N$. A full-blown discussion is after Theorems (ref) and (ref) below; here we offer a heuristic preview of the rationale of ((ref)). The statistic $V_{N}\left( \left\lfloor Nu\right\rfloor \right) $ could be sensitive to outliers occurring a few periods after the start of the sample, which could inflate $ V_{N}\left( \left\lfloor Nu\right\rfloor \right) $ and lead to a spurious rejection of the null of no changepoint; the weights $\left( u\left( 1-u\right) \right) ^{2}$ reduce the impact of outliers close to the sample endpoints on $V_{N}\left( \left\lfloor Nu\right\rfloor \right) $. On the other hand, this weighing scheme also reduces power in the presence of a genuine break at the beginning/end of the sample. The further weight $\left( u\left( 1-u\right) \right) ^{-\alpha }$ is designed to \textquotedblleft pick up\textquotedblright\ the test statistic at sample endpoints, boosting power versus breaks located close to the sample endpoints. The process in ( (ref)) can be compared with the weighted CUSUM process employed in the changepoint detection literature ( csorgo1997); as we show in Section (ref), larger values of $\alpha $ result in having nontrivial power versus changepoints closer to the beginning/end of the sample.

Asymptotics under the null

Let

Define also the process

where $\left\{ \Gamma (u,t),0\leq u\leq 1,t\in \mathcal{T}\right\} $ is an $ r $-dimensional Gaussian process with $E\Gamma (u,t)=0$ and covariance kernel $E\left( \Gamma (u,t)\Gamma ^{\top }\left( u^{\prime },t^{\prime }\right) \right) $ $=$ $\left( \min \left\{ u,u^{\prime }\right\} -uu^{\prime }\right) \mathbf{D}(t,t^{\prime })$. \newline Our first main result provides the functional weak limit of weighted versions of $V_{N}$ under $H_{0}$.

Theorem (ref) is the building block to carry out changepoint detection. We note that, heuristically, the Gaussian process $\Gamma (u,t)$ is a Brownian bridge at each \textquotedblleft slice\textquotedblright\ across $t$ . As we show in Lemma (ref), this is a consequence of the fact that, under the null, $V_{N}\left( k\right) $ and its weighted versions are well approximated by the (weighted) squared CUSUM process, modulo some extra terms that, in the limit, either vanish or enter the expression as constants.

A natural approach to test for changepoints is to use the max-type statistic

Under $H_{0}$, it follows by Theorem (ref) and continuity\footnote{ The Law of the Iterated Logartihm for Gaussian processes entails that the limit in ((ref)) is a.s. finite - see also the proof of Theorem (ref) for details. Indeed, having $\alpha <1$ is crucial to this argument, since it also holds that

} that

From a practical point of view, the limiting law of $T_{N}$ contains several nuisance parameters, such as the covariance kernel $\mathbf{D}(t,t^{\prime }) $ defined in ((ref)), and the variance $\sigma _{0}^{2}$ defined in ((ref)). In Section (ref), we discuss the computation of critical values for tests based on $T_{N}$. From a technical point of view, one of the main ingredients to show Theorem (ref) is the weak invariance principle for partial sums of dependent functional time series, shown in berkes:horvath:rice:2013. However, in our case we consider a weighted version of the partial sum process, which requires a nontrivial extension of the arguments in berkes:horvath:rice:2013. The case $\alpha =1$ is also of interest, and it corresponds to the standardised CUSUM\ (csorgo1997); studying this would require a strong invariance principle for partial sums of functional time series which, to our knowledge, is not available in the literature.

Asymptotics under the alternative

Consistency under a single break and asymptotic power function

We begin by considering the case of a single break - i.e., $R=1$ in ((ref)) - in the presence of a changepoint of size

and we also defined its rescaled counterpart as

Define

$a_{N}=N^{1/2}\Vert \mathcal{\delta }\Vert /\left( 2\sigma (\theta (1-\theta ))^{\frac{3}{2}-\alpha }\right) $, and let $\mathcal{N}$ denote a standard normal random variable.

Theorem (ref) states that tests based on $T_{N}$ have power in the presence of changepoints of possibly vanishing magnitude - i.e. $\left\Vert \mathcal{\delta }\right\Vert =o\left( 1\right) $ - and occurring close to sample endpoints - i.e. $k^{\ast }=o\left( N\right) $ or $N-k^{\ast }=o\left( N\right) $. In order to understand the result in Theorem (ref) , some examples may be helpful. Considering the case of a mid-sample break - with $k^{\ast }=cN$ for some $0<c<1$ - equation ((ref)) boils down to requiring $\lim_{N\rightarrow \infty }N\left\Vert \mathcal{\delta } \right\Vert ^{2}=\infty $. Hence, mid-sample breaks can be detected even when the magnitude $\left\Vert \mathcal{\delta }\right\Vert $ drifts to zero as $N\rightarrow \infty $, as long as $\left\Vert \mathcal{\delta } \right\Vert $ shrinks at a rate slower than $N^{-1/2}$. Conversely, consider the case of a non-vanishing break, i.e. $\left\Vert \mathcal{\delta } \right\Vert >0$, and a changepoint located close to the beginning of the sample, viz. $k^{\ast }=o\left( N\right) $. In such a case, changepoints can be detected as long as they occur at least $N^{\left( 1-\alpha \right) /\left( 2-\alpha \right) }$ periods from the beginning of the sample. In the unweighted case $\alpha =0$, this reflects that breaks occurring $o\left( N^{1/2}\right) $\ periods from the sample endpoints cannot be reliably detected; on the other hand, increasing $\alpha $ makes tests more able to detect breaks occurring closer to the beginning/end of sample. We note however that, upon inspecting our proofs, the rates of asymptotic approximation deteriorate as $\alpha $ approaches $1$, thus reflecting the size/power trade-off. Finally, equation ((ref)) describes the asymptotic power function in the case of a changepoint occurring \textquotedblleft not too close\textquotedblright\ to the sample endpoints.

Estimation of the breakdate: consistency and limiting distribution

We now consider, in greater depth, the case where the changepoint occurs mid-sample, viz.

for some $0<\theta <1$. The max-type statistic $T_{N}$ defined in ((ref)) gives the estimator

from which the estimated breakdate can be computed as $\widehat{k} _{N}=\left\lfloor N\widehat{\theta }_{N}\right\rfloor $. In the next theorem, we state the consistency of the break fraction estimator $\widehat{ \theta }_{N}$, and derive its asymptotic distribution in the (customarily studied) case where the size of the break drifts to zero as $N\rightarrow \infty $.\footnote{ The fixed break case, i.e. $\left\Vert \delta \right\Vert >0$, can be studied along similar lines as the proof of Theorem (ref); however, in this case the limiting distribution of the estimated changepoint depends on many nuisance parameters, thus being of scarce practical use. In the large break case $\Vert \delta \Vert \rightarrow \infty $, it can be shown that $P( \widehat{k}_{N}=k_{*})\rightarrow 1$.} We define the drift function

with $m_{\alpha }\left( 0\right) =0$, where $I\left( \cdot \right) $ is the indicator function; and the two-sided standard Wiener process $\widetilde{W} \left( u\right) =W_{1}\left( -u\right) I\left( u\leq 0\right) +W_{2}\left( u\right) I\left( u\geq 0\right) $, where $\left\{ W_{1}\left( u\right) ,u\geq 0\right\} $ and $\left\{ W_{2}\left( u\right) ,u\geq 0\right\} $\ are two independent standard Wiener processes.

According to Theorem (ref), the estimator of the break fraction $\theta $ is consistent; the estimated breakdate $\widehat{k}_{N}$ is also consistent in the sense that $\widehat{k}_{N}-{k^{\ast }}=o_{P}\left( N\right) $. Theorem (ref) refines the consistency of $\widehat{\theta } _{N}$ in the case of a break of vanishing magnitude, stating, in essence, that $\widehat{k}_{N}-{k^{\ast }}=O_{P}\left( \Vert \mathcal{\delta }\Vert ^{-2}\right) $. The limiting distribution is the same as one would have when using the maximally selected weighted CUSUM process - this (again) reinforces the conclusion from Lemma (ref) that $ V_{N}\left( k\right) $ is related to the CUSUM process. In principle, it would be possible to construct confidence intervals for ${k^{\ast }}$, by simulating the percentiles of the (nuisance free) random variable $\xi _{\alpha }$ calculated at $\widehat{\theta }_{N}$, and using the means $ \widehat{k}_{N}^{-1}\sum_{i=1}^{\widehat{k}_{N}}X_{i}\left( t\right) $ and $ \left( N-\widehat{k}_{N}\right) ^{-1}\sum_{i=\widehat{k}_{N}+1}^{N}X_{i} \left( t\right) $ to estimate $\mathcal{\delta }\left( t\right) $ and the long run variance $\sigma ^{2}=\int \int^{\top }\Vert \mathcal{\delta }\Vert ^{-2}\mathcal{\delta }^{\top }\left( t\right) \mathbf{D}(t,t^{\prime }) \mathcal{\delta }\left( t^{\prime }\right) dtdt^{\prime }$.

The case of multiple breaks: binary segmentation

We now consider the case of multiple breaks. Recalling that $I\left( \cdot \right) $ is the indicator function, this case corresponds to

We consider the case of \textquotedblleft well-separated\textquotedblright\ breaks of non-vanishing magnitude.

In this case, it is possible to show that our tests have power, by marginally adapting the proof of Theorem (ref). Here, we discuss in greater detail how to estimate the number of changepoints, $R$, in addition to the locations thereof. Whilst the literature has developed several techniques, we focus on the binary segmentation approach proposed by vostrikova1982detection. The algorithm can be described as follows (see also Algorithm (ref) in the Supplement for pseudocode). Starting from the whole sample, we apply our test using a fixed $0<\alpha <1$, to check whether there is at least one changepoint. If a break is detected, we estimate its location using ((ref)), and then split the sample around the estimated breakdate. The procedure is then iterated on each subsample, until either no changepoint is detected, or a stopping rule (typically based on the length of the sub-sample) is triggered.

Formally, consider a subsample with starting and ending points $1\leq \ell <u\leq N$, under the constraint that $u-\ell > 4$; define the weighted statistic $\left( \left( k-\ell \right) \left( u-k\right) /\left( u-\ell \right) ^{2}\right) ^{-\alpha }V_{N}^{\left( \ell ,u\right) }\left( k\right) $, where

and let its maximally selected counterpart be $T_{N}^{\left( \ell ,u\right) } $. The interval $(\ell ,u)$ is marked to have a changepoint if $ T_{N}^{\left( \ell ,u\right) }$ exceeds a (user-chosen) threshold $\tau _{N}$ . Practically, the choice of the threshold $\tau _{N}$ can be based on any slowly diverging sequence satisfying mild growth constraints (see expression (ref), below), and we refer to Section (ref) in the Supplement for examples. Hence, the corresponding changepoint estimator in the interval $(\ell ,u)$ can be defined as $\widehat{k}$ $=$ $ \operatorname*{sargmax}_{l\leq k\leq u}\left( \left( k-\ell \right) \left( u-k\right) /\left( u-\ell \right) ^{2}\right) ^{-\alpha }V_{N}^{\left( \ell ,u\right) }\left( k\right) $, where \textquotedblleft $\operatorname*{sargmax}$\textquotedblright\ denotes the smallest integer that maximizes the expression. The sample is then split around $\widehat{k}$, and the procedure iterated until it comes to a stop. The final output is a set of estimated changepoints $\widehat{ \mathcal{H}}=\left\{ \widehat{k}_{1},...,\widehat{k}_{\widehat{R}}\right\} $ sorted in increasing order, and the estimate $\widehat{R}$.

Theorem (ref) stipulates the consistency of $\widehat{R}$ and of the breaks locations, $\left\{ \widehat{k}_{r}\right\} _{r=1}^{R}$. Heuristically, this is because, under the alternative, $V_{N}\left( k\right) $ is equal to the (squared) CUSUM process plus a \textquotedblleft small\textquotedblright\ term, thus having the same properties as the CUSUM. Importantly, our results require $\alpha >0$, which reinforces the importance of considering weighted statistics.

Computation of critical values

By the multivariate KL expansion (happ), the $\mathbb{R}^{r}$ -valued Gaussian process $\Gamma (u,t)$ in (ref) can represented as

where $\left\{ B_{\ell }(u),0\leq u\leq 1\right\} $ is a sequence of independent, standard univariate Brownian bridges, and the eigenvalue/eigenfunction pairs $(\lambda _{\ell },\phi _{\ell })$ satisfy

where the eigenfunctions $\phi _{\ell }$ are $r$-valued and form an orthonormal basis. Hence

In view of (ref), inference based on functionals of $ \Delta (u)$ requires an estimate of $\sigma _{0}^{2}$, and of the eigenvalues $\lambda _{\ell }$ in (ref). As far as the latter is concerned, note that $\mathbf{D}(t,s)$ is the long-run covariance of the sequence $\{X_{j}\}$. Therefore, a standard weighted-sum-of-covariances estimator can be employed for the consistent estimation of $\mathbf{D}(t,s)$, which in turn leads to estimates for the eigenvalues $\lambda _{\ell }$. We describe this procedure below. For a kernel function $\mathcal{K}:\mathbb{R}\rightarrow \mathbb{R}$ (see Assumption (ref) below), we define

where $h>0$ is a bandwidth parameter,

$\overline{X}_{j}\left( t\right) =X_{j}(t)-\widehat{\mu }_{N}(t)$, and $ \widehat{\mu }_{N}(t)=N^{-1}\sum_{j=1}^{N}X_{j}(t)$. (Above, we set $ \widehat{\boldsymbol{\gamma }}_{\ell }(t,s)\equiv 0$ for $\ell\geq N$). Note that, in ((ref)), we estimate the mean function $\widehat{\mu } _{N}(t)$ using the full sample. Under the null, this does not pose any problems given that $\mu (t)$ is constant. However, under the alternative $ \mu (t)$ is not estimated consistently; the bias in the estimation of $\mu (t)$ would enter $\widehat{\mathbf{D}}_{N}(t,s)$, making it diverge at a rate $h$. This is well-known in the literature on changepoint detection, and it has been associated with a decrease in power and the phenomenon known as \textquotedblleft non-monotonic\textquotedblright\ power (see casini2021prewhitened). This can be ameliorated by implementing a \textquotedblleft piecewise demeaning\textquotedblright , where the mean function is estimated by splitting the sample around each candidate changepoint $k$; however, unreported simulations show that using \textquotedblleft piecewise demeaning\textquotedblright\ yields some improvements in the power, but the test becomes (sometimes massively) oversized in small samples.

Assumptions (ref) and (ref) characterise the kernel $ \mathcal{K}\left( \cdot \right) $ and the bandwidth $h$, respectively; many of the customarily employed kernels satisfy Assumption (ref).

Lemma (ref) stipulates the consistency (in Frobenius norm) of $\widehat{ \mathbf{D}}_{N}(t,s)$. The lemma immediately entails that, for every fixed $ 1\leq \ell \leq N-1$,

where $\widehat{\lambda }_{1}\geq \widehat{\lambda }_{2}\geq \ldots $ are the eigenvalues of the operator $\phi \mapsto \int \widehat{\mathbf{D}} _{N}(t,s)\phi (s)ds$, $\phi \in L^{2}(\mathcal{T})$, suggesting that $ \widehat{\lambda }_{\ell }$\ is a good estimate of $\lambda _{\ell }$. Further, by the ergodic theorem (breiman:1968), under $H_{0}$

Hence, we can approximate the distribution of functionals of $\Delta (u)$ with functionals of

for sufficiently large $\widehat{M}$, $N$, using standard Monte Carlo techniques.

Testing for distributional change

We consider an extension of the testing procedure defined above to detect changes in the distribution of functional observations. Our approach is based on testing for the equality of the characteristic function, i.e., ultimately, on comparing expectations of a transformation of the data. Given that the data undergo a transformation, but the test statistics are the same, it can be expected that all the theory developed above can still be applied with no changes required. Indeed, compared with approaches based on using ((ref)) with $\eta <2$, our methodology has three distinct advantages. Firstly, the limiting distribution, in our case, involves the integral of the square of the Gaussian process (ref) , which greatly simplifies our computations. This is a consequence of having $ \eta=2 $; using $\eta <2$ would preclude this result ( biau:bleakley:mason:2016). Secondly, the binary segmentation algorithm discussed in Section (ref) can be applied also in this case, with no modifications required. This is a consequence of the fact that our test statistics for the detection of distributional changes are based on comparing expectations; conversely, as matteson:james:2014 put it, when using ((ref)) with $\eta <2$, binary segmentation \textquotedblleft cannot be applied in this general situation because it assumes that the expectation of the observed sequence consists of a piecewise linear function, making it only suitable for estimating changepoints resulting from breaks in expectation.\textquotedblright Thirdly, although we consider the empirical characteristic function, our tests can be immediately generalised to to include weighted empirical characteristic functions, or other transformations that may characterize the underlying distribution, such as e.g. moment generating function, or the Mellin transform, among other possibilities. This is a consequence of the fact that the theory in Section (ref) can be applied to test for the constancy of the expectation of any (univariate or multivariate) weakly dependent functional-valued time series, including transformations of functional-valued series.

Let $\mathbf{i}$ denote the imaginary unit, i.e. $\mathbf{i=}\sqrt{-1}$. Given a sequence of $Y_{\ell }=\{Y_{\ell }(s),0\leq s\leq 1\}$, $\ell =1,\ldots ,N$ of $L^{2}([0,1];\mathbb{R})$-valued functional observations, we consider the following null and alternative hypotheses

for $i=1,\ldots ,R$, again with the convention $k_{0}=1$ and $k_{R+1}=N$. Testing $H_{0}^{\prime }$ versus $H_{A}^{\prime }$ can be done with substantively weaker assumptions on the (moments of the) sequence $Y_{\ell }$ than what is required by Assumption (ref).

Inspired by berkes2009detecting, we pre-process the infinite dimensional data $Y_{\ell }(t)$ by projecting them into a finite dimensional vector

where $\{\psi _{\ell },\ell \geq 1\}$ is an orthonormal basis of $L^{2}[0,1]$ . Thence, we define the corresponding $\mathbb{C}$-valued random functions

where $t=(t_{1},\ldots ,t_{d})^{\top }\in \lbrack -1,1]^{d}$, and $d$ is user-chosen. Heuristically, $EX_{\ell }(t)$ is (an approximation of) the characteristic functional of $Y_{\ell }$, and therefore comparing averages of $X_{\ell }(t)$ before and after a point in time $k$ is a natural way of checking whether the distribution of $Y_{\ell }$ changes or not. Viewing each $\left\{ X_{\ell }(t),t\in \lbrack -1,1]^{d}\right\} $ in (ref) as an $\mathbb{R}^{2}$-valued random function $(\text{Re }\,X_{\ell }(t),\text{Im}\,X_{\ell }(t))^{\top }$, we may apply the test statistics proposed in Section (ref) to test the hypotheses $ H_{0}^{\prime }$ versus $H_{A}^{\prime }$.

Let $\xi _{j,\ell }^{(m)}=\int_{0}^{1}Y_{\ell }^{(m)}(s)\psi _{j}(s)ds$, and define $X_{\ell }^{(m)}(t)=\exp \left( \mathbf{i}\sum_{j=1}^{d}t_{j}\xi _{j,\ell }^{(m)}\right) $. We show that $\left\{ X_{\ell },-\infty <\ell <\infty \right\} $ is a Bernoulli shift sequence which satisfies Assumption (ref).

In order to construct the auxiliary functions $X_{\ell }(t)$ defined in (ref), one must first choose a basis $\{\psi _{\ell }\}$. Though any orthonormal basis of $L^{2}([0,1];\mathbb{R})$ will suffice, when the observations $Y_{\ell }$ satisfy $E|Y_{1}(s)|^{2}<\infty $, $s\in \lbrack 0,1]$, Principal Component Analysis (PCA) based approaches are among the most popular choices for selecting $\{\psi _{\ell }\}$, and typically lead to good finite-sample performance. Under the assumption that $ E|Y_{1}(s)|^{2}<\infty $, $s\in \lbrack 0,1]$, define

According to the PCA approach, the $\psi _{\ell }$ in ((ref)) are chosen as the eigenfunctions of $C(t,s)$

where $\chi _{1}>\chi _{2}>\ldots $, and $\{\psi _{\ell },\ell \geq 1\}$ are orthonormal - note the requirement that eigenvalues are well-separated, which is typical of (functional) PCA\ (see e.g. horvath:kokoszka:2012). With this choice of basis, typically the approximation $Y_{\ell }(t)\approx \sum_{i=1}^{d}\xi _{j,\ell }\psi _{\ell }(t)$ requires only a small number $d$ of projections for good finite-sample performance. We estimate the covariance function $C(t,s)$ in (ref) as

where $\widehat{\mu }_{Y,N}(t)=N^{-1}\sum_{j=1}^{N}Y_{j}(t)$ is the sample mean. If Assumption (ref) holds with $\beta >2$, then by the ergodic theorem it holds that

Thus, if $(\widehat{\chi }_{\ell },\widehat{\psi }_{\ell })$ are the eigenvalue-eigenfunction pairs defined by

where $\int \left\vert \widehat{\psi }_{\ell }\left( t\right) \right\vert ^{2}dt=1$, then for each fixed $1\leq \ell \leq N-1$, the eigenfunctions are estimated consistently modulo a sign - i.e., it holds that $E\Vert \widehat{ \psi }_{\ell }(t)-\iota _{\ell }\psi _{\ell }(t)\Vert ^{2}=o_{P}(1)$, where $ \iota _{\ell }$ is a random sign (see Theorem 2.8 in horvath:kokoszka:2012). Since the variables $X_{\ell }(t)$ do not depend on the sign of $\psi _{\ell }$, one can then construct $X_{\ell }(t)$ in (ref) based on

With the PCA-based choice (ref), it can be verified that $ \widehat{X}_{\ell }\left( t\right) =\exp \left( \mathbf{i}\sum_{j=1}^{d}t_{j} \widehat{\xi }_{j,\ell }\right) $ still satisfies Lemma (ref). Hence, all the results of Section (ref) hold when using the empirical energy function $V_{N}\left( k\right) $ based on $\widehat{X} _{\ell }\left( t\right) $ to test for $H_{0}^{\prime }$ in ((ref)) versus $H_{A}^{\prime }$ in ((ref)).

Simulations

We provide some Monte Carlo evidence on the performance of our test statistics, and some guidelines on how to implement the tests; further details and results (including a set of experiments on binary segmentation) are reported in Section (ref) in the Supplement. We use the following Data Generating Process (DGP), inspired by happ, based on a truncated multivariate KL representation

for $1\leq i\leq N$, where: $t\in \mathcal{T=}\left[ 0,1\right] $, $X_{i}(t)$ is univariate, $\mathcal{Z}_{\ell ,i}$ are $N\left( 0,1\right) $ and { uncorrelated across $\ell $, $\left\{ \phi _{\ell }(t),1\leq \ell \leq M\right\} $ form an orthonormal basis, and $\nu _{i}\left( t\right) $ is an i.i.d. Gaussian measurement error with mean zero and scale $ E\left\Vert \nu _{i}\left( t\right) \right\Vert ^{2}=\sigma _{\nu }^{2}$. As far as $\nu _{i}\left( t\right) $ is concerned, we consider two designs: a benchmark one with no measurement error (i.e., $\sigma _{\nu }^{2}=0$), and one with $\sigma _{\nu }^{2}=0.25$. We allow for serial dependence in the $ X_{i}(t)$'s through an AR($1$) structure in the $\mathcal{Z}_{\ell ,i}$ across $i$, viz. $\mathcal{Z}_{\ell ,i}=\rho \mathcal{Z}_{\ell ,i-1}+e_{\ell ,i}^{\mathcal{Z}}$ for all $1\leq \ell \leq M$, with $e_{\ell ,i}^{\mathcal{Z }}\sim i.i.d.N\left( 0,1\right) $ across $i$ and $\ell $. Under the null, we set $\mu _{i}(t)=0$ for all $1\leq i\leq N$, for simplicity and with no loss of generality. The eigenvalues $\lambda _{\ell }$ in ((ref)) are generated as

unreported simulations show that using different schemes (e.g. a linear, or a Wiener one) does not alter the results. The observations $X_{i}(t)$ are sampled on an equispaced grid of $S=128$ points. As is typical in FDA, a possible approach would be to pre-process and smooth the data, converting the discretely observed $X_{i}(t_{j})$, $1\leq j\leq S$ into functional objects by projecting them onto a suitably chosen basis; in our case, this would only help with dimensionality reduction, since the coefficients of the expansion are not required by any of our procedures. However, our test statistics are not particularly computationally demanding, and therefore pre-processing is not strictly required. Indeed, as hormann2022consistently put it \textquotedblleft for the processing of real data we will most often use the discretised curves anyway\textquotedblright . In our case, for example, the integral in equation ((ref)) will be computed numerically, and the most natural choice of nodes in the numerical computations are the discretised sampling points $t_{j}$, $1\leq j\leq S$ - hence, we suggest as a guideline that no data smoothing/pre-processing is carried out, at least for \textquotedblleft reasonable\textquotedblright\ values of $S$. As far as other specifications are concerned, we compute $\widehat{\sigma }_{N}^{2}$ and $\widehat{\mathbf{D}}_{N}(t,s)$ as described in Section (ref). We have used the Parzen kernel, and we have selected the bandwidth $h$ according to the optimal rules derived in andrews1991. We simulate $ \Delta _{N,M}(u)$ over a grid with exactly $N$ points, which we recommend in practical applications. All results are based on using an estimate $\widehat{ M}$ of $M$, chosen so that the first $\widehat{M}$ eigenvalues of $\widehat{ \mathbf{D}}_{N}(t,s)$ explain a prespecified amount of the total variability (we set this to $0.95$, which is a bit higher than in other papers, but still comes with a great dimensionality reduction). Critical values for weighted functionals of $\Delta _{N,M}(u)$ are computed using $500$ replications. All simulations are carried out with $1,000$ replications; all routines have been written using GAUSS 21.0.6. }

Empirical rejection frequencies under the null, at a nominal $5\%$ level, are reported in Table (ref); see also Section (ref) in the Supplement for further cases. In the i.i.d. case, our tests have excellent size control in all cases: the empirical rejection frequencies lie in the confidence interval $\left[ 0.036,0.064\right] $ even for sample sizes as small as $N=50$, and for all the values of $\alpha $ considered in our simulations. In general, our tests are almost never oversized, suggesting that spurious break detection is highly unlikely. When serial dependence (especially) and/or measurement errors are present, the tests are somewhat conservative for small samples and large $\alpha $, but this improves as $N$ increases. Upon closer inspection, this is due to the fact that the bandwidth $h$ employed in ((ref)) seems too high, and reducing it would increase the size; in turn, this suggests that, prior to implementing the tests, some qualitative considerations based on the presence of measurement error, and a bandwidth selection rule based on $ \alpha $, may yield improvements.

Turning to power, we compute empirical rejection frequencies under the at-most-one-change alternative where, in ((ref)), for simplicity in illustration, we shift the data by a constant after the breakpoint, namely:

where $\delta(t)\equiv C_\delta>0$, and the size of the change set to $ \|\delta\|=C_\delta \in \left\{ 0.4,0.8,1.2,1.6,2\right\} $. We consider two scenarios: a mid-sample break, with $k^{\ast }=\left\lfloor N/2\right\rfloor $ (Figure (ref)), and a late-occurring changepoint, with $ k^{\ast }=\left\lfloor 0.9N\right\rfloor $ (Figure (ref)). Results are obtained for $N=200$, and using $200$ replications to save computational time; for brevity, in Figures (ref) and (ref) we report only results for the i.i.d. case and the case of serial dependence and measurement error.\footnote{ Further results, with $N=100 $, confirm the findings reported in this section, and are reported in Section (ref) in the Supplement.} Figures (ref) and (ref) confirm that using higher $ \alpha $ is beneficial when the breakdate $k^{\ast }$ is close to sample endpoints, whereas, in the presence of mid-sample breaks, the test generally has good power, which tends to be lower as $\alpha $ increases (the discrepancy increases as $\left\Vert \delta \right\Vert $ declines).

In a second set of experiments, we explore the performance of our methodology to test for changes in the distribution proposed in Section (ref) via a small Monte Carlo exercise. We generate the one-dimensional functional data ${Y_{\ell }}(t)$\ using ((ref)) with no measurement error, viz. ${Y_{\ell }}(t)=\mu _{\ell }(t)+\epsilon _{\ell }(t)$ with $ \epsilon _{\ell }(t)$

We project $Y_{\ell }(t)$ onto its first Principal Component - that is, we use $d=1$ in ((ref)). We do this merely for computational simplicity; when computing the eigenvalues of the long-run variance matrix associated with $X_{\ell }(t)$, we use the algorithm in Section 3.2 in happ, based on the multivariate KL expansion. In Table (ref), we report the empirical rejection frequencies under the null, showing that our methodology has excellent size control for $ N\geq 100$; when $N=50$, tests appear to be mildly oversized.

We separately consider the following alternative hypotheses:

with $\left\Vert \mathcal{\delta }\right\Vert =1$, to consider changes in the mean function;

with $\mu _{\ell }(t)=0$ for all $1\leq \ell \leq N$, to consider a change in the (unconditional) variance which is helpful to understand whether our methodology can detect heteroskedasticity; and lastly

where $t_{m,\ell }^{\left( 3\right) }$\ are i.i.d. random variables, independent across $m$ and $\ell $, with a Student's t distribution with $3$ degrees of freedom, and $\mu _{\ell }(t)=0$ for all $ 1\leq \ell \leq N$, as a more general alternative where the data, after a period of \textquotedblleft normal\textquotedblright\ fluctuations, exhibit heavy tails. Results in Table (ref) show that our tests - even when using $d=1$ - have excellent power under all cases in the presence of a mid-sample break, which is also estimated correctly. Hence, the test developed in Section (ref) can be used to detect shifts in the mean, in the variance, or in the tails - of course, the test is an omnibus test, and therefore it is non-constructive in that, upon rejecting the null, it does not indicate a specific alternative. In the case of end-of-sample breaks, Table (ref) shows that the test is sensitive, as expected, to the choice of $\alpha $, and that as $\alpha $ approaches $1$ the power increases, as does the accuracy in estimating the changepoint. The test performs very well, even for small sample sizes ($N=100$) in the presence of changes in the mean and in the tails (i.e., under ((ref)) and ((ref)) respectively), whereas its performance is less good in the presence of shifts in the variance (i.e., under ((ref))) although it picks up as both $\alpha $ and $N$ increase.

Changepoint detection in high-frequency financial data

We apply our tests for changepoint detection in the mean and in the distribution of intraday return patterns (on a month-on-month basis) in high-frequency trading of the S&P 500 index.\footnote{ In Section (ref) in the Supplement, we study temperature data, which is another classical application of FDA (see e.g. berkes2009detecting).} High-frequency trading data lend themselves to being studied through the lenses of FDA, as they typically contain a huge amount of data for which a parsimonious representation is necessary; furthermore prices change continuously on a daily basis, and therefore daily prices are genuinely functional objects, whose sampling points are the observed prices recorded over several points in time each day. Examples of applications of FDA to high frequency financial data include e.g. genccay2001introduction, and muller2011functional; kokoszka2012functional consider an alternative definition of return (know as cumulative intraday returns, or CIDRs) which seems to be particularly suited for predictions using FDA.\footnote{ In Section (ref) in the Supplement, we complement our analysis by considering changes in the mean and in the distribution of CIDRs.}

Returns are calculated from closing prices that are recorded at equispaced 5-min intervals between $00:00$ and $23:55$ each day, corresponding to a sampling frequency $S=276$. We have used the period spanning from January 3rd, $2022$, until September 3rd, $2023$. In order to balance the sample and ensure that there are $S=276$ sampling points each day, we have removed $21$ trading daily curves in which some data were missing,\footnote{ A list of the relevant days is available upon request.} for a total of $414$ functional datapoints. Denoting prices at day $i$ as $P_{i}\left( t\right) $ , we construct month-on-month log returns as $Y_{i}\left( t\right) =\ln P_{i}\left( t\right) -\ln P_{i-21}\left( t\right) $, having used $21$ lags as the average amount of trading days in a month; hence, the resulting sample size is $N=393$, effectively starting from February 3rd, $2022$. We have implemented our test using the guidelines and specifications suggested in Section (ref).\footnote{ Critical values for weighted functionals of $\Delta _{N,M}(u)$ are computed using $500$ replications, due to the reduced computational times in the empirical exercise; we note however that results do not differ in any significant way upon altering this specification.} In particular, when using the KL expansion, we employ a number of bases $\widehat{M}$, chosen so that the first $\widehat{M}$ eigenvalues of the estimated long-run variance explain $95\%$ of the total variability; we note that, in unreported experiments, altering this specification did not change any of the final results. Contrary to Section (ref), we use an estimate of the covariance kernel\ $\mathbf{D}\left( t,s\right) $ computed using the estimator of the covariance functions ((ref)) with pre- and post-break demeaning; we note that we tried to use $\widehat{\mathbf{D}} _{N}(t,s)$ defined in ((ref)) - that is, without demeaning before and after the candidate breakdate - but results do not change in any way. \footnote{ As far as other specifications of $\widehat{\mathbf{D}}_{N}(t,s)$ are concerned, we have used a Parzen kernel and the bandwidth $h$ chosen according to andrews1991.} Upon inspection, this is due to the fact that, across all exercises, the bandwidth is selected as at most $h=1$; this, in turn, suggests that the data have only little serial dependence. We tried to assess the sensitivity of our results by varying $h$, but virtually no changes were noted. Finally, results are reported at a nominal level of $ 5\%$; however, we have also carried out - by way of sensitivity analysis - detection at $1\%$ and $10\%$ levels, and results are discussed in the notes of our tables.

We begin by applying our test for a changepoint in the mean. In Table (ref), we report results for $\alpha =0.5$, but we also tried $ \alpha =0$ and $\alpha =0.99$, obtaining exactly the same outcome: there is only one break, located at October 20th, $2022$, whose estimate appears to be remarkably robust (note also the discrepancy between the daily averages). Whilst it is difficult to associate a particular event to that date, in general the common wisdom among financial analysts is that S&P 500 index started recovering, after a turbulent year and after hitting its low in October $2022$, around the second half of that month.\footnote{ A qualitative description can be e.g. found at https://www.usbank.com/investing/financial-perspectives/market-news/is-a-market-correction-coming.html } No further breaks in the mean function were found.

It is well known (see e.g. kim2004more) that the behaviour of financial markets is characterised by not being adequately described by the Gaussian distribution; hence, it is important to check if there are changes not merely in the mean (or in the variance), but in the whole distribution. Thus, after finding the presence of a break in the mean, we demean the data in each of the two segments around October 20th, $2022$, and carry out the test for distributional changes, on the demeaned data, discussed in Section (ref). We use exactly the same specifications as in Section (ref), using only one principal component (i.e. $d=1$) in the construction of $X_{i}(t)$ in ((ref)). Tests are carried out at a nominal level of $5\%$ by default (we also tried $1\%$ and $10\%$, see the notes to Table (ref)); when using binary segmentation, we use, as threshold, $\tau _{N}=c_{a}\sqrt{\ln N}$, where $a$ is the nominal level of the test; results are generally robust to this (we tried $\tau _{N}=c_{a}\ln \ln N$, and $\tau _{N}=c_{a}{\ln N}$, and no changes were noted). We report our findings using weights $\alpha =0$, $0.5$ and $0.99$; we used binary segmentation, so the case $\alpha =0$ is not reliable per se, as it may lead to overestimation of the number of regimes, but we use it as a benchmark for the other two sets of results.

The results in Table (ref) show a much richer picture that changes in the mean alone. Interestingly, the same changepoints are found across all values of $\alpha $, including $\alpha =0$; the only difference is in the date of the break estimated between July and October $2022$, which appears to be estimated $8$ trading days later when using larger values of $ \alpha $. Otherwise, results are exactly the same; indeed, we also experimented with other values, but results were the same even in those cases. The same robustness was found when altering other specifications of the procedure, e.g. the estimation the covariance kernel or of the number of terms $\widehat{M}$ in the KL expansion. The estimated breakdates are, at least in some cases, highly suggestive; interestingly, in all cases, a change in regime corresponds to a change in the sign of our measure of skewness, which confirms the stylised fact that skewness is time-varying ( alles1994regularities; bekaert1998distributional). The first break, recorded at March 20th, $2022$, corresponds to a peak in the S&P 500, after which the market entered a bear phase to stay below that peak until July $2023$; after removing the first 21 observations, our month-on-month return series starts effectively in February, so that the first regime - characterised by a strongly negative measure of skewness - reflects the uncertainty due to the war in Ukraine and its impact on the global economy. The second regime, between March 20th and April 20th, is characterised by a positive skewness, possibly indicating that the market - after the stalling of the Russian offensive - was expecting an upward price movement. This did not materialise, and in April $2022$ the market experienced a strong correction, partly also due to inflation expectation and underperformance of high-tech firms.\footnote{ https://www.marketwatch.com/story/the-stock-market-swoon-just-sent-the-s-p-500-into-its-second-correction-of-2022-11651265882 } After April 20th the market entered a bear phase, chracterised by negative skewness, until a turning point was reached on July13th, $2022$ on account of the FED ending (temporarily) its rate hiking. A correction occurred after August 30th, $2022$, with a slump that lasted until approximately the second half of October (the changepoint was recorded on October 20th, $2022$). From thereon, the market started a rebound which lasted for the remainder of our sample period; during this long horizon, the market was again characterised by positive skewness.

Discussion and conclusions

In this paper, we propose a family of weighted statistics to detect changepoints possibly dependent, multivariate functional data. Although we focus our exposition on the well-studied case of changes in the mean, our tests can be applied to much more general changepoint problems, such as detecting changes in the whole distribution. We base our test statistics on the notion of energy distance, a recently proposed measure of proximity between distributions; we use a version of the (empirical) energy distance which is particularly suited to determining the equality of the first moment of random variables, showing that, under the null of no breaks, this is related to the familiar CUSUM\ process. Our statistics can be applied under very general forms of (weak) serial dependence, thus being suitable for the analysis of several datasets, including meteorological, financial and economic time series. By using a set of weights which place more emphasis on observations occurring close to the sample endpoints, we are able to detect changepoints occurring very close to the beginning/end of the sample. Also, our approach is sufficiently flexible to allow for generalisations to e.g. testing for changepoints in the (marginal) distributions of a sequence. In particular, our approach is based on checking whether expectations of functions of our data remain constant over time; consequently, we can use all the technology available in the literature, such as e.g. binary segmentation in order to detect (and estimate the number and location of) multiple changepoints. An important feature of our procedures is its computational simplicity: critical values can be derived with arbitrary precision, and this requires only the eigenvalues of the covariance operator of the data, which can be quickly computed via any available statistical package. Our simulations show that our statistics have excellent finite sample performance even for small samples, thus making their use possible in virtually all contexts involving FDA.

This work leads to several possible future directions, including extensions to energy distances for functional data beyond the case $\eta =2$, and more broadly further exploration of generalized energy distances. The use of the characteristic function in Section (ref) can be viewed as a finite-dimensional approximation of the characteristic functional; an interesting direction would be to more deeply explore finite-dimensional approximations of the characteristic functional and similar transformations in the context of functional time series.

\setcounter{section}{0} \setcounter{subsection}{-1} \setcounter{subsubsection}{-1} \setcounter{equation}{0} \setcounter{lemma}{0} \setcounter{theorem}{0}

Further Monte Carlo evidence and guidelines

Empirical rejection frequencies under the null: further results

We complement the results in Table (ref) by considering the cases of i.i.d. data with measurement error, and the case of serially dependent data without measurement error.

Empirical rejection frequencies under the alternative: further results

We begin by reporting the power against one changepoint, with the same design as in equation ((ref)) using $N=100$. As can be seen in Figures (ref) and (ref), the results are similar, although the test is less powerful compared to the results in Figures (ref) and (ref), which is expected due to the smaller value of $N$. In particular, in the presence of an end-of-sample break, power is ensured only for large values of $\alpha $.

We now report the median values of the estimated breakdate $k^{\ast }$, in the case of a single changepoint (when this is detected), under the same set-up as in Section (ref) - see equation ((ref)) in particular. Results in Tables (ref)-(ref) should be read in conjunction with Figures (ref)-(ref) and (ref)-(ref), and broadly confirm the theory spelled out in Theorem (ref). In the case of mid-sample breaks, the estimator of $k^{\ast }$ is usually very good when a changepoint is detected (Tables (ref)-(ref)), even for small break sizes like $\left\Vert \mathcal{\delta }\right\Vert =0.4$; this is true across all values of $\alpha $, although, in the case of small breaks ($\left\Vert \mathcal{\delta }\right\Vert =0.4$), the performance of $ \widehat{k}$ when $\alpha $ gets closer to $1$ seems to worsen. Conversely, when breaks occur close to the end of the sample ($k^{\ast }=0.9N$), results in Table (ref) and (ref) differ dramatically across $\alpha $: as expected, when $\alpha $ increases, $\widehat{k}$ performs better, and it performs very well when $\left\Vert \mathcal{\delta } \right\Vert \geq 0.8$ (and even more so when $N=200$). Interestingly, in this case $\widehat{k}$ appears to have a downward bias, which vanishes as $ N\left\Vert \mathcal{\delta }\right\Vert $ increases.

Binary segmentation: pesudocode and Monte Carlo evidence

We begin by reporting some pseudocode for the practical implementation of the algorithm. Let, for short

The pseudocode is in Algorithm (ref) below.

We now report a small Monte Carlo exercise to assess the performance of the binary segmentation procedure discussed in Section (ref). In particular, we consider the following DGP

where the random part $\sum_{\ell =1}^{M}\lambda _{\ell }^{1/2}\mathcal{Z} _{\ell ,i}\phi _{\ell }(t)+\nu _{i}\left( t\right) $ is generated in the same way as in Section (ref), and we consider two mid-sample changepoints (i.e., $R=2$) in the mean function $\mu _{j}(t)$, viz.

with $k_{1}^{\ast }=\left\lfloor 0.35N\right\rfloor $ and $k_{2}^{\ast }=\left\lfloor 0.7N\right\rfloor $, and $\delta_i(t)\equiv C_{\delta_i}$ are constants, $C_{\delta_1}=0$, $C_{\delta_2}=2$ and $C_{ \mathcal{\delta } _{3}} =3$. When using Algorithm (ref), we select the threshold

where $c_{a}$ is the critical value at nominal level $a=0.05$.

Results in Table (ref) contain measures of location of the estimator of the number of changepoints $\widehat{R}$, and the average values, across simulations, of the estimated breakdates, using $500$ simulations. Results are obtained for $N=200$ and with $\mathcal{Z}_{\ell ,i} $ generated as i.i.d. across $1\leq i\leq N$; we consider the presence of measurement errors, setting $\sigma _{\nu }^{2}=0.25$ as in Section (ref), but in unreported experiments we noted that having $\sigma _{\nu }^{2}=0$ does not change the results in any significant way.

In addition to investigating the performance of binary segmentation in the presence of shifts in the mean as per ((ref)), we also explore its performance in the presence of changes in the distribution. In particular, we consider the \textquotedblleft epidemic\textquotedblright\ alternative in a model with zero mean

and

where, as in ((ref)), $t_{m,\ell }^{\left( 3\right) }$\ are i.i.d. random variables, independent across $m$ and $\ell $, with a Student's t distribution with $3$ degrees of freedom; all the other specifications are the same as above. We use $N=200$ and, as above, $ k_{1}^{\ast }=\left\lfloor 0.35N\right\rfloor $ and $k_{2}^{\ast }=\left\lfloor 0.7N\right\rfloor $. Alternative ((ref)) represents a case, relevant in practice, where the data experience a period of turbulence characterised by heavy tails, after which they revert to normal. Results are in Table (ref); we found $\tau _{N}=c_{\alpha }(\ln N)^{1/2}$ to be a better choice in this case, and we suggest this choice of threshold when testing for changes in the distribution.

\setcounter{subsection}{-1} \setcounter{subsubsection}{-1} \setcounter{equation}{0} \setcounter{lemma}{0} \setcounter{theorem}{0}

Further empirical evidence

Further empirical evidence: changepoint analysis of cumulative intraday returns

We complement our findings in Section (ref) by applying our tests for changes in the mean and in the distribution to cumulative intraday returns (CIDRs henceforth), whose usefulness is demonstrated in a contribution by kokoszka2012functional. We use the same dataset as in Section (ref), having removed the same curves consisting of partial trading days. CIDRs are defined as

where $P_{i}\left( t\right) $ is the daily price evaluated at $t$, and $ t_{0} $ is, for each day $1\leq i\leq N$, the start of trading for the day (in our case, midnight). Contrary to the use of month-on-month returns, in this case we can use the whole sample of daily curves consisting of $N=414$ functional datapoints.

Tests have been applied with the same specifications as in Section (ref) in the main paper. We did not find any changepoints in the mean, irrespective of the value of $\alpha $. Conversely, applying the test for distributional changes to the demeaned data, several changepoints are found, summarised in Table XXX. In the table, as in the rest of the paper, we have computed the measures of skewness and kurtosis as

Further empirical evidence: changepoint analysis of temperature data

In this section, we illustrate our approach using temperature data, where FDA is applied \textquotedblleft naturally\textquotedblright ; more broadly speaking, recent contributions in the area of climate science show that using time series methods can be beneficial (see diebold2022probability; diebold2023will; and ditlevsen2023warning).

Following berkes2009detecting, we use a sample of $N=251$ yearly curves, recorded on a daily basis between $1772$ and $2022$. Each curve contains average daily temperatures (in degrees Celsius) recorded in Central England; apart from removing the data corresponding to February 29th in leap years in order to balance the sample, no further transformation is applied to the data.\footnote{ The data have been downloaded from \url{https://www.metoffice.gov.uk/hadobs/hadcet/}, where a brief description of the dataset can also be found. A more complete description of the data can be found in parker1992new, to which we refer for details.} Our techniques are particularly suited to this dataset for a number of reasons: firstly, we do not need to invert any large-scale matrix, contrary to horvath1999testing, and therefore we can use the daily sampling frequency as opposed to transforming it into monthly averages; secondly, temperature data might exhibit linear or nonlinear serial dependence (see e.g. bowers), which our tests are designed to take into account, unlike those proposed in berkes2009detecting; and, finally, our weighted test statistics are also suited to detect changepoints occurring close to the end of the sample, thus allowing to shed light on the presence and extent of changes in average temperatures in recent years. We apply our tests for changepoints in the mean using our tests with $\alpha \in \left\{ 0,0.5,0.65\right\} $, by way of comparison and robustness check; we note that using different values of $\alpha $ does not alter the main conclusions, although higher values of $\alpha $ seem to estimate the changepoint date later and later. In order to take into account the possible presence of multiple changes, we apply binary segmentation. We have implemented our test using the same specifications as described in Section (ref) in the main paper, also carrying out the same robustness checks with no noticeable changes in the results.

Results are in Table (ref); in Figure (ref), we also report the average temperature functions between each of the estimated changepoints for the various values of $\alpha $. With small and medium values of $\alpha $ (i.e., $\alpha =0$ and $\alpha =0.5$), we identify three changepoints. The first one to be identified (corresponding to the \textquotedblleft strongest\textquotedblright\ break) is estimated to have occurred in $1919$. This result is essentially in agreement with the findings in berkes2009detecting (and also in horvath1999testing), where the first changepoint is found around $ 1926$. This estimated date corresponds to the so-called Early Twentieth Century Warming (hegerl2018early), a well-documented phenomenon which partly coincides with the well-known phenomenon of warming of the Arctic (bengtsson2004early), and which \textquotedblleft still defies full explanation\textquotedblright\ (bronnimann2009early, p. 735). Our estimated date is ealier than that of berkes2009detecting, which could be ascribed to the estimation error, but also to the \textquotedblleft pull\textquotedblright\ effect, on the data, of the UK heatwave of $1911$. We also estimate a breakdate at $ 1842 $, which could be ascribed to the anthropogenic effect of the Industrial Revolution, and again it is similar to the estimate of $1850$ in berkes2009detecting. On the other hand, berkes2009detecting also find one changepoint in $1808$. None of our statistics finds evidence of a changepoint around this date, even at $10\%$ nominal level; on account of the lack of serial correlation, our data could be roughly interpreted as falling into the \textquotedblleft i.i.d. with measurement error\textquotedblright\ category, for which our simulations indicate no undersizement and excellent power even in sample sizes. Finally, we find clear evidence of a changepoint in $1987$ ($1988$ when using $\alpha =0.5$, and $1989$ when using $\alpha =0.65$). This is significant even at $1\%$ nominal level, which corresponds to the beginning of (rapid) global warming - the late $1980$'s date confirms the statement, in the \textquotedblleft State of the World 1989\textquotedblright\ WorldWatch report ( koutaissoff1989state), that the 90's would be the \textquotedblleft turnaround decade\textquotedblright\ as far as climate change is concerned. Indeed, Figure (ref) shows very clearly the presence of a pronounced increase in average daily temperatures between the first and the fourth subsamples. When using $\alpha =0.65$, essentially the same results are found, but the first changepoint is estimated at $1932$ , i.e. one decade later than the other two changepoint estimates. This can be read in the light of the results in Table (ref), which suggest that, as $\alpha $ increases, the estimated changepoint may be increasingly biased.

Finally, we also considered the possible presence of changes in the distribution of temperature data, applying the test developed in Section (ref). To this end, we demeaned the data in each segment, and applied the test in Section (ref) using $d=1$. Results are in Table (ref). No changepoints were detected for $\alpha =0$ and $\alpha =0.5$, even at $1\%$ nominal level, suggesting that, if a changepoint is present, this is located towards the sample endpoints; indeed, when using $ \alpha =0.65$, the presence of one changepoint emerges even at $10\%$ nominal level, with estimated date $1996$. Comparing descriptive statistics, skewness and kurtosis in the two subperiods seem very similar (in both cases suggesting Gaussianity); conversely, the estimated variances seem to indicate that there is a changepoint in the variability of temperatures between the two subperiods.

Descriptive statistics for the whole sample, and for the segments identified when using $\alpha =0$, are in Table (ref). Results using other values of $\alpha $ (and, therefore, other estimated changepoints) are available upon request.

\setcounter{subsection}{-1} \setcounter{subsubsection}{-1} \setcounter{equation}{0} \setcounter{lemma}{0} \setcounter{theorem}{0}

Preliminary lemmas

Henceforth, unless stated otherwise, we carry out our proofs for the case $ r=d=1$, for simplicity and without loss of generality. We use the following notation: $\left\lceil \cdot \right\rceil $ is the ceiling function, that is the function that rounds a number to the nearest, largest integer; $C$ denotes a generic constant independent of $N$, $k$, $t$ that may change from line to line.

We begin by recalling some results in berkes:horvath:rice:2013.

For convenience, throughout the remainder of this section we use the notation $S_{k}(t)=\sum_{i=1}^{k}X_{i}(t),\quad t\in \mathcal{T}$, and without loss of generality under $H_{0}$ we assume $\mu _{i}=0$ for all $i$, so that $EX_{i}(t)=E\epsilon _{i}(t)=0$. We now establish some intermediate weighted approximation results related to the partial sums $S_{k}$.

The next lemma provides an approximation for $V_{N}(k)$ in terms of $S_{k}$. For each $1\leq k\leq N$, recall $u=k/N$ and define